Preparation of the mapping family
The L. vannamei used in 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
A total of 284 shrimp (average body weight: 20.78 g) were randomly selected 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, an acute ammonia stress test was performed. 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 NH4Cl stock solution (prepared by dissolving analytically pure NH4Cl 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, NH4Cl 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 every 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 re-amplified using PCR. SLAF sequencing was carried out on an Illumina HiSeq system, following the Illumina-recommended procedure. To assess the accuracy of library construction, the same library-construction and sequencing steps using the genome of Oryza sativa japonica as a control was performed. 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 dual-index 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 (chi-square 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].
Transcriptome sequencing, candidate gene identification and quantitative real-time PCR (qRT-PCR) verification
To identify differentially expressed genes (DEGs) between ammonia-tolerant and ammonia-sensitive L. vannamei, the transcriptomes of four L. vannamei families were sequenced: the mapping family (LV-N) and three other randomly chosen common families (LV-A, LV-C, and LV-F) with different genetic backgrounds. Our previous analysis indicated that the 24-h median lethal concentration of NH4Cl was 140.96 mg/L, 189.19 mg/L, 117.88 mg/L, and 137.26 mg/L for families LV-A, LV-C, LV-F and LV-N, respectively (Supplementary Material, Table S1). Two hundred shrimp from each family were randomly selected, and subjected to the acute ammonia stress test (345.94 mg/L ammonia-N), as described above. In each family, 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, specimens 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 an ammonia-tolerant sample and an 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]. unigenes were considered 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 sequences (NR) (ftp://ftp.ncbi.nih.gov/blast/db/), Protein family (Pfam) (http://pfam.xfam.org/), Clusters of Orthologous Groups (http://www.ncbi.nlm.nih.gov/COG/), Swiss-Prot (http://www.uniprot.org/), Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/), and Gene Ontology (GO) (http://www.geneontology.org/).
After obtaining DEGs, candidate genes among the DEGs were identified. We consider candidate genes associated with ammonia tolerance when (1) candidate genes located within the QTL interval; (2) candidate genes differentially expressed between the most ammonia-tolerant and the most ammonia-sensitive individuals in the mapping family (LV-N); (3) the regulation pattern (up- or down-regulated expression) of candidate genes between the most ammonia-tolerant and the most ammonia-sensitive individuals was consistent across the four families of Litopenaeus vannamei (LV-A, LV-C, LV-F, and LV-N).
qRT-PCR was used to validate the RNA-seq results by quantifying the expression of the candidate gene (LOC113809108) in the ammonia-tolerant and ammonia-sensitive pooled samples from the four L. vannamei families (LV-A, LV-C, LV-F, and LV-N). RNA-seq and qRT-PCRs analyses were carried out using the same samples. qRT-PCRs were performed using SYBR Premix Ex TaqTM II kits (TaKaRa, Japan), according to the manufacturer's instructions. The primer sets used to detect LOC113809108 gene expression levels were designed using the Primer Premier software (version 5.0) [40] as follows: 5’-ACTTGGGTGCTGTAGCTCAA-3’ and 5’-CTCGACAGCAACCAGGGTAT-3’. L. vannamei 18S RNA was used as the internal reference gene; this gene was amplified using the primer sets as follows: 5’-GCCTGAGAAACGGCTACCACATC-3’ and 5’-GTAGTAGCGACGGGCGGTGTGT-3’ [41]. The qRT-PCR cycling program was as follows: preheating at 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. The qRT-PCR was carried out at 95 °C for 40 s, 95 °C for 5 s, and 62 °C for 30 s for 40 cycles. Three parallel qRT-PCRs were carried out for each sample. Relative gene expression levels were calculated using the 2-ΔΔCT method[42].