Interaction between SNCA, LRRK2 and GAK increases susceptibility to Parkinson's disease in a Chinese population

PD is a complex disease, and may result from gene–gene and gene–environment interactions. There are limited studies on gene–gene interactions in PD. We and others have previously shown that SNCA rs356219, LRRK2 (rs2046932 and rs7304279) and GAK (rs1564282) are risk factors in sporadic PD. Since the expression of SNCA and neurotoxicity of alpha-synuclein are affected by LRRK2 and GAK, we hypothesize that their genetic risk variants may interact with each other. Here we investigated the interaction of SNCA rs356219, LRRK2rs7304279 and rs2046932 and GAK rs1564282 using the Multifactor Dimensionality Reduction (MDR) in a Chinese PD patient–control series (534 patients and 435 controls) and the cumulative risk effect of SNCA, LRRK2 and GAK. The MDR analysis showed a significant gene–gene interaction between the rs356219 of SNCA, rs2046932 of LRRK2 and rs1564282 of GAK. Moreover, individuals with increasing numbers of variants had an increasing likelihood of having PD, compared with those carrying none of the variants. The estimated OR for developing PD in individuals carrying 3 variants was 5.89. We demonstrated for the first time that SNPs in SNCA, LRRK2 and GAK interacted with each other to confer an increased risk of PD. In addition, PD risk increased cumulatively with the increasing number of variants.

Introduction
Parkinson's disease (PD; OMIM #168600) is one of the most common neurodegenerative diseases, affecting about 2% of people over the age of 65 [1]. The disease is characterized by resting tremor, rigidity, bradykinesia and postural instability, and associated with selective loss of dopamine neurons (DA) and formation of Lewy bodies (LBs) [2]. About 90% of cases are apparently sporadic, monogenetic mutations are now estimated to cause about 10% of PD cases [1]. Over the past two decades, numerous studies have been conducted to explore the genetic basis of PD. A recent large-scale meta-analysis of genomewide association data in Europe has replicated 22 loci for Parkinson's disease (including SNCA, MAPT, LRRK2, BST1, GAK), while it has also reported 6 new risk loci (SIPA1L2, INPP5F, MIR4697, GCH1, VPS13C and DDRGK1) [3]. However, most reports just analyzed the association between single nucleotide polymorphism (SNP) and sporadic PD with small attributable risk. To date, few studies have evaluated gene-gene and gene-environment interactions.
SNCA is one of the causative genes of PD associated with alphasynuclein in Lewy bodies [4,5]. Alpha-synuclein, encoded by SNCA, is the major fibrillar component of Lewy bodies. It is a presynaptic phosphoprotein with an intrinsic propensity to aggregate and has been demonstrated to have a role in both inherited and idiopathic PD [6]. Based on the critical role that alpha-synuclein may play in PD and functional connections with other identified genes, some efforts have been made in exploring the gene-gene interactions focusing on SNCA [7][8][9][10].
Besides alpha-synuclein, the LRRK2 proteins have also been identified in the Lewy bodies by immunohistochemical studies [11]. Mutations in LRRK2 are generally considered the most common genetic determinant of familial and sporadic PD. A recent study reported a strong interaction between LRRK2 and alpha-synuclein at both over-expressed and endogenous levels [12]. Moreover, LRRK2, which is thought to be an upstream factor in the neurodegenerative pathway, have also been shown to enhance alpha-synuclein-mediated cytotoxicity [13]. MAPT and SNCA may modify LRRK2-related risk for PD [14].
Another newly identified gene, named GAK, was one of the 137 genes differentially expressed in the substantia nigra pars compacta of PD patients when compared with controls, with a 1.56-fold change [15]. GAK was recently reported to have a modification on alphasynuclein expression and toxicity [16]. In addition, GAK and two other proteins have been identified as binding partners of LRRK2 and these eNeurologicalSci 1 (2015) [3][4][5][6] proteins form a complex that promotes clearance of Golgi-derived vesicles through the autophagy-lysosome system both in vitro and in vivo [17]. Based on these biological evidences, we hypothesize that the genetic risk variants of these genes may interact with each other.
In 2001, Ritchie et al. first developed a multifactor dimensionality reduction (MDR) method for detecting and characterizing high-order gene-gene and gene-environment interactions in case-control and discordant-sib-pair studies with relatively small samples [18]. The MDR method is nonparametric (i.e., no hypothesis about the value of a statistical parameter is made), is model-free (i.e., it assumes no particular inheritance model), and is directly applicable to case-control and discordant-sib-pair studies 3. So far, this new method has been successfully used for genetic studies of common complex multifactorial diseases [19].
In our previous studies, we have demonstrated that SNCA rs356219, LRRK2 (rs7304279 and rs2046932) and GAK rs1564282 increased the risk of sporadic PD in a Han population from mainland China [20][21][22].
Here we examined the interaction using the MDR method and cumulative risk effect of SNCA, LRRK2 and GAK. To our knowledge, this is the first attempt to explore the interactions between SNCA, LRRK2 and GAK.

Subjects
A total of 969 ethnic Han Chinese study subjects comprising 534 independent sporadic PD patients and 435 neurologically healthy control individuals were recruited from the Department of Neurology at West China Hospital. All patients (the mean age 58.30 ± 11.05, range from 30 to 86, 41.1% women) were examined and evaluated longitudinally by two movement disorders neurologists and diagnosed with idiopathic PD according to the UK Parkinson's Disease Society Brain Bank Clinical Diagnostic Criteria [23]. None of the patients had a positive family history. All control individuals (the mean age 52.30 ± 14.11, range from 30 to 91, 43.9% women) were healthy volunteers without any neurological or psychiatric diseases and were recruited from the same ethnic group. Written informed consent was obtained from all participants. DNA was extracted from blood leukocytes by standard procedures. This study was approved by the Ethics Committee of Sichuan University.

Genetic analysis
All participants were genotyped by using a matrix-assisted laser desorption/ionization time-of-flight (MALDI-TOF) mass spectrometry on a MassArray system (Sequenom, San Diego, USA). Approximately 15 ng of genomic DNA was used to genotype each sample. Locusspecific polymerase chain reaction (PCR) and detection primers were designed using the MassArray Assay Design 3.0 software (Sequenom, San Diego, USA). The sample DNAs were amplified by primers flanking the targeted sequence, followed by dephosphorylation and allele-specific primer extension. Cleaned extension products were loaded into a 384format Spectro-Chip, and subjected to MALDI-TOF mass spectrometry. The resultant data were analyzed by the Sequenom MassArray Typer software (Sequenom, San Diego, USA). The methods were carried out in accordance with the approved guidelines.

Statistical analysis
We assessed Hardy-Weinberg equilibrium (HWE) in cases and controls with a Fisher's exact test. The frequencies of the alleles and genotypes in the patients and control groups were analyzed using the Chi-square test. The clinical data were analyzed by student T test or Mann-Whitney U test. A two-tailed P-value ≤ 0.05 was considered statistically significant. Statistical analysis was performed using the Statistical Package for the Social Sciences version 16.0 (SPSS, Chicago, IL, USA) for Windows.
We used the MDR to analyze the gene-gene interactions between SNCA, LRRK2 and GAK. MDR is a powerful method for analyzing interactions. It pools genotypes into "high-risk" and "low-risk" groups to reduce multidimensional data into one-dimension [18]. The specific methods and procedures of MDR are introduced by Ritchie et al. [18]. The classification errors and the prediction errors are estimated by 10-fold crossvalidation. Additionally, 10-fold cross-validation consistency (CVC) and permutation test are used to select the final best model.
Chi-square test was used to compare the clinical manifestations with different numbers of associated variants.

Characteristics of participants
Data from a total of 969 individuals including 534 sporadic PD patients and 435 healthy controls were analyzed. All polymorphisms were in Hardy-Weinberg equilibrium for PD patients and controls. The mean age of PD patients, which consisted 41.1% women, is 58.30 ± 11.05 years, ranging from 20 to 80. The mean age of the healthy controls is 52.30 ± 14.11 years, a little younger than the PD groups.

Associations between each single nucleotide polymorphism with PD susceptibility
Some of the data for each single nucleotide polymorphism (rs356219, rs7304279, rs2046932 and rs1564282) have been reported in our previous papers which showed that these are risk variants [20][21][22]. Here, we included 969 subjects who have full genotype information of the four SNPS, to do the gene-gene interaction analysis. Subjects with GG + GA genotypes of rs356219 have an increased risk compared to those with AA genotype. So did the TT + CT genotypes of rs7304279, TT + TC genotypes of rs2046932 and TT + TC genotypes of rs1546282 (Table A.1). 3.3. MDR analysis of gene-gene interactions between SNCA (rs356219), LRRK2 (rs7304279, rs2046932) and GAK (rs1564282) Gene-gene interactions were investigated for PD using MDR method, and two significant interactions were found (Table A.2). Compared with the other models, a two-locus model incorporating rs356219 and rs2046932 was the best with the maximum testing accuracy of 0.6091 and cross-validation consistency of 10 out of 10. However, the best three-locus model including rs356219, rs2046932 and rs1564282 also had a maximum cross-validation consistency of 10 and higher testing accuracy of 0.6059, secondary to the best two-locus model. The best two-and three-locus models were both significant at P b 0.05. Though subtle difference existed, to pinpoint the polymorphisms that are of particular interest, the three-locus model was regarded as the overall best MDR model in our study.

Clinical characteristics of PD carrying different risk variants
The clinical features (gender, age, age at onset, disease duration, onset symptoms, and Hoehn-Yahr stage) were similar between PD patients carrying none of the variants and PD patients carrying one to three variants (Table A.4).

Discussion
Here we found that the four SNPs (rs356219, rs2046932, rs7304279 and rs1564282) of three genes (SNCA, LRRK2 and GAK) are associated with PD susceptibility. Utilizing MDR analysis, we demonstrated for the first time that compared with those carrying none of the variants, the estimated OR for developing PD in individuals carrying 3 variants (rs356219 of SNCA, rs2046932 of LRRK2 and rs1564282 of GAK) was 5.89 (OR = 5.89, 95% CI = 2.08, 16.68, P b 0.000). This suggests that these variants have a cumulative risk effect on increased susceptibility to PD.
SNCA gene encodes alpha-synuclein a major component of Lewy bodies. Lewy body contains not only alpha-synuclein, but also MAPT and LRRK2 proteins [11,24]. MAPT induces fibrillization of alphasynuclein, such that fibrillization is considered to be an initial step in the generation of PD-pathogenic alpha-synuclein aggregates [25]. In 2007, Goris et al. observed a significant interaction between the MAPT H1 haplotype and SNCA 30 SNP variants, with the combination of both risk genotypes doubling the risk for PD [4]. However, Mamah et al. failed to find any joint effects for MAPT H1 haplotype and SNCA REP1 variants with PD susceptibility in their study two years before [6]. Another study including Italians again did not find interaction between the two genes [7]. Wider et al. failed to discover the gene-gene interactions between SNCA and MAPT genes in their SNCA, MAPT and GSK3β association study, which was based on the phosphorylation of GSK3β on tau, encoded by MAPT [8].
Alpha-synuclein expression may be influenced by LRRK2 [26][27][28]. Whereas tau is a microtubule-associated protein, LRRK2 is known to interact with microtubules [29]. MAPT is functionally linked to both LRRK2 and alpha-synuclein with implications for the pathogenesis of PD. One study revealed no interactions between SNCA, MAPT and LRRK2 [30], while another study showed that MAPT and SNCA may influence LRRK2-related risk for PD in a Chinese cohort [14]. Others found that MAPT and LRRK2 genes interacted to increase the susceptibility to PD in Chinese individuals [31], but this was not replicated in another study [32].
Cyclin G-associated kinase (GAK) is a 160 kDa serine/threonine kinase, located on 4p16.3 [33]. It is a ubiquitously expressed protein, containing highly conserved serine/threonine kinase, PTEN and J-domains [34]. The clathrin binding C-terminal domain of GAK has been shown to bind pre-cathepsin (CTSD), which was implicated as the main lysosomal enzyme involved in alpha-synuclein degeneration [35]. In addition, Dumitriu et al. showed that the T allele of the rs1564282 GAK SNP was nominally associated with decreased expression of GAK exon 28 in frontal cortex, and the decreased expression of GAK might be a possible mechanism by which this SNP is linked to increased risk for PD. They found that GAK changes alpha-synuclein expression levels and toxicity, suggesting that the GAK protein and a-synuclein interact in a pathway involved in PD pathogenesis [13].

Conclusions
In conclusion, we demonstrated for the first time that SNPs in SNCA, LRRK2 and GAK interacted with each other to confer an increased risk of PD. In addition, PD risk increased cumulatively with the increasing number of variants. Further studies to explore the biological basis of the gene interactions will be interesting.

Author contributions
W. Y. and L. C. designed the research, collected the patients, genotyped the SNPs and wrote the paper. N. L. and L. W. collected the patients. E. T. and R. P. designed the research. All authors reviewed the manuscript.

Conflict of interest
None.

Acknowledgments
The study was supported by Duke-NUS Graduate Medical School, Singapore Millennium Foundation, National Medical Research Council, and West China Hospital of Sichuan University. We gratefully acknowledge Professor Dong Zhou, Li He, Guanggu Yuan, Yingru Gou, Tao Li and Qiang Wang. We would also thank Dr. Wenjun Chen, Yan Wu, Xingkai An, ZijuanZhang, Xueli Chang, Xueye Mao and all staff in the laboratory for their help.  Appendix A


extfieldsofstudy: Biology Medicine

text:
Construction of Ultradense Linkage Maps with Lep-MAP2: Stickleback F2 Recombinant Crosses as an Example

High-density linkage maps are important tools for genome biology and evolutionary genetics by quantifying the extent of recombination, linkage disequilibrium, and chromosomal rearrangements across chromosomes, sexes, and populations. They provide one of the best ways to validate and refine de novo genome assemblies, with the power to identify errors in assemblies increasing with marker density. However, assembly of high-density linkage maps is still challenging due to software limitations. We describe Lep-MAP2, a software for ultradense genome-wide linkage map construction. Lep-MAP2 can handle various family structures and can account for achiasmatic meiosis to gain linkage map accuracy. Simulations show that Lep-MAP2 outperforms other available mapping software both in computational efficiency and accuracy. When applied to two large F2-generation recombinant crosses between two nine-spined stickleback (Pungitius pungitius) populations, it produced two high-density (∼6 markers/cM) linkage maps containing 18,691 and 20,054 single nucleotide polymorphisms. The two maps showed a high degree of synteny, but female maps were 1.5–2 times longer than male maps in all linkage groups, suggesting genome-wide recombination suppression in males. Comparison with the genome sequence of the three-spined stickleback (Gasterosteus aculeatus) revealed a high degree of interspecific synteny with a low frequency (<5%) of interchromosomal rearrangements. However, a fairly large (ca. 10 Mb) translocation from autosome to sex chromosome was detected in both maps. These results illustrate the utility and novel features of Lep-MAP2 in assembling high-density linkage maps, and their usefulness in revealing evolutionarily interesting properties of genomes, such as strong genome-wide sex bias in recombination rates.

Introduction
Recombination and linkage disequilibrium are two inextricably bound facets of the forces driving haplotype formation, on which natural selection can work (Hill and Robertson 1966;Feldman et al. 1996;Gessler and Xu 2000;Otto and Lenormand 2002;Posada et al. 2002). Thus, assessing the extent of recombination and linkage in the genome of an organism is important for understanding the structural, functional, and evolutionary characteristics of the genome (Wang et al. 2009;Hohenlohe et al. 2011;Kai et al. 2011). A genetic linkage map provides not only the relative order of the markers, but also a direct measure of the extent of recombination and linkage disequilibrium across chromosomes. In sexually reproducing organisms, it also allows identification of the role of each sex in creating novel haplotypes (Broman et al. 1998;Sakamoto et al. 2000;Lenormand 2003;Hedrick 2007). From a structural genomic standpoint, a linkage map provides the data necessary to analyze the presence, location, and relative size of chromosomal rearrangements, such as inversions (Tanksley et al. 1992;Agresti et al. 2000;Bansal et al. 2007). As such, linkage maps can also facilitate de novo genome assembly and validation by enabling the identification of chimeric scaffold constructs (Rastas et al. 2013;Fierst 2015). Furthermore, because haplotype length and persistence are strongly influenced by natural selection, comparisons of the levels of linkage surrounding a locus across multiple populations can provide clues to understand the evolutionary history of the locus in examination (Birky and Walsh 1988;Kreitman and Hudson 1991).
Apart from providing insights into the genomes of target species and populations, high-density linkage maps of nonmodel species can set a strong foundation for comparative genomics and the analysis of synteny across species, providing vital clues for our understanding of genome evolution and speciation (Kulathinal et al. 2009;Larkin et al. 2009;Michalak de Jimenez et al. 2013;Zhang et al. 2014). Additionally, these maps would allow linking phenotypes to specific regions of the genotype (through quantitative trait locus or association mapping analyses), and thus hold the key to understand the genetics of complex phenotypic traits (Paterson et al. 1988;Flint and Mackay 2009;Goddard and Hayes 2009).
Modern development in high-throughput sequencing technologies, such as restriction site associated DNA tags (RADtags; Miller et al. 2007), allow a cost-effective detection of several thousands of single nucleotide polymorphism (SNP) markers in the genome of nonmodel organisms. Compared with microsatellite markers, SNPs have the potential to substantially simplify the creation of linkage maps because they can be potentially genotyped with a greater accuracy and genome coverage than microsatellites (Kruglyak 1997;Slate et al. 2009). However, by substantially increasing the sample space containing the true marker order, large marker data sets increase the computational burden involved with linkage map construction. A number of different approaches have been devised to tackle this problem (van Os, Plet, et al. 2005;van Os, Stam, et al. 2005;Margarido et al. 2007;Tong et al. 2010;Van Ooijen 2011), but effective solutions are still few (Rastas et al. 2013;Liu et al. 2014). For instance, although the linkage mapping program Lep-MAP (Rastas et al. 2013) is in principle capable of creating linkage maps containing tens of thousands of markers, the computational time required for mapping a very large number of markers (!2,000 per chromosome) becomes unfeasible, as it increases cubically with the number of markers per chromosome. Likewise, Lep-MAP does not allow modeling sex-specific recombination rates. Hence, there is a need for improved linkage-mapping software to make efficient use of the high-throughput data provided by new sequencing technologies.
The main aims of this study were 2-fold. First, to introduce and benchmark a substantially improved version of the Lep-MAP (Rastas et al. 2013) software (henceforth Lep-MAP2) capable of creating ultra-high-density linkage maps. Second, to use Lep-MAP2 to construct two high-density linkage maps for nine-spined sticklebacks (Pungitius pungitius) based on large SNP panels obtained using a RAD sequencing approach. The nine-spined stickleback is a nonmodel teleost, closely related (Kawahara et al. 2009;Guo et al. 2013) to the three-spined stickleback (Gasterosteus aculeatus) whose genome has been sequenced (Kingsley and Peichel 2007). Both these species are important models for an increasing amount of evolutionary biology and genetics research (Bell and Foster 1994;McKinnon and Rundle 2002;Kingsley and Peichel 2007;Wootton 2009;Merilä 2013), including the study of sex chromosome evolution (Peichel et al. 2004;Kitano et al. 2009;Ross et al. 2009;Shikano, Herczeg, et al. 2011;Shikano, Natri, et al. 2011;Natri et al. 2013). Hence, we were interested in comparing the degree of synteny and collinearity between nine-spined stickleback linkage maps and the three-spined stickleback genome in order to infer the frequency of inverted and transposed genomic regions, which are suspected to play an important role in both speciation (Flaxman et al. 2014) and local adaptation (Yeaman 2013). Specifically, we were interested in exploring possible heterogeneity in sex-specific recombination rates across the different linkage groups, as well as identifying possible structural rearrangements and recombination heterogeneity in the sex chromosomes.

A Brief Description of Lep-MAP2
Lep-MAP2 software for constructing ultradense linkage maps is based on Lep-MAP (Rastas et al. 2013) with the following novel features and improvements: 1) It takes into account achiasmatic meiosis (recombination in one sex only) and models sex-specific recombination rates, 2) the marker ordering algorithm scales to a much larger number of markers than that in Lep-MAP, 3) it can utilize and gain speed using multicore processors, and 4) the data analyzing pipeline has been improved to ease the map construction. Furthermore, it is largely automated and requires minimal user interaction. It can analyze multiple outbred families simultaneously as well as typical inbred crosses, and can handle all types of genetic marker data (e.g., SNPs, microsatellites).
The input of Lep-MAP2 consists of genotypes of one or several full-sib families (parents and their offspring), given in pre-makeped LINKAGE (Lathrop et al. 1984) pedigree format. The format gives the pedigree information on columns 1-6 and genotypes starting on column 7 onward. Only full-sib type pedigree structure is supported, but data from several types of crosses (e.g., backcrosses) can be treated as full-sibs (supplementary file S1, Supplementary Material online).
The data workflow of Lep-MAP2 with descriptions of five modules included into the program are given in supplementary file S1, Supplementary Material online. Lep-MAP2 software is publicly available together with its source and documentation at http://sourceforge.net/projects/ lepmap2/ (last accessed December 22, 2015).

Simulated Data
Simulations were used to compare the performance of Lep-MAP2 with that of TMAP (Cartwright et al. 2007), JoinMap (Van Ooijen 2011), and HighMap (Liu et al. 2014). To this end, the accuracy (i.e., ability of the software to recover the correct marker order and map length) as well as computational time taken by different programs for one, two, and five family data sets with three different genotyping error rates (0%, 1%, and 5%) were assed. However, we did not evaluate the influence of genotyping errors on the five-family data sets as the run time of TMAP and JoinMap became too prohibitive. Simulated data were created as explained in supplementary file S2, Supplementary Material online. In short, it consisted of 100 individuals from a full-sib family with 10 chromosomes and 300 markers per chromosome. The recombination probability between adjacent markers was set to 0.333% and 0.167% for the father and mother, respectively. The parents were informative (heterozygous) with a probability of 0.5.
Lep-MAP2 was run ten times on each simulated data set, and both the results of the first run (LM1) and the run with the highest likelihood (LM10) are reported. For other software, only one run was conducted. For each run, we computed the Kendall tau (Kendall 1938) between the found and the correct order on the subset of informative markers with detectable recombinations. We also measured the time of each run using the Linux command "time." The timings of LepMAP2 and TMAP were measured with a desktop computer running Linux and having 24 GB of memory and four Intel Core i7-4790 central processing units (CPUs) running at 3.60 GHz frequency. JoinMap was run on a Windows 7 Enterprise computer with 128 GB memory and dual Xeon E5-2640 v3 CPUs running at 2.60 GHz frequency. HighMap was run by the developers of the program itself, because at the time of this study, HighMap was not available for general use. One specific limitation of JoinMap was that it is a 32 bit binary with an obligatory graphical front end. Thus it is difficult to efficiently run multiple jobs or to time them, and the need for direct user input proved to be quite high. Hence, run time comparisons between JoinMap and other software were limited to single-family comparisons and multiple-family comparisons with zero error rate. HighMap runtimes could not be clocked (see above).

Stickleback F 2 Recombinant Crosses
Adult nine-spined sticklebacks were collected from a marine population in Southern Finland (Helsinki, 60 13 0 N, 25 11 0 E) and from two pond populations in northeastern Finland (Rytilampi 66 23 0 N, 29 19 0 E and Pyö reä lampi 66 15 0 N, 29 26 0 E) in 2006 and 2011, respectively. Two F 1 hybrid generations were created by mating a marine female to a pond male from each pond population, and the F 2 generations were generated from the repeated mating of a single full-sib F 1 pair for each hybrid cross. In the case of the first cross (Helsinki-Rytilampi, henceforth HR cross), the resulting 283 F 2 offspring are the same as used in Shikano et al. (2013) and Laine et al. (2013). In the case of the second cross (Helsinki-Pyö reä lampi, henceforth HP cross), 284 F 2 offspring were obtained for the purpose of this study. More details about crossing and rearing procedures used to create HR cross can be found from Shikano et al. (2013). The procedures for setting up and rearing the HP cross were mostly identical to those used for the HR cross. Sex of all the F 2 offspring in both crosses was identified by genotyping all individuals for a sex-linked microsatellite marker (Stn19) as detailed in Shikano, Herczeg, et al. (2011).
This study did not involve human subjects, and our experimental protocols were approved by the National Animal Experiment Board, Finland (permission numbers: ESLH-STSTH223A and STH037A).

DNA Extraction and RAD Library Construction
Genomic DNA was extracted from ethanol-preserved fin clips using the phenol-chloroform method (Taggart et al. 1992). RAD library construction and sequencing were performed by BGI HONGKONG CO., Ltd. Briefly, DNA was fragmented by the restriction enzyme PstI, and DNA fragments of 300-500 bp were gel purified. Illumina sequencing adaptors and libraryspecific barcodes were ligated to the digested DNA fragments, and barcoded RAD samples were then pooled and sequenced on the Illumina HiSeq2000 platform with 45-bp single-end strategy. Twenty-four lanes were used for the HR cross and 30 for the HP cross sequencing. For each cross, grandparents and parents were sequenced in one lane (i.e., four individuals per lane) to increase their sequence coverage, and thereby also the number of mappable SNPs. Adapters and barcodes were eliminated from reads and quality was checked using FastQC (Andrews S. FastQC, http://www.bioinformatics.bbsrc.ac.uk/ projects/fastqc/, last accessed December 22, 2015).

RAD-Tag Sequencing and Genotyping
An in-house pipeline was used to obtain the genotype calls from the raw single-end reads for each individual as follows. As the read length varied between 41 and 45 bp, reads were cropped by keeping only the first 41 bp. The reads of parental individuals without any missing nucleotides in both populations were pooled and identical reads were grouped together. These grouped reads were processed in the descending order of the number of occurrence. Each processed read was added to the sequence list and all its neighbors within edit (Levenshtein) distance of two were removed from the order. Only sequences occurring between 10 and 1,500 times were kept and taken as reference sequences: The sequences occurring less than 10 times were considered sequencing errors and sequences with >1,500 occurrences were likely repeat regions, thus these sequences were discarded. In the end, a total of 712,005 reference sequences, each 41 bp in length, were recorded. We verified that the number of reference sequences agreed with the number of restriction sites in the genome.
All raw reads (cropped to 41 bp) were mapped against the reference sequences with BWA (Li and Durbin 2009), together with SAMtools ), producing a single bam file for each individual.
Individual genotype posterior probabilities, taking into account the read and mapping qualities, were obtained by inhouse scripts (Kvist et al. 2015, their Appendix S1) from the output of SAMtools ) mpileup on the bam files giving multiple alignment of reads for each reference position. Only positions with 2 or more alleles and no more than 20 indels among all individuals were considered. Furthermore, minimum read coverage of 3 was required for more than 158 individuals, which assured that for most (>60%) individuals it would be possible to call genotypes (or at least give some informative posterior) on the remaining markers. The parental genotypes were called by maximizing the likelihood of offspring and parent genotypes using module ParentCall, and then each offspring was called with respect to its called parental genotypes. Parental genotypes were called only if their likelihood was 100 times higher than the second best parental genotype combination. The offspring genotypes were called similarly as the parents.
Genotype calling thus identified 41,730 potential SNP markers shared between the 2 crosses. These potential SNPs were then independently quality checked for linkage map construction and low-quality SNPs were discarded.
Additional genotype data on 226 microsatellite markers were added to the HR data set. These are a subset of markers used in Shikano et al. (2013), but for the purposes of this study 28 markers (listed in supplementary table S1, Supplementary Material online) were omitted from linkage map construction because of high segregation skew or high error estimates (see section on linkage map construction and supplementary table S1, Supplementary Material online).
As a validation of the quality of the RAD-generated SNPs, and also to confirm that sample identities (including sex identity) had not been mixed between sampling and sequencing, for each cross we re-extracted DNA from the original tissue samples and genotyped all individuals for a subset of SNPs (64 for HP and 33 for HR) using the Sequenom (San Diego, USA) platform at the Finnish Institute for Molecular Medicine (FIMM). After discarding SNPs that failed to amplify, we were able to directly compare 58 SNPs for HP and 29 for HR cross between RAD and FIMM calls. These comparisons revealed that no mixing of samples had occurred between sampling and sequencing.

Linkage Map Construction with Lep-MAP2
Linkage group (LG) assignment was obtained using Lep-MAP2 software. First, the SeparateChromosomes module was executed with an logarithm of odds (LOD) score limit of 20 and minimum LG size of 10. Second, singular markers were added to the found LGs using the JoinSingles module with an LOD score limit of 10. Markers with more than 40 (about 14%) missing genotypes were removed from the LGs. Over 6,000 markers in LG12, which corresponds to the sex chromosomes , were informative only for the paternal side in both crosses. To reduce the computational burden, only a common subset of 1,548 of these (paternally informative) markers with at most 10 missing values were kept in this chromosome.
Lep-MAP2 filters out markers by comparing the offspring genotype distribution and the expected Mendelian proportions (segregation distortion test). The default value of dataTolerance = 0.01 was used to filter out highly segregated markers ( 2 test, P < 0.01), thus 1 out of 100 markers should be removed by chance alone. This filtering removed 2,238 and 1,601 SNPs from HP and HP data sets, respectively, and also 28 microsatellite markers from the HR data set as described above. Marker order was determined by allowing different recombination probabilities in both sexes. Ten independent runs were conducted and the marker order with the best likelihood was kept. Only one of the exactly identical markers was used in marker ordering. Furthermore, if there were two markers with identical genotypes but one had more missing genotypes, only the one with less missing genotypes was kept. Lep-MAP2 marks unused markers as duplicated markers and takes their position from the corresponding nonduplicate marker. The final map also included an estimate for genotype error for each marker. The error parameters correspond to the hidden Markov model (HMM) used to model recombinant haplotypes in Lep-MAP2. The recombination rates correspond to the transition parameters in the HMM, whereas the emission probabilities define the error parameters (supplementary file S1, Supplementary Material online). Finally, markers with genotype error rate estimate >0.1 were removed. Few (<0.1%) markers from the ends of LGs were also removed with the criteria that 1) they contributed over 10 cM (per marker) to map length and 2) the parental coverage from their corresponding sequence was above 500 (likely repeat) or below 20 (likely haplotype or sequencing error) on all markers to be removed. The number of removed markers in this last step was 179 SNPs for HP and 174 SNPs plus 8 microsatellites for HR. The HR map was re-evaluated without microsatellites by running Lep-MAP2 on the final order as the initial marker order.
The number of initial LGs found was 23 for HP and 21 for HR. By comparing the linkage maps between HR and HP crosses, it became clear that the two smallest LGs in the HP cross were parts of LGs 2 and 9 in the HR cross. These groups were added to the corresponding LGs and the maps were reevaluated. The reason why these parts were initially separated is in the large gaps shown in figure 2. However, it was also clear that part of HR LG4 was missing from the maps. This part was found by inspecting (SeparateChromosomes with dataTolerance = 0.0001) the markers filtered out based on the segregation distortion test. This part was also separated by a large gap from other markers in LG4, visible in figure 2, and added to the map as described above. Finally, all the maps (and also HR without microsatellites) were polished by running ten independent runs of Lep-MAP2 using the found marker order as the initial marker order.

Comparisons with Three-Spined Stickleback Genome
To compare genomic synteny between nine-and three-spined sticklebacks, the reference sequences with SNP makers in the linkage map were mapped onto the three-spined stickleback genome (Ensembl release-75) by BLAST (Altschul et al. 1997) with an e-value cut-off at 1 Â 10 À5 by considering sequence divergence between nine-and three-spined sticklebacks (Guo et al. 2013). In order to infer whether interchromosome rearrangements have occurred in the nine-or three-spined stickleback lineage, we conducted a BLAST search against the genome sequences of medaka (Oryzias latipes), which is the closest ancestor of sticklebacks with a sequenced genome. The genomic synteny was visualized using CIRCOS (Krzywinski et al. 2009).

Statistical Analyses
LG lengths were log 10 transformed to achieve normality and homoscedasticity between groups. To partition variance in LG lengths to effects of cross and sex, we used the R package lme4 to perform a two-way random-effects ANOVA with interaction. Cross, sex, and their interaction were also used as fixed terms in an ANOVA to determine the significance of their effect on the log 10 of the LG lengths. Correlations between LG lengths were estimated using nonparametric Kendall's rank correlation (tau) to account for the nonnormal distribution of LG lengths. Fisher's exact test was used to test for associations between marker genotypes and sex. Because the two only terms in these analyses were sex and marker genotype, Fisher's exact test is equivalent to a logistic regression. The large number of markers tested makes the use of a P value threshold not useful, because a large number of tests would exceed it purely by chance. Thus, the significance of association between marker genotypes and sex was assessed by comparing the distribution of the observed P value with the P value expected under a null hypothesis of no association.

Lep-MAP2 Performance
Based on the simulations, it is clear that Lep-MAP2 can produce very accurate linkage maps. In single-family simulations, Lep-MAP2 outperformed all other software, both in recovering the correct marker order and the actual map length (150 cM; table 1). Although JoinMap was equally good as Lep-MAP2 in finding the correct marker order, it appeared to be sensitive to map-length inflation due to genotyping errors (table 1). In terms of computational time, LEP-MAP2 was substantially faster than TMAP at all genotyping error rates in single-family comparisons, but slower than JoinMap ( fig. 1). For finding the correct marker order in multiple-family mapping simulations, the differences among Lep-MAP2, TMAP, and JoinMap were negligible: The different programs produced maps of roughly equally high quality (table 1). However, Lep-MAP2 was much faster than TMAP or JoinMap ( fig. 1). Moreover, the speedup obtained by utilizing multicore processors to run Lep-MAP2 was very closely linear to the number of cores used (1 core: 5 h 13 min; 2 cores: 2 h 37 min; 4 cores: 1 h 21 min). Using 4 cores and the same desktop computer as used in the simulations, linkage map construction for 5,000 markers and 16 families was completed in 5 days and 12 h (parameters filterWindow = 10 and polishWindow = 100 were used for extra speedup). The accuracy (jtauj) of obtained solution was 0.999. Based on this, we estimate that the maximum data set size for Lep-MAP2 to analyze in 1 week is about 10,000 markers and 10,000 individuals on a fast computer (with 32 or more cores). The computational time in Lep-MAP2 scales linearly with the number of individuals, and quadratically with the number of markers. Based on the single-family runs, the performance of HighMap was similar to that of TMAP in terms of map order, but better in terms of map length (table 1). However, its performance in comparison with Lep-MAP2 and JoinMap was poor on all fronts (table 1).

Sex-Averaged Linkage Maps for Nine-Spined Sticklebacks
After quality control (QC) thinning we identified an excess of 15,000 markers in both crosses. In the HP cross, we identified 20,054 markers, of which 13,060 (65%) were uniquely informative in building the linkage map (table 2). In the HR cross, 18,691 markers were identified, of which 14,998 (80%) were uniquely informative in building the linkage map (table 2). Overall, 22,761 markers were mapped, 15,984 of which were common to both crosses. The difference in marker numbers is due to the fact that even though all markers were initially chosen from the RAD calls as common between crosses, some markers successfully mapped in only one cross due to a high rate of missing genotypes or a questionable segregation pattern in the other cross.  LGs in both maps (table 2), a number that matches the expected number of chromosomes in the nine-spined stickleback (2n = 42; Ocalewicz et al. 2008). After applying the LG size corrections suggested by Tripathi et al. (2009), the sexaveraged HP and HR cross maps spanned a total of 1980.74 and 2528.96 cM, respectively (table 2). The average genomewide marker density for the sex-averaged HP and HR cross maps were 6.73 and 5.99 markers/cM, respectively (table 2  and

Sex-Specific Linkage Maps
Males and females exhibited a substantial difference in map lengths, with female map length exceeding that of males by a factor of 1.5 (HP) to 1.9 (HR; table 2). The average male map lengths for HP and HR crosses were 1,584.26 and 1,724.27 cM, whereas those for females were 2,335.72 and 3,344.09 cM, respectively (table 2). Accordingly, the marker densities in both male maps (HP 8.43 markers/cM; HR 8.80 markers/cM) were 1.5-2 times higher than those in female maps (HP: 5.83 markers/cM; HR: 4.48 markers/cM; table 2). There were no significant differences in the number and distribution of sex informative markers in the two sexes in either   cross (HP paired t-test, t 20 = À1.25, p = 0.17; HR paired t-test, t 20 = À0.06, P = 0.55), indicating that the map length differences between sexes are not explainable by differences in the number of male only versus female only informative markers. In terms of the relative importance of sex and cross identity, sex explained 58% of variance in LG length (ANOVA: F 1,80 = 97.6, P < 10 À14) whereas the cross effect was much smaller (7% of variance explained; F 1,80 = 18.3, P < 10 À4 ). The sexby-cross interaction was also significant (F 1,80 = 7.4, P = 0.008), but it explained only 8% variance in the data. Hence, in spite of the clear similarities in length of the different LGs across the two crosses, there were also some differences ( fig. 2, table 2, and supplementary table S3, Supplementary Material online). In males, the LG length was uncorrelated between the two crosses (t 19 = 0.24, P = 0.14), whereas in females the similarity was higher (t 19 = 0.44, P = 0.005), indicating that in females the LG length order in one cross provides a reasonable approximation of the LG length order in the other cross. For the sex-averaged maps the concordance in the LG length distribution was good, albeit far from perfect (t 19 = 0.57, P < 0.001).

Association between Markers and Sex
Sex determination in Gasterosteidae varies across species, with evidence of recent evolution of sex chromosomes (Ross et al. 2009). Association analyses allowed assessment of evidence for sex-associated loci outside the sex chromosomes (LG12). After determining each marker/sex association P value, LogQQ plots were used to visualize the results at the genome-wide level ( fig. 3a-d). Although a plot of the P values including all LGs shows an impressive deviation from the null hypothesis of no association in both crosses ( fig. 3a and b), removal of LG12 totally erases any signature of sex-marker associations ( fig. 3c and d). These results suggest that although the markers on LG12 are strongly associated with sex, none of the markers on any of the other LGs are sex linked. Of particular note is the pattern of association between markers and sex on LG12 in both crosses ( fig. 3e and f). When plotting the Àlog 10 of the P value of association with given SNP and sex against map position, two interesting observations can be made. The first is that the pattern of distribution of the P values segregates at three different levels, shown in figure 3e and f as "lines." These lines correspond to the three possible SNP genotype arrangements in the F 1 parents: The "0 line" corresponds to the loci where the F 1 male had the same allele on the X and the Y chromosome, and thus male and female F 2 are equally divided across the genotypes; the "middle line" corresponds to the loci where both F 1 parents are heterozygotes, and thus 50% of the male and 50% of the female F 2 are identified by their genotype; and the "top line" corresponds to the loci where the female is a homozygote and the male is a heterozygote with a distinctive allele on the Y, giving a unique genotype to the males and to the females in F 2 . Thus the F 1 genotype completely explains the P value pattern we observe on LG12. The second observation is that, despite a strong association between sex and genotype, the part of LG12 that is syntenic with the three-spined stickleback Chr7 (see below) does not show an association between markers and sex in its distal part. This finding suggests that while the translocated part is still recombining, recombination in the ancestral part has almost (but not entirely) ceased ( fig. 3e and  f). This inference is supported by comparing linkage map lengths and marker numbers (a proxy of physical size of the chromosomes) between the ancestral and translocated parts of LG12: In spite of being shorter in terms of linkage map length in males ( fig. 4), the ancestral part of the chromosome has a much higher number of markers (HP male 2,077 markers; HR male 785) than the novel translocation syntenic with the three-spined stickleback Chr7 (HP male 483 markers, HR male 424 markers). Accordingly, the marker density on the ancestral part of LG12 is much higher than that in the translocated part (HP males 93.8 vs. 7.7 markers/cM, HR males 30.2 vs. 5.4 markers/cM). Note that recombination on the ancestral portion of LG12 is not completely suppressed, because this region does not collapse into one single fixed haplotype. Interestingly, recombination in this region is lower than in the novel part in females (marker density: HP 21.1 vs. 8.2 marker/cM; HR 9.7 vs. 5.4 marker/cM), despite the fact that the ancestral part of the LG is actually bigger (in terms of recombination) than the novel part ( fig. 4).

Genomic Synteny between Nine-and Three-Spined Sticklebacks
BLAST searches of the 34,015 nine-spined stickleback reference sequences with 41,730 SNPs against the three-spined stickleback genome sequence indicated a high degree of genomic synteny. A total of 11,030 of the 34,015 nine-spined stickleback reference sequences yielded high-scoring BLAST hits and 10,320 of these had unique high-scoring BLAST hits on the three-spined stickleback genome, with 6,506 located in coding regions. In all further comparisons, only these 10,320 sequences were utilized.
In the sex-average linkage map of the HP cross, 5,229 reference sequences with 5,732 SNPs mapped to the threespined stickleback genome (table 3) indicating a high degree of genomic synteny between nine-and three-spined sticklebacks ( fig. 5a). Although most of the 5,229 reference sequences were located in the syntenic LG pairs, 244 (4.3%) showed interchromosomal rearrangements between the species (table 3). For example, 120 reference sequences located on three-spined stickleback Chr7 spanning around 10 Mb (range 125,492-10,364,802 bp) were mapped to LG12 in the nine-spined stickleback, together with 601 reference sequences that mapped on Chr12 in the three-spined stickleback genome (fig. 5a).
The sex-average nine-spined stickleback linkage map of the HR cross was very similar to that of the HP cross in respect to the syntenic relationships with the three-spined stickleback genome ( fig. 5b). For instance, 4,791 reference sequences with 5,241 SNPs mapped to the three-spined stickleback genome mostly in perfect synteny, and only 206 (4.3%) reference sequences showed interchromosomal rearrangements between nine-and three-spined sticklebacks. As in the case of the HR cross, half (103) of these reference sequences were located on three-spined stickleback Chr7 covering a 10 Mb (range 86,906-10,364,802 bp) region of this chromosome.
Of the reference sequences showing interchromosomal rearrangements, 174 were found in both HP and HR maps, and 96 (55%) were located on Chr7 in the three-spined stickleback genome while mapping on LG12 in both nine-spined stickleback linkage maps. In addition, we found that both HP and HR maps harbored sequences which remain unassembled in the three-spined stickleback genome. The HP map included 103 such scaffolds, and the HR map included 96 similar scaffolds (supplementary table S4, Supplementary Material online). Considering the high genomic synteny between nine-and three-spined sticklebacks, these findings might help to improve the three-spined stickleback genome assembly.

Discussion
High-density linkage maps, such as those constructed here using Lep-MAP2, provide means to gain insights to genome-wide linkage and recombination patterns and thereby the structural, functional, and evolutionary characteristics of the genome itself (Wang et al. 2009;Kai et al. 2011;Hohenlohe et al. 2011). A quick survey of the current literature (supplementary table S5, Supplementary Material online) shows that although marker numbers in recent linkage mapping studies continue to increase (median = 3,677 markers), they are still modest compared with the maps we produced with the Lep-MAP2 ( fig. 6a). Also in terms of marker density, our maps are among the highest published to date ( fig. 6b). Apart from suggesting genome-wide recombination suppression in male nine-spined sticklebacks, these maps support the suggestion ) that a translocation of an autosomal chromosome arm to a sex chromosome has taken place after the nine-spined stickleback diverged from the three-spined stickleback. The results further suggest that the translocated part of this neo-sex chromosome is still recombining, whereas the ancestral part has nearly, but not entirely, ceased to do so. Furthermore, the comparative genomic analyses revealed a high degree of synteny between three-and nine-spined stickleback genomes, with some evidence of infrequent interchromosomal rearrangements. In the following, we will discuss these findings and their implications to our understanding of stickleback genome evolution, highlighting the value of these new ultra-high-density linkage maps as genomic resources of broad utility. We also discuss the utility and advantages of Lep-MAP2 in the construction of high-density linkage maps. The two high-density second-generation SNP-based linkage maps, constructed using the RAD-seq approach, provide substantial improvements over the previously available microsatellite-based maps for nine-spined sticklebacks Shikano et al. 2013). The basic structure of the SNP-based HR map was in agreement with the microsatellite-only HR linkage map: In both maps the microsatellites mapped to the same LGs, and their overall order was comparable between the two maps (Wilcoxon rank-sum test, 198 matched observations, W = 19,678, P = 0.34). However, the new maps increased the overall coverage of the LGs, which is not surprising in the view that the new maps contained 83-105 times more markers than the previous maps , illustrating the effectiveness of the RAD-seq approach in SNP discovery, genotyping, and linkage mapping in a nonmodel organism (Etter and Johnson 2012). Hence, the high marker densities and more even distribution of markers across the different LGs in the SNP-based maps yielded a far more refined image of the genetic landscape of the ninespined stickleback genome than that provided by first-generation microsatellite-based linkage maps. Both these features-increased marker density and coverage-helped us to not only verify the high degree of synteny in genomes of nine-and three-spined sticklebacks, but also to detect genomic rearrangements that have occurred during the evolutionary history of sticklebacks. Namely, earlier comparative genomic analyses between the three-and nine-spined sticklebacks have provided some preliminary insights into the genome evolution between these two model species (Shikano et al. 2010;Guo et al. 2013), and our results refine this picture. Based on their estimated divergence of around 13 Ma (Bell et al. 2009), a high degree of synteny was expected. This expectation was fulfilled: >5,000 SNPs were uniquely mapped onto the three-spined stickleback genome in each of the maps with a high degree of synteny. About 65% of these SNPs were located within coding regions, possibly because substitution rates are higher in noncoding than in coding regions of genomes (Guo et al. 2013). We also discovered that many scaffold sequences which remained unassembled to chromosomes in the three-spined stickleback genome assembly (Jones et al. 2012) showed strong linkage to specific LGs in both of our nine-spined stickleback linkage maps, suggesting that the utility of the high-density linkage maps in genome assembly is not limited to target species, but could also aid in assembling genomes of closely related species.
An earlier study identified a possible rearrangement between an autosome (LG7) and the sex chromosome (LG12) in the nine-spined stickleback . Results of this study confirmed this finding with two independent crosses, and provided higher resolution information about the size and location of this rearrangement. Namely, it appears that a chromosomal segment corresponding to 36% (ca. 10 Mb) of Chr7 in three-spined sticklebacks has fused to LG12 in the nine-spined stickleback. Frequent chromosomal rearrangements in fish genomes are well known (Mank et al. 2006;Mank and Avise 2009), and especially sex chromosome fusions with autosomes appear to be common (Mank and Avise 2009; Kitano and Peichel 2012). Based on comparative analyses of the three-spined stickleback and medaka genomes, no interchromosomal rearrangement were found in these LGs. Therefore, it appears that the chromosome rearrangement between LG7 and LG12 occurred in the ninespined stickleback after it diverged from the three-spined stickleback. In the linkage map from North American ninespined sticklebacks, the segmental part of LG7 that linked to LG12 in our study did not show linkage to either LG12 or to the remaining part of LG7 ). This discrepancy could be explained if this rearrangement is not present in the North American populations of nine-spined sticklebacks. Alternatively, the small number (=120) of individuals utilized in the North American study ) might have rendered the power to detect the rearrangement low. Further cytogenetic analyses and/or more refined genetic maps based on larger sample sizes from North American and other ninespined stickleback clades can clarify this issue. Irrespective of the situation in other populations, the occurrence of this rearrangement among the eastern European nine-spined stickleback populations is now an undisputed fact.
Apart from this translocation event, interchromosomal rearrangements were infrequent (<5%) in both maps. The marker order rearrangements in each LG varied substantially, from almost complete homology to cases where several inversion and translocation events between the two species were indicated to have happened. These findings align with the conjecture that although synteny among species is often conserved, the gene order within syntenic blocks is frequently changed (Woods et al. 2005;Kasahara et al. 2007). In addition, our findings support the observations that intrachromosomal rearrangements are fixed more frequently than translocations among nonhomologous chromosomes (Woods et al. 2005;Sé mon and Wolfe 2007). Furthermore, studies utilizing physical genetic maps of closely related species or populations show that most detected inversions tend to be of small size (<1 kb; Kirkpatrick 2010; see also Feuk et al. 2005;Jones et al. 2012). In our study, the detected inversions were all relatively large, but this may be partly explained by our approach, which did not allow the detection of small inversions. Single SNP rearrangements are hard to interpret in our data as the linkage map is based on recombination frequencies, rather than on physical genetic positions. Therefore, a more refined estimate of inversions that are fixed among these two stickleback species cannot be performed until the whole genome sequence is available for the nine-spined stickleback. Also, we wish to emphasize that although our linkage maps were de novo assembled, all the comparative genomic analyses rely heavily on the three-spined stickleback reference genome sequence, which has been shown to contain assembly errors (Kasahara et al. 2007;Glazer et al. 2015). Although we corrected for known errors in our analyses, it is possible that some of the rearrangements we have discovered (or overlooked) might still owe to problems in the reference sequence. Although we have no a priori reason to believe that this would have biased our results-especially given that three-spined stickleback studies have found evidence for frequent inversions even at the intra-specific level (Deagle et al. 2012)only access to the physical map of the nine-spined stickleback can help to eliminate doubts about the fixed inversions among these two species.
We discovered evidence for dramatic sexual dimorphism (SD) in recombination frequency at the genome-wide level: In both crosses, female maps were substantially (1.5-1.9 times) longer than male maps, and these differences were not confined to the sex LG but also occurred in autosomal LGs. Likewise, linkage blocks were on average shorter in females than in males, both in autosomal and sex chromosomal LGs. These observations suggest genome-wide recombination suppression in males, which are the heterogametic sex in this species. Similar observations are common in a wide variety of taxa (Burt et al. 1991;Lenormand and Dutheil 2005;Hedrick 2007;Brandvain and Coop 2012), including many fish species showing even more extreme SD in recombination frequency (Onchorynchys mykiss: 3.2, Sakamoto et al. 2000;Danio rerio: 2.7;Singer et al. 2002). Although the ultimate reasons for sex-specific recombination rates remains elusive (Lenormand and Dutheil 2005;Hedrick 2007;Brandvain and Coop 2012), the fact that it occurs implicates an important role for female meiosis in generating genetic variability in many species, including the nine-spined stickleback. The higher female recombination frequency also means that for some specific gene-mapping applications, the choice of study design with emphasis on either segregation of male or femalespecific variability may be desirable (Singer et al. 2002).
Sex chromosomes in teleost fish show rapid evolution (Charlesworth 2004;Volff et al. 2007;Kondo et al 2009;Natri et al 2013;Shikano et al. 2013), to the point that closely related stickleback (Gasterosteidae) species have different sex chromosome systems (Ross et al. 2009). Our results confirmed that LG12 is the sex chromosome in the nine-spined stickleback, and that there are no sex-linked markers on any of the other LGs. Hence, sex determination in this species is caused by a well-defined sex chromosome, and not by polygenic factors spread across the genome. Of particular note is that the part of LG12 (the neo-sex chromosome) that is syntenic with the three-spined stickleback Chr7 showed a clearly different pattern of association with sex as compared with the rest of this LG. In particular, the P value pattern in this part of LG12 shows the presence of a transition between the ancestral LG12 and the syntenic portion. The gradual decay of association between sex and genotype across LG12 is clearly indicative of total recombination suppression in the ancestral part of LG12, whereas most of the syntenic part of the LG is still recombining. Hence, the data suggest that the sex chromosomes in the nine-spined stickleback are undergoing a rapid evolution making it an ideal model to understand both the evolution of sex determination systems, as well as patterns and processes occurring in early stages of sex chromosome evolution.
Although marker numbers and densities in our maps were very high in comparison with most earlier linkage maps-even in comparison with those created with RAD-seq approach (e.g., 8,257 SNPs, Gonen et al. 2014; 1,622 SNPs, 1.16 markers/cM, Kakioka et al. 2013; 755 SNPs, 0.5 markers/cM, Recknagel et al. 2013; fig. 6)-we did not observe a perfectly uniform distribution of markers across maps in either of the crosses or sexes. The presence of large intermarker gaps observed in the maps is especially interesting because the marker densities in our maps were higher than what could be fully resolved with the number of recombination events associated with an F 2 -cross, and many markers mapped to the same linkage positions despite the fact that these markers were individually resolved in the sequence assembly. Hence, one possible explanation for these large gaps is that they signal the presence of recombination hotspots across the P. pungitius genome. This is not implausible because recombination is known to be uneven across the genome and to be mostly affected by recombination hotspots (Auton and McVean 2007). When compared with an exponential distribution these intermarker gaps do not seem to arise by chance alone (HP female map: D = 0.072, p = 0.003; HP sex-averaged map: D = 0.0838, p = 0.0002; HR sex-averaged map: D = 0.0644, p = 0.002), but at present, given the lack of a physical map, it is difficult to determine whether recombination or another cause, such as a localized absence of RAD-seq derived loci, is the cause of these large intermarker gaps. Because Lep-MAP2 removes markers that show a significant deviation from the Mendelian proportions expected in an F 2 , and because we were able to select a high number of markers that fulfilled this condition for our maps, our data are not suited to study segregation distortion. We note though that for three LGs (LG2 and LG9 in the HP cross and LG4 in the HR cross) the markers on one end of the LG, while mapping together, mapped almost completely independently compared with the rest of the LG and were placed on maps through multiple mapping iterations. We cannot currently offer an explanation for the difficulties in mapping markers at these particular LGs.
Finally, being a modern implementation of a linkage mapping software applicable to sequencing-based data sets, Lep-MAP2 aims to directly tackle the software limitations which are holding back our ability to create ultradense maps from high-throughput data sets. The increased data output by genotyping-by-sequencing technologies have rendered many linkage mapping software obsolete, and there is an increasing demand for high-density linkage maps as tools for genome assembly validation (Fierst 2015). In addition to its ability to deal with large data sets, our simulations on smaller data show that Lep-MAP2 produces results that are both highly reliable and precise, and much more so than those of other available linkage mapping software. Because Lep-MAP2 is implemented fully in JAVA, it is truly machine independent, and thus compatible with any computing choices or needs of the user. It can utilize multiple cores of typical CPUs and is easily run on a computing cluster without any direct user intervention. These features together with its ability to handle high-throughput data-such as the one analyzed in this article-place Lep-MAP2 among the top state-of-the-art linkage mapping software currently available.
In conclusion, apart from giving novel insights into the genomic architecture of the nine-spined stickleback, the ultra-high-density linkage maps described in this study illustrate the power and utility of Lep-MAP2 software in handling large marker data sets, as well as in modeling sexspecific recombination. The homology analyses revealed a high degree of interchromosomal synteny between threeand nine-spined sticklebacks, but also that inversions have frequently occurred during their divergence. The results also confirmed ) the presence of an interchromosomal rearrangement that has led to the formation of a neo-sex chromosome in the nine-spined stickleback, as well as fairly strong genome-wide recombination suppression in male nine-spined sticklebacks. The constructed maps should also provide useful resources for further QTL mapping and comparative genomic analyses, as well as aid in the assembly of the nine-spined stickleback genome sequence. We note that our maps can also prove to be a valuable resource for improving the three-spined stickleback genome assembly: Given the high degree of synteny between the genomes of these two species, our finding of unmapped scaffolds in the three-spined stickleback assembly suggests their likely location in the three-spined stickleback genome. Hence, we envision that the results, insights, and resources created in this study will not only be useful for future genomic studies of nine-spined sticklebacks, but also for those of other closely related taxa.


extfieldsofstudy:
Computer Science
Biology
Medicine

text:
Systematic noise degrades gene co-expression signals but can be corrected

Background In the past decade, the identification of gene co-expression has become a routine part of the analysis of high-dimensional microarray data. Gene co-expression, which is mostly detected via the Pearson correlation coefficient, has played an important role in the discovery of molecular pathways and networks. Unfortunately, the presence of systematic noise in high-dimensional microarray datasets corrupts estimates of gene co-expression. Removing systematic noise from microarray data is therefore crucial. Many cleaning approaches for microarray data exist, however these methods are aimed towards improving differential expression analysis and their performances have been primarily tested for this application. To our knowledge, the performances of these approaches have never been systematically compared in the context of gene co-expression estimation. Results Using simulations we demonstrate that standard cleaning procedures, such as background correction and quantile normalization, fail to adequately remove systematic noise that affects gene co-expression and at times further degrade true gene co-expression. Instead we show that a global version of removal of unwanted variation (RUV), a data-driven approach, removes systematic noise but also allows the estimation of the true underlying gene-gene correlations. We compare the performance of all noise removal methods when applied to five large published datasets on gene expression in the human brain. RUV retrieves the highest gene co-expression values for sets of genes known to interact, but also provides the greatest consistency across all five datasets. We apply the method to prioritize epileptic encephalopathy candidate genes. Conclusions Our work raises serious concerns about the quality of many published gene co-expression analyses. RUV provides an efficient and flexible way to remove systematic noise from high-dimensional microarray datasets when the objective is gene co-expression analysis. The RUV method as applicable in the context of gene-gene correlation estimation is available as a BioconductoR-package: RUVcorr. Electronic supplementary material The online version of this article (doi:10.1186/s12859-015-0745-3) contains supplementary material, which is available to authorized users.

2 encephalopathy. Thus, a good set of negative controls would be a set of genes that can be assumed to play no role in the disease process.
One possible set of negative control genes would be so-called housekeeping genes. Housekeeping genes are involved in the basic maintenance of the cell, and their expression levels tend to be constant across a wide range of conditions. As noted in the main text, published lists of housekeeping genes are available, and one of the methods used to discover these genes is to identify genes whose expression levels are observed to be fairly constant [2].
In selecting a set of negative controls for the analyses in this paper, we elected to not use the housekeeping genes and simply used a set of genes that we observed to be fairly constant within our own datasets. Selecting our own set of negative controls in this manner also allowed us to ensure that we selected a set of genes with a representative range of expression levels. More specifically, we selected our negative controls as follows: First, we calculated the inter-quantile range of every gene. We then binned the genes according to their mean expression levels. Finally, we selected a fraction of the genes from each bin with the lowest IQR and designated these genes as our negative controls.
In the analyses of this paper, this method for selecting negative control genes seems to have produced good results. However, this method for selecting negative control genes may not be best.
In particular, we have not investigated the performance of alternative sets of negative controls (e.g. published lists of housekeeping genes); such an investigation is beyond the scope of this paper. Housekeeping genes may also not be ideal as they have unknown correlation structure which may not be representative of background co-expression.
We would also like to point out some potential pitfalls associated with using this method for selecting negative controls. The primary concern is that the genes we select may not actually be negative controls. Recall our model, Y = Xβ + W α + . We note that there are three sources of variation that contribute to the overall observed variation in gene expression levels the variation of interest (Xβ), the unwanted variation (W α), and random noise ( ). By selecting genes that show an overall low level of variation, we hope to exclude genes that exhibit any variation of interest. However, it is also possible that we might select genes that simply have an unusually low level of unwanted variation or random noise. A more pernicious problem may arise if X and W are negatively correlated. In such cases, the variation of interest and the unwanted variation may partially cancel out and thus may prevent the effective removal of the systematic noise.
Unfortunately, it is impossible to assess whether these caveats apply in our real datasets.
Finally, note that RUV-random assumes a mixed-effects model, in which the W α term is mod-3 eled as random. The negative controls are used to estimate the covariance of this term. However, by systematically selecting genes observed to have a low level of variation as our negative controls, we may bias our estimate of the covariance of W α.

RLE Plots for Large Data Sets
RLE plots are valuable graphical tools for estimating parameters required for the application of RUV-random. Unfortunately, for large data sets RLE plots loose their appeal, as it is impractical to compare hundreds or even thousands of boxplots. Here, we present a new RLE plot that neatly summarizes the information for all samples into one plot and additionally offers the possibility of adding information on known sources of systematic noise.
The original RLE plot is a boxplot of the devation from the median gene expression calculated for each gene in a sample. Using the width of the boxplot against the median deviation, it is possible to plot all samples in the same plot. While this still allows the easy identification of outliers, it has the added benefit of providing a study-wide perspective. In particular, samples can be colored according to the levels of a known source of systematic noise (compare Figure 1). This is implemented in the function RLEPlot in the package RUVcorr. The left-hand panel shows novel RLE plot for RUV-random treated data and the right-hand panel displays the untreated data. The samples are colored by which batch they belong to. The corresponding legend can be found on the far left. For a more thorough description of the novel RLE plot see the text.

FIVE LARGE DATASETS ON GENE EXPRESSION IN THE HUMAN BRAIN
We downloaded the raw (i.e. non-normalized) data for each study. Data for the Colantuoni et al study [3] (acession number: GSE30272) and the Hernandez et al study [4]  Due to the differences in study design and technology preprocessing protocols varied. However, for all datasets we applied background correction [8] and then quantile-normalization [9] or in case of the Colantuoni et al dataset loess-normalization [10]. In the following we briefly describe the preprocessing steps taken for each study. Both studies relied on Agilent technology to measure gene expression (see [11] for more details). Analysis was conducted using the R package limma.  S1 Table for the names of the removed arrays). These were consequently removed from further analysis. We then proceeded to apply the standard normalization approaches. Since the datasets in this study included some replicate samples. We averaged the samples for these. We then removed transcripts that were nonhuman or not annotated as well as all technical controls. Transcripts with multiple measurements were represented by the measurement with the largest median value across all samples. Finally, we log 2 transformed all gene expression values.

Dataset by Kang et al
In this study, gene expression was measured using the Affymertix Human Exon 1.0 ST array.
Analysis was conducted with the R package aroma.affymetrix [13]. First we checked array density plots and spot images identifying 11 abnormal or failed arrays that we removed from further analysis (see S1 Table for the names of the removed arrays). We then normalized the data using the approaches outlined above. After this probe level modelling was applied to the differently treated datasets including the untreated data. We then removed transcripts that were non-human or not annotated as well as all technical controls introduced by the manufacturer. Transcripts with multiple measurements were represented by the measurement with the largest median value across all samples. Finally, we log 2 transformed all gene expression values.

Dataset by Colantuoni et al
In this study, gene expression was measured using a custom Illumina two-color microarray with a common reference. Using the R package limma [14] we calculated the M values for further analysis, which is standard procedure when dealing with two-color microarrays with a common reference. These were used to create array density plots as well as PCA plots. In combination with MA-plots for each sample, these plots helped us to identify three failed or abnormal arrays which were consequently removed from further analysis (see S1 In this study, gene expression was measured using an Illumina HumanHT-12 V3.0 beadchip. Analysis was conducted with the R package limma. First we checked array density plots, MAplots and PCA-plots identifying 5 abnormal or failed arrays that we removed from further analysis (see S1 Table for the names of the removed arrays). We then normalized the data using the approaches outlined above. We then removed transcripts that were non-human or annotated as 'bad matches' by the manufacturer. Transcripts with multiple measurements were represented by the measurement with the largest median value across all samples. Finally, we log 2 transformed all gene expression values.

PRIORITZING EE CANDIDATE GENES
In 2014, Oliver et al provided new insights into the mechanism of EE using gene prioritization [15]. In this paper, we follow their approach with few modifications and some updates concerning known EE genes and candidates.

Gene Prioritization Method
Candidate genes are prioritized when they are co-expressed with at least one known EE gene.
Thus, for each candidate and each known EE genes combination, co-expression needs to be established. Oliver et al used the PCC (or a weighted version introduced by Horvath [16] in case of multiple samples per genes) in order to establish the amount of co-expression. In order to avoid declaring spuriously correlated genes co-expressed, they used a thresholding approach. We also make use of the PCC combined with a thresholding approach. Here, we set the threshold for an absolute value of the PCC that corresponds to the top 0.2 of ranked random genes. These random genes were ranked according to their maximum absolute correlation with any of the known EE genes. This process was repeated 1000 times and the average maximum absolute correlation was calculated and set as a threshold. Note that the size of the randomly selected genes corresponded to the number of candidates in the various datasets.

Known EE Genes and Candidates
The list of known EE genes used in the Oliver et al paper was updated from 29 known EE genes to 33 confirmed EE genes (as known in September 2014). Candidate genes included all candidates used in the Oliver et al paper, well-established genes for other types of epilepsies and recent significant findings from the ILAE meta-analysis of genome-wide association studies [17].
Genes 500kb up-or downstream of the location of these findings were included in the list. Table   I gives the number of candidates and known EE genes found in each dataset. Note that for the application of RUV-random, all known EE genes and candidates were excluded from being selected as negative control genes.

SIMULATION
The simulation framework is based on the RUV linear model, which is described in the Section Methods and Materials of the accompanying paper. For the purposes of investigating the performance of different cleaning procedures in the context of gene-gene correlations, it is necessary to known the true underlying correlation structure between the simulated genes. Therefore, the fundamental concept of this simulation is that Σ = Cor(Xβ) represents the genuine gene-gene correlation structure induced by the biological signal. This setup allows control of the average absolute value of the gene-gene correlations through the dimensionality of X and β, p. In the following we explain the simulation setup in detail.

Simulation Setup
For the simulation the parameters of the RUV linear model where simulated in the following manner: • X is a m × p matrix with x i,j ∼ N (0, 1) for i = 1, ..., m and j = 1, ..., p.
• β is a p × n matrix with where n c is the number of negative control genes.
Note that we define Σ i,j = 0 when either i = (p − n c + 1), ..., p or j = (p − n c + 1), ..., p. To simulate correlation between X and W define X = X|W is m × (k + p) matrix with rows X i ∼ N (0, Σ X ) for i = 1, .., m. The covariance matrix is defined by Σ X = with g denoting the dimensionality of the shared subspace of X and W . Hence, the parameter g controls the average correlation between the columns of X and W . If g is increased then the average correlation increases.
In our R-package RUVcorr the sizes of all parameters presented above can be varied. For the simulation scenarios presented here we used the following values, unless they were varied for specific simulations: • Number of samples m = 1000 • Number of genes n = 2500 of which n c = 2000 are negative control genes • Level of systematic noiseα = 2 • Variance of random noise˜ = 0.01 • Dimensionality of systematic noise k = 10 • Dimensionality of factor of interest p = 5 • No correlation between X and W Note that in the application of RUV-random generally all 2000 simulated negative control genes were used. We also assumed that the values of k and ν were known.

Further Simulation Scenarios
We investigated the performance of RUV-random when the systematic noise level was varied (compare Table II). Interstingly, the performance of RUV-random in the simulations is absolutely unchanged no matter how much systematic noise is present in the data. When combining RUV-random with other cleaning approaches we had to find the optimal choice for the parameters k and ν first. Interestingly, the optimal choices for k were generally very large indicating that the application of QN or BC introduced new systematic noise that RUVrandom struggles to remove. Generally, RUV-random on its own did much better than combining approaches (compare Table III). However, it needs to be remembered that we did not simulate measured background noise.

Incorrectly Specified Negative Control Genes
We simulated the effect of incorrectly specified negative control genes. For this different proportions of our 2000 negative controls were selected to be genes that were simulated to have co-expression with other genes. Figure 2 demonstrates that even 50% incorrectly specified negative control genes had little effect on the performance of RUV-random. Indeed, it seems as long as there are enough genuine negative control genes RUV-random perfoms well. Thus, we investigated how many negative control genes are necessary for RUV-random to perform well. Figure 3 shows that with a number twice as large as k correlations estimated from the RUV-random corrected simulated data are already very close to the truth. However, we would still strongly recommend using as many negative controls as possible.

REAL DATA
In this Section we briefly comment on the application of RUV-random to the real datasets and we provide some more information on the results of the prioritzation of the candidate EE genes.

Application of RUV-random to the Real Datasets
As discussed earlier, the application of RUV-random requires the researcher to choose various parameters. For this paper we limit ourselves to giving our parameter choices for each dataset (compare Table IV), but do not provide detailed explanation for our choices. Briefly, we tried to selectk and ν in a way to avoid over-correction and removal of the factor of interest, while at the same trying to ensure that all systematic noise from known sources was removed. To this end we used the methods discussed in the Section Methods and Materials of the paper and Section 1 1.1. For all datasets we used 4000 empirically chosen negative control genes for RUV-random. Note that the performance of RUV-random is determined by the absolute number of negative control genes and not by the proportion of negative control genes to all genes in the dataset.