QTL and transcriptomic analyses for ammonia tolerance in the Pacific white shrimp ( )

Background: Ammonia is one of the most common toxicological environment factors affecting shrimp health. Although ammonia tolerance in shrimp is closely related to successful industrial production, few genetic studies of this trait are available. Results: In this study, we constructed a high-density genetic map of the Pacific white shrimp ( Litopenaeus vannamei ) using specific length amplified fragment sequencing (SLAF-seq). The constructed genetic map contained 807,505 polymorphic markers spanning 44 linkage groups, with a total distance of 6,360.12 centimorgans (cM) and an average distance of 0.37 cM. Using this genetic map, we identified a quantitative trait locus (QTL) that explained 7.41–8.46% of the phenotypic variance in L. vannamei survival time under acute ammonia stress. We then sequenced the transcriptomes of the most ammonia-tolerant and the most ammonia-sensitive individuals from each of four genetically distinct L. vannamei families. We found that 7546 genes were differentially expressed between the ammonia-tolerant and ammonia-sensitive individuals. Using QTL analysis and the transcriptomes, we identified one candidate gene (annotated as an ATP synthase g subunit) associated with ammonia tolerance. Conclusions: In this study, we constructed a high-density genetic map of L. vannamei and identified a QTL for ammonia tolerance. By combining QTL and transcriptome analyses, we identified a candidate gene associated with ammonia tolerance. Our work provides the basis for future genetic studies focused on molecular marker-assisted selective breeding.

nucleotide polymorphisms (SNPs) associated with ammonia tolerance in L. vannamei using marker-trait correlation analyses [15], while Lu  Finally, several studies identified transcriptomic changes and differentially expressed genes in L. vannamei after ammonia stress [9,18]. However, no studies have investigated the QTLs associated with ammonia tolerance in shrimp.
QTL analysis effectively identifies molecular markers or candidate genes associated with economically-important traits in plants and animals [19]. QTL analyses usually require high-density genetic linkage maps. To date, genetic linkage maps have primarily been constructed using high-throughput sequencing technologies, such as restriction site-related DNA sequencing (RAD-seq), genotyping sequencing (GBS), and specific length amplified fragment sequencing (SLAF-seq) [20]. In particular, SLAF-seq efficiently identifies and genotypes large-scale SNPs [20]. SLAF [29]. SLAF-seq has also been successfully applied to L. vannamei [ 30].
Therefore, we used SLAF-seq to construct a high-density genetic map of L. vannamei. We then performed a QTL analysis of ammonia tolerance in L. vannamei.
We also compared the transcriptomic differences between ammonia-tolerant and ammonia-sensitive individuals across several L. vannamei families to identify candidate genes within QTLs.

Preparation of the mapping family
The L. vannamei used in all experiments were obtained from the shrimp-breeding center at the Guangxi Academy of Fishery Sciences (Nanning, Guangxi, China). The L. vannamei family used for mapping was constructed using artificial insemination.
In brief, a male shrimp from a family with a relatively high ammonia-tolerance (obtained via 10 consecutive generations of breeding) was mated with a female shrimp from a common family. The hatched offspring were reared for about 1 year.
Then, a male and female shrimp were randomly selected from the year-old offspring and mated. The F1 progeny were used for mapping (LV-N).

Measurement of ammonia tolerance
We randomly selected 284 shrimp (average body weight: 20.78 g) from the LV-N family. Selected shrimp were transferred to a 2 m × 4 m × 1 m indoor pool and allowed to acclimate for one week. Aquatic conditions during the acclimation and experimental periods were kept constant: temperature of 27.0 ± 0.5°C, pH of 8.1 ± 0.2, salinity of 30.2‰, and dissolved oxygen of 6-8 mg/L; culture water was kept aerated, and shrimp were fed formulated pellets (Zhengda Feed, China) daily at a ratio of 5% body weight. Following acclimation, we performed an acute ammonia stress test. The ammonia-N concentration used for the acute stress test was 345.94 mg/L, based on the results of a preliminary experiment. This was the concentration at which half of the experimental shrimp died in 72 hours under stress. The ammonia-N concentration of the water in the experimental pool was controlled by adding NH 4 Cl stock solution (prepared by dissolving analytically pure NH 4 Cl in filtered seawater). The concentration of ammonia-N in the water was measured daily using standard methods [31]. To keep the ammonia-N concentration constant, NH 4 Cl stock solution was added if the ammonia-N concentration was <345.94 mg/L, and seawater was added if the ammonia-N concentration was >345.94 mg/L. During the experiment, shrimp heath was observed once per hour, and dead shrimp were removed immediately. Shrimp were considered dead when lying motionless on the bottom of the pool and not responding to external stimuli. Collected dead shrimp were immediately frozen in liquid nitrogen and stored at −20°C for DNA extraction.
The survival time of each shrimp was used as a proxy for ammonia tolerance. The experiment ended when all shrimp had died.

