A high-density genetic map construction and sex-related loci identification in Chinese Giant salamander

The Chinese giant salamander Andrias davidianus is an important amphibian species in China because of its increasing economic value, protection status and special evolutionary position from aquatic to terrestrial animal. Its large genome presents challenges to genetic research. Genetic linkage mapping is an important tool for genome assembly and determination of phenotype-related loci. In this study, we constructed a high-density genetic linkage map using ddRAD sequencing technology to obtain SNP genotyping data of members from an full-sib family which sex had been determined. A total of 10,896 markers were grouped and oriented into 30 linkage groups, representing 30 chromosomes of A. davidianus. The genetic length of LGs ranged from 17.61 cM (LG30) to 280.81 cM (LG1), with a mean inter-locus distance ranging from 0.11(LG3) to 0.48 cM (LG26). The total genetic map length was 2643.10 cM with an average inter-locus distance of 0.24 cM. Three sex-related loci and four sex-related markers were found on LG6 and LG23, respectively. We constructed the first High-density genetic linkage map and identified three sex-related loci in the Chinese giant salamander. Current results are expected to be a useful tool for future genomic studies aiming at the marker-assisted breeding of the species.


Background
The Chinese giant salamander Andrias davidianus, historically widely distributed in China, is the world's largest extant amphibian. Due to its habitat destruction and hunting by man, it has been classified as endangered by the International Union for Conservation of Nature and Nature Resources since the 1980s [1]. To conserve the wild resources of giant salamander, artificial breeding technology was developed and succeeded [2,3]. Because of its palatability and nutritional value, the Chinese salamander is commercially cultured and amphibians are now allowed to reach the market [4]. It is considered an important species to study because scientists see it as a living fossil present on Earth since 350 million years ago and linking aquatic to terrestrial organisms, with complex sex determination and differentiation characteristics. Until now, most studies on A. davidianus focused on the immune mechanisms [5][6][7][8], and only few papers were reported for the sex-differentiation issues [9][10][11].
Previous studies reported that the Chinese giant salamander has a huge genome (about 50 Gb) which was a big challenge to assemble [12], and to this direction the construction of a genetic linkage map will be a valuable tool. In the previous study, Simple sequence repeats (SSR) was the mainly basic tools to generate a linkage map in various species. With the development of biotechnology, single nucleotide polymorphisms (SNP) has become the most popular tools to construct the highdensity linkage maps, because of their widely distribution in the genome and development of next genomic sequence technology. Recently, high-density genetic linkage maps have become increasingly important and widespread in cultured aquatic species. In fish, genetic linkage maps have been reported in Salmo salar [13,14], Cyprinus carpio [15], Paralichthys olivaceus [16], Megalobrama amblycephala [17], Larimichthys polyactis [18,19], Oreochromis niloticus [20], Ictalurus punctatus [21], and Lates calcarifer [22], all of which possess relatively small genomes. In amphibian, genetic linkage maps have been constructed in several species including Rana temporaria [23], Hyla arborea [24], and Lissotriton newts [25]. However, high-density genetics maps are not yet available in amphibians.
Quantitative trait locus (QTL) mapping is an effective means to relate phenotype to the genome. In previous study, QTL mapping has been used in various species for different traits location including loci related to ecomic phenotype [26] and resistance to disease [27,28]. The composite interval mapping (CIM) is a common method for mapping QTLs based on linkage map [29] and it was used to located binary trait such as sex [30] and metamorphosis-related traits loci in fish [31]. To identify the sex chromosome of A. davidianus without a genome, the CIM method was effective approach. In the previous study, we explored the female-specific marker to identify the genetic sex of the salamander which was important to detect the sex reversal salamander treated by high temperature and sex hormone in the study of sex determination mechanism [32,33]. we also characterized sex related gene in the sex differentiation [34,35]. Studies on the genetic level are necessary to identify sex related regions and ascertain the sex determination mechanism. In the present study, we used a single cultured A. davidianus full-sib family to construct a highdensity linkage map based on the ddRAD sequencing. Using the constructed linkage map. We identified three sex-related loci in different linkage groups. These results provide an effective genetic tool for A. davidianus genome assembly and marker-assisted breeding.

Library construction and sequencing
Ninety-seven RAD-seq libraries from two parents and 95 offspring were constructed and sequenced on an Illumina HiSeq X platform to generate raw reads. Based on the data trimming, 10.49 billion reads comprising ∼1516.07 Gb of sequencing data were individually partitioned into RAD tags according to their multiplex identifiers (GenBank accession No. SRP155453). Finally, female and male parent datasets which contained 210, 000,000 and 190,000,000 clean reads, respectively, were correspondingly partitioned into 1,887,649 and 6,378,461 RAD tags. A total of 1458.71 Gb of data (mean 15.35 Gb per individual, 51 M reads pair per sample) were produced and divided into 302,909,245 RAD tags (ranging from 1,934,524 to 6,800,866 with a mean of 3,188, 518.37) from the 95 offspring, for individual SNP discovery (Table S1).

SNP discovery and genotyping
RAD tags from each individual were clustered and compared, after stringent selection by the above-described method. A custom k-mer matching algorithm that excluded exact sequence matches (monomorphic loci) per 41 bp sequence was used to identify the SNP of the parents contained the RAD tags. A panel of 10,896 highfidelity SNPs in both parents was identified, and alleles for each marker were assigned to their respective parent donor. The SNPs were analyzed across the 95 offspring with genotyping data and classified into three categories using an in-house script: maternal heterozygous (3343 SNPs), paternal heterozygous (6066 SNPs), and heterozygous in both (1487 SNPs) (Table S2).

High-density linkage map construction
In the present study, Join-map and Lep-Map were used to construct the high-density genetics linkage map of A. davidianus. In both the maternal and paternal map, a total of 10,896 segregating SNPs were successfully classified into 30 linkage groups which was consistent with the karyotype 2n = 60 ( Fig. 1). The maternal map contained 4830 SNPs with a total genetic distance of 2580.04 cM., The length of LGs ranged from 0 cM (LG27) to 298.26 cM (LG1) with an average genetic length of 86.00 cM (Table 1). Seven thousand five hundred fifty-three SNPs representing a total length of 2011.80 cM ranging from 0 cM (LG21, LG27, and LG30) to 244.33 cM (LG1) with an average genetic length of 67.06 cM were consisted the corresponding paternal map. The longest linkage group was LG1 which was 244.33 cM and 298.26 cM in maternal and paternal map, respectively. However, the shortest was 0 cM distributed in LG21, LG27, and LG30 in the paternal map and LG27 in the maternal map.
Finally, the integrated map consisting of 30 LGs was constructed by 10,896 effective markers. The total length of the integrated map was 2643.10 cM. The average number of effective loci was 363 distributed in genetic distance of 88.10 cM with an average inter-locus distance of 0.24 cM. The genetic distance of LGs ranged from 17.61 cM (LG30) to 280.81 cM (LG1), and a mean inter-locus distance ranged from 0.11 (LG3) to 0.48 cM (LG26). The maximum density of SNP marker in linkage group was LG3 containing 878 effective loci and minimum density was LG30 containing 105 effective loci (Table S3).

Discussion
Genetic linkage mapping is an important tool in genomic and genetic studies [36,37]. Abundant molecular markers with an established order of arrangement provide an excellent resource for breeding research. Most amphibians possessed a large genome [38,39], which exhibited a challenge in constructing a high-density genetic linkage map. Research exploring SNPs in A. davidianus has been reported, but identified markers were not used to build a linkage map [40]. Sessions et al. assembled thirty linkage groups in maternal, paternal, and integrated maps, showing consistency with the karyotype [41]. Andrias davidianus possesses a genome of~50 G [12], hampering genomic and genetic study because of lack of large genome assembly technology and linkage map. The present study exhibited the first high-density genetic linkage map using SNPs of A. davidianus.
Thirty linkage groups were constructed with the average marker-interval of 0.24 cM, which was the highest density linkage map in amphibians. In common frog, 107 SSR markers was distributed in 15 linkage groups with total length 1698.8 cM [23]. In Lissotriton newts, total of 1146 markers was distributed in 12 linkage groups with an average marker interval of 1.29 cM [25]. The highest density linkage map of fish species was constructed for channel catfish, with the average marker interval of 0.22 cM [21]. Hence, the A. davidianus SNP linkage map displays high density even compared to fish species with relatively small genomes. In this study, we constructed 0 cM linkage groups in maternal and paternal map respectively. The main reason may be the   specificity of A. davidianus genome. Some types of markers with different genotype distributed on different location in one chromosome. However, the other types of markers obtained the same genotype from same individual distributing on same location in linkage group (0 cM). These results indicated the LD differ in maternal and paternal linkage groups. The high-density linkage map provided the means of fine mapping A. davidianus sex-related QTL. Three loci were identified on the maternal map LG6 (qs-1, qs-2) and LG11 (qs-3). The genetic distance of the QTL locus ranged from 0.1 cM to 3.5 cM suggested a narrow regions in the linkage group. This may provide valuable information for sex-related gene detection. While the same method was used to identify the loci on the maternal and paternal maps, all observed sex-related loci were on the maternal map. It is possible that some key sexrelated markers were homozygous on the paternal map. Sex chromosomes have a greater number of repeat sequences, making the assembly of sex chromosomes and the construction of a sex linkage group more difficult, and sex is readily influenced by factors other than chromosome [42]. Some sex-related QTLs have been detected in other species, seven in the half-smooth tongue sole with five sex-linked markers [43], and 175 sexrelated regions in turbot, most distributed among three LGs, most likely to be the sex chromosome [44]. In the present study, the major loci were detected only on linkage group LG6, providing evidence for the sex chromosome of A. davidianus. We analyzed sex related markers in this region and found correlation coefficient which ranged from 57.81 to 68.18% between phenotype and genotype in eight maternal SNP marker. Possibly due to low sequencing depth, some key SNPs in sex-related loci were not mapped on the linkage group. This region was considered as sex-related loci. So many SNPs in this region was related with sex too [17]. Verification could be carried out to confirm sex-linked genes and markers in A. davidianus through expanded single-linkage data in future. Further research requires expanded knowledge of the A. davidianus genome.

Conclusions
We constructed a high-density genetic linkage map using ddRAD sequencing technology to obtain SNP genotyping data from one full-sib family. A total of 10,896 markers were grouped and oriented into 30 linkage groups (LG), representing 30 chromosomes of  .hk_hk_96246) in the sex-related loci were significantly correlated with phenotype which indicated the four markers were sex-related markers. We constructed a High-density genetic linkage map and mapped 3 sex-related loci, which provides the linkage map is a useful tool for genomic study and provides a genetic basis for marker-assisted breeding of Chinese giant salamander.

Mapping population
Total of 97 apparently healthy A. davidianus including two parents (8 year old) and 95 salamander larvae (body weight 34.24 ± 12.13 g; body length 17.85 ± 1. 76 cm) from one F3 full-sib family were collected from Zhejiang Yongqiang Chinese Giant Salamander Limited Company (Jinhua, Zhejiang province, China). The salamanders were soaked in 4 L MS222 with the concentration of 0.5 g/L for 10mins to anaesthetize. After all the salamanders were unconscious, the vertebra was broken from the neck, then the skin of the neck was removed, and muscle tissue was collected. All operations were carried out as per Yangyze River Fisheries Research Institute Care Committee (no. 2013001). The muscle tissue was collected from each salamander using the tissue DNA extracting kit (Omega, USA) to extract genomic DNA following manufacturer's instruction. The concentration and quality of DNA were detected by Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and 1% agarose gel electrophoresis. Genetic sex of each offspring was identified using a sex-specific marker explored by our lab [32].

RAD library construction, sequencing, and genotyping
The RAD library construction, sample indexing, and pooling followed the natural population [43,45].  [46] had been used to cluster the trimmed reads by sequence similarity. A maximum base-pair mismatch of one was allowed in this step for the genetic mapping population. Stacks were used under default parameters to collapse the RAD-tags into clusters for SNP calling. High-Quality and Biallelic SNPs were identified under following criteria which were 1) > =6 and < =100 reads for each loci, otherwise this loci will be considered as missing data, 2) present in at least 60% (< 39 genotyping-deletion rate) of the individuals.

Linkage map construction
RAD-based SNPs were first tested against the expected segregation ratio at the beginning of the linkage analysis. Markers were filtered out by a chi-square test for separation ratio at the intersection sites with P < 0.001. The linkage map was constructed using joinmap and Lep-Map. Joinmap was used to divide markers into groups using a grouping tree function in CP model with LOD = 8.0 [47] and Lep-MAP was used to determine marker order [48]. Kosambi mapping function was used to estimate of each linkage group distances.

Sex-related loci mapping
CIM method of WinQTL Cart 2.5 software was used to map the sex-related loci on the male and female linkage maps respectively [49]. Model six with three parameters as 10-cM window size, five control markers, and a step size of 1-cM was used to analyze the forward and backward stepwise regression. The location of each QTL was determined based on its LOD peak location and surrounding region with the threshold value of 3.0. Singlemarker association tests with WinQTL Cart 2.5 was used to analysis the association between phenotype and genotype for each marker [49]. All maternal SNPs in the located sex-related region were analyzed to detect correlation with the phenotype to find sex-related SNPs.