DNA extraction
DNA was collected from the 284 F1 (LV-N) shrimp and the two parent shrimp. Marine animal genomic DNA extraction kits (Tiangen Biotech, China) were used to extract DNA from the tail muscle of each shrimp. DNA was quantified using a NanoDrop spectrophotometer and 1% agarose gel electrophoresis with a lambda DNA standard.

SLAF library preparation and sequencing
First, we predicted the digestion of the L. vannamei genome (https://www.ncbi.nlm.nih.gov/genome/?term=Vannamei) [32] using self-developed software. We digested the extracted genomic DNA of all LV-N shrimp using the endonucleases identified by the predictive software. Then, dual-index sequencing adaptors were ligated to the DNA fragments obtained by digestion with T4 ligase, and the fragments were amplified using polymerase chain reactions (PCRs). PCR products (314-414 bp including the adaptor sequences) were purified and reamplified using PCR. SLAF sequencing was carried out on an Illumina HiSeq system, following the Illumina-recommended procedure. To assess the accuracy of library construction, we performed the same library-construction and sequencing steps using the genome of Oryza sativa japonica as a control. Library construction and sequencing were performed by Biomarker Technologies Corporation (Beijing, China).

SLAF-seq data analysis and genotyping
The raw sequencing reads were quality controlled by removing reads with a quality score <20. The remaining raw reads were grouped by individual based on the dualindex adaptor sequences. The dual-index adaptor and 5-bp end sequences were then trimmed to obtain clean reads. The clean reads were mapped to the L.
vannamei genome (https://www.ncbi.nlm.nih.gov/genome/?term=Vannamei) [32] using BWA [33]. Reads mapped to the same position with >95% identity were considered the same SLAF. SNP-based polymorphic SLAF markers were identified by aligning reads from the same SLAF sequence. These polymorphic SLAF markers were then filtered by removing those with a parental sequencing depth less than 10-fold; those where the number of SNPs was >5; those where the proportion of genotypes covering offspring was <70%; and those with significant segregation distortion (chisquare test P<0.05). The remaining polymorphic SLAFs were classified into eight separate patterns: aa × bb, ab × cd, cc × ab, ab × cc, ef × eg, hk × hk, nn × np, and lm × ll. Because the mapping population used in this study was an F1 population, the polymorphic SLAF with the pattern aa × bb was removed, and the remaining polymorphic SLAFs were used for the construction of the genetic map.

Genetic map construction and QTL analysis
After coding the genotypes of the polymorphic SLAF markers, the genetic map was constructed using the single-chain clustering algorithm in HighMap [34], with the probability log threshold set to ≥5.0 and a maximum recombination rate of 0.4. The Kosanbi mapping function was used to convert percent recombination to genetic distance (cM). QTL analysis was conducted using the R/qtl software package [35].
The logarithm of odds (LOD) threshold was determined based on 1,000 permutations (P < 0.05). The phenotypic variance explained by the QTL was estimated using the formula 1 - 10 -2LOD/n , where n was the sample size [36].  Table S1). We randomly selected 200 shrimp from each family, and subjected all shrimp to the acute ammonia stress test (345.94 mg/L ammonia-N), as described above. In each family, the 20 shrimp with the longest survival times (i.e., the most ammonia tolerant) were collected, as were the 20 shrimp with the shortest survival times (i.e., the most ammonia sensitive). When collecting the ammonia-sensitive shrimp, shrimp that were out of balance and lying on the bottom of the pool were judged to be dying, and were collected immediately, without waiting for death. The hepatopancreas of each shrimp was extracted, and hepatopancreases were pooled to form one ammonia-tolerant sample and one ammonia-sensitive sample per family.
Total RNA was extracted from each pooled sample using TRIzol reagent (Invitrogen, USA), following the manufacturer's instructions. Residual genomic DNA was removed with DNase I. RNA purity (OD260 / 280), concentration, and absorption peak were measured using a NanoDrop 2000. RNA integrity was assessed using an RNA Nano 6000 Assay Kit with an Agilent Bioanalyzer 2100. The isolated mRNA was divided into 100-400 bp fragments using an RNA fragment reagent (Illumina, USA). cDNA libraries were then constructed using NEBNext Ultra RNA Library Prep Kits (Illumina, USA), following manufacturer's recommendations, and sequenced on an Illumina HiSeq system (Illumina, USA). Library construction and sequencing were performed by Biomarker Technologies Corporation (Beijing, China).
Raw sequencing reads were trimmed and filtered using in-house Perl scripts to remove adaptor sequences and low-quality reads; the Q20, Q30, GC-content, and sequence duplication levels of the clean data were calculated. Clean reads were then aligned to the L. vannamei genome (https://www.ncbi.nlm.nih.gov/genome/? term=Vannamei) [32] using Hisat2 2.1.0 (http://ccb.jhu.edu/software/hisat2/index.shtml) [37]. Matched reads were counted to determine gene expression levels using the fragments per kilobase of transcript per million mapped reads (FPKM) method [38]. DEGs were identified using edger [39]. We considered unigenes differentially expressed when the false discovery rate (FDR) was ≤ 0.01 and the fold change between groups was > 2. DEGs were functionally annotated against the following databases: Non-Redundant protein After obtaining DEGs, we identified candidate genes among the DEGs. We consider candidate genes associated with ammonia tolerance when (1) candidate genes located within the QTL interval; (2)  Relative gene expression levels were calculated using the 2-ΔΔCT method [42].

Phenotypic variation
We developed an ammonia-tolerant shrimp family (designated LV-N) for mapping, and subjected 284 LV-N shrimp to an acute ammonia stress test. All shrimp died within 2-98 hours, with a mean survival time of 65 hours. Individual survival times were normally distributed and thus suitable for QTL detection.

SLAF-seq and genotyping
Based on our digestive enzyme prediction using the reference genome of L. Across all reads, the average Q30 was 95.81%, the average GC content was 40.60%, and the GC distribution was normal (Table 1). We used the rice (Oryza sativa japonica) genome as a control to estimate the validity of our library construction.  Table S2). These results indicated that SLAF library construction and sequencing were normal.
After filtering and clustering all reads, 807,505 SLAFs were identified. The average sequencing depth of these SLAFs was 42.8-fold for the male parent, 42.14-fold for the female parent, and 12.43-fold for the progeny (Table 1). Of the 807,505 highquality SLAFs detected, 293,415 (36.34%) were polymorphic (Table 1). After further filtering, the remaining 115,973 SLAF markers were successfully classified into eight genotypic patterns: ab × cd, cc × ab, aa × bb, ab × cc, ef × eg, lm × ll, hk × hk, and nn × np. The most common pattern was aa × bb, followed by nn × np and lm × ll (Fig. 1). Because the mapped population was an F1 population, we eliminated aa × bb as a valid marker.

Characteristics of the genetic map
Linkage analysis labeled 17,338 SLAF markers on the genetic map: 11,512 on the male map, 10,293 on the female map, and 17,338 on the sex-average map (Fig. 2).  Table S3, Table S4, and Table S5

QTL mapping of ammonia-tolerance
A QTL analysis of the ammonia-tolerance trait in the LV-N L. vannamei family was performed based on the genetic maps. The LOD threshold was 4.75 (1000 permutations, P < 0.05). Thus, QTLs with LOD scores >4.75 were considered effective QTLs. Using this criterion, we identified a QTL within LG19 for ammonia tolerance (Fig. 3). The phenotypic variation explained by this QTL was 7.41-8.46%, the LOD score was 4.75-5.45, and the confidence interval was 12.42-29.43 cM.  Table S6).

Transcriptome sequencing, candidate gene identification and qRT-PCR
The numbers of DEG annotations recovered in the databases searched were similar across the four L. vannamei families. For instance, the COG terms mainly enriched in the DEGs from all four families were posttranslational modification, protein turnover, chaperones, and general function prediction only (Fig. 4); the GO terms primarily enriched in the DEGs from all four families were binding, catalytic activity, cellular process, metabolic process, cell, cell part, single-organism process and membrane functions (Fig. 5).
By aligning the DEGs with the QTL region in LG19, we identified 107 DEGs located in the QTL interval. The expression levels and annotations of these DEGs were listed in Supplementary Material Table S7. Of these DEGs, only one gene (LOC113809108) met the criterion we used to determine candidate genes associated with ammonia tolerance. This gene was annotated as an ATP synthase g subunit. LOC113809108 was located in the QTL interval, and was significantly upregulated in the most ammonia-tolerant shrimp compared to the most ammonia-sensitive shrimp from families LV-N and LV-C (Fig. 6). This gene was also upregulated in the most ammonia-tolerant shrimp from families LV-A and LV-F, but this difference in expression was not significant ( Table 2).
The qRT-PCR analysis showed that the patterns of LOC113809108 gene expression in ammonia-tolerant and ammonia-sensitive pooled samples from the families LV-A, LV-C, LV-F, and LV-N were similar to the patterns determined using RNA-seq: LOC113809108 gene expression was upregulated in the ammonia-tolerant shrimp as compared to the ammonia-sensitive shrimp across all four families (Fig. 6).

Discussion
The genome of L. vannamei is large (~2.45 Gb) [32]. Whole-genome deep resequencing is relatively costly for large genomes and is often not necessary for gene/QTL mapping [43,44]. In this study, we constructed a high-density genetic map of L. vannamei using SLAF-seq, which is an effective method by which to discover large numbers of SNPs and to perform large-scale genotyping [20].
Compared to traditional methods of genetic map construction (e.g., random amplified polymorphic DNA (RAPD), amplified fragment length polymorphism (AFLP), and simple sequence repeat (SSR), the SLAF-seq method has several advantages for large-scale SNP discovery and genotyping: high density, high throughput, high efficiency, and low cost [45]. Previous studies have developed genetic maps of L. vannamei using RAPD, AFLP, and SSR, but in these maps, the average distance between adjacent markers was 1-5 cM [46][47][48][49]. The average distance between adjacent markers in the SLAF-seq genetic maps of L. vannamei in this study was substantially shorter (0.34 cM). Notably, the average distance between adjacent markers found here was also less than in a previously reported SLAF-seq genetic map of L. vannamei (0.75 cM) [30], possibly because we used a larger sample size and a greater sequencing depth. However, the number of LGs in the genetic map of L. vannamei in this study (44) was consistent with the number of LGs in the previously reported genetic map of L. vannamei [ 30]. This indicated that L. vannamei had 44 chromosomes.
As most animals cannot self-fertilize, it is difficult to develop common populations for genetic mapping (e.g., F2, recombinant inbred line (RIL), and nearlyisogenic line (NIL) populations). Therefore, we used an F1 population of L. vannamei to construct the genetic map, relying on a pseudo-testcross strategy. This strategy was based on the selection of single-dose markers present in one parent and absent in the other, and carried at a 1:1 ratio by the F1 offspring [50]. Therefore, gamete separation in each individual can be directly analyzed. The pseudo-testcross strategy has been widely used to construct animal F1 populations for genetic mapping [51][52][53][54]. In this study, we developed an F1 population using one ammonia-tolerant male parent (the result of 10 generations of selective breeding) and one female shrimp from a common family.
Previous studies have suggested that the size of the mapped population might affect the accuracy of the genetic map and the QTL analysis, and have shown that genetic map accuracy increases with the size of the population used [55].
Specifically, populations of >200 individuals are considered sufficient for the construction of accurate genetic maps [55]. We thus used 284 randomly-selected individuals from the F1 population to construct the linkage map. However, the determination of shrimp survival time during the ammonia stress experiment depended on human observation, and thus may not have been perfectly accurate.
To reduce the possible impacts of measurement inaccuracies on the QTL analysis, we used a relatively large F1 population. This larger population increased the accuracy of the QTL mapping, compensating for any instances of human error in survival time measurement.
Ammonia stress is one of the biggest challenges facing shrimp aquaculture.
Ammonia not only has a direct lethal effect on shrimp [7], but also inhibits the shrimp immune system and increases sensitivity to pathogens [56]. Breeding new varieties of ammonia-tolerant shrimp is therefore an important target of the shrimp aquaculture industry. Here, we identified a QTL for ammonia-tolerance, located on LG19 at 169.09-169.49 cM, that explained 7.41-8.46% of the phenotypic variation in ammonia tolerance. To the best of our knowledge, this is the first QTL for ammoniatolerance reported in shrimp. However, this QTL spans a large chromosomal region and may contain hundreds of genes. Therefore, in order to identify functional genes associated with ammonia tolerance, we sequenced the transcriptomes of four L.
vannamei families (LV-N, LV-A, LV-C, and LV-F) with different genetic backgrounds.
Combining QTL mapping and gene expression analysis, we identified a single DEG (LOC113809108), located in the QTL interval, that was annotated as an ATP synthase g subunit and was significantly upregulated in the ammonia-tolerant LV-N and LV-C shrimp. ATP synthase is a double-motor enzyme that is involved in ATP synthesis, ATP hydrolysis-dependent processes, and the regulation of the proton gradient across some membrane-dependent systems [57]. Several studies have shown that ammonia excretion in aquatic animals is associated with Na + /K + -ATPase, which is mainly located on the basolateral membrane of branchial cells; NH 4 + is excreted into the environment when K + is replaced by NH4 + via the Na + /NH4 + exchanger [58][59][60]. Indeed, a previous study suggested that in L. vannamei, high ammonia tolerance was mainly a result of improved ammonia excretion and detoxification, as well as an accelerated energy metabolism [17]. Therefore, we speculate that ATP synthesis might affect the ammonia tolerance of L. vannamei by regulating ATP synthesis and controlling cellular ammonia excretion.

Conclusions
In this study, we constructed a high-density genetic map of L. vannamei and identified a QTL for ammonia tolerance. By combining QTL and transcriptome analyses, we identified a candidate gene associated with ammonia tolerance. Our work provides the basis for future genetic studies focused on molecular markerassisted selective breeding.     Expression of LOC113809108 gene from the transcriptomic analysis validated by qRT-PCR. Ex

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download. Supplementary