Juveniles from 140 shrimp families were evaluated for their tolerance to high-salinity stress. Briefly, a total of 100 juveniles from each family were raised in separate tanks with seawater at a salinity of 25 ppt. In order to reduce the tank effect during the experiment, all the families were cultured in the same room, and the culturing seawater in each tank was the same. Moreover, the increase of the salinity in each tank was checked every 6 h per day. The salinity of seawater in each tank was raised from 25 to 55 ppt gradually over 3 days by adding seawater with a salinity of 100 ppt. Briefly, the seawater salinity was adjusted 4 times at 6-h intervals each day, and around 2.5 ppt was raised each time. Dead shrimp were collected every 4 h. The final survival rates of these 140 families after high-salinity treatment were calculated and presented in Supplementary Figure 1 . Based on the survival rate of each family during high-salinity stress, two families (fam1 and fam2) were chosen as tolerant families (designated as tFam1 and tFam2, respectively), and another two families (fam103 and fam139) were selected as susceptible families (designated as sFam1 and sFam2, respectively). More information about these four families is shown in Table 1 . In order to exclude the effects of shrimp size, two paired comparison groups (Groups I and II) were chosen for further experiments and sampling. Group I included tFam1 and sFam1, and Group II included tFam2 and sFam2. Shrimp in each group had similar body lengths but different tolerance.

Another 60 healthy juveniles from above four selected families cultured under normal salinity of 25 ppt were collected and put into four separate tanks. They were kept in seawater with a salinity of 25 ppt at a temperature of 26°C ± 1°C for 12 h before further high-salinity treatment. The cephalothoraces of nine individuals were collected from each family as normal salinity samples. The entire cephalothoraces of three individuals were put together as one sample, and three samples were taken for each family and preserved in liquid nitrogen for further RNA extraction. The samples of tolerant families (tFam1 and tFam2) were named tFam1-25 and tFam2-25, and the samples of susceptible families (sFam1 and sFam2) were named sFam1-25 and sFam2-25. After the above samples were collected, seawater with a salinity of 100 ppt was added into each tank gradually to raise the salinity of seawater to 45 ppt during the next 48 h. The procedure to raise the salinity of seawater was the same as described in Section 2.1. When the salinity of seawater reached 45 ppt, shrimp were collected from each family. The sampling method was the same as described above. The samples were designated as tFam1–45, tFam2–45, sFam1–45, and sFam2–45 and preserved in liquid nitrogen for further use.

Total RNA was extracted using RNAiso reagent (Takara, Maebashi, Japan) following the manufacturer’s protocol. The quality and integrity of RNA were evaluated by electrophoresis on 1% agarose gel. RNA content was quantified by NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The cDNA was synthesized by PrimeScript™ RT reagent Kit with gDNA Eraser kit (Takara, Japan) following the manufacturer’s protocol.

The mRNA was enriched from the total RNA using Oligo (dT) beads and Ribo-ZeroTM Magnetic Kit (Epicentre, Madison, WI, USA); then the enriched mRNA was fragmented into short fragments, which were reverse transcribed into cDNA. After purification with QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end reparation, and poly(A) tailing, the cDNA fragments were ligated to Illumina sequencing adapters. The ligation products at suitable sizes were sequenced using Illumina HiSeq2500 by Gene Denovo Biotechnology Co. (Guangzhou, China). The raw reads from the sequencing platform were filtered by fastp (Chen et al., 2018) with the following parameters: Q quality score ≤ 20, reads with undetermined nucleotides larger than 10%, and reads containing adapters. The rRNA-mapped reads were removed, and the remaining reads were mapped to the genome sequence of L. vannamei (http://www.shrimpbase.net/lva.download.html) using Bowtie 2 (Langmead and Salzberg, 2012) and TopHat2.2.4 (Kim et al., 2015). The mapped reads were assembled with StringTie v1.3.1, and the gene expression levels were calculated using the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) method (Pertea et al., 2015; Pertea et al., 2016).

Differential expression analyses of unigenes were carried out using DESeq2 (Love et al., 2014) software. The differentially expressed genes (DEGs) were selected by the false discovery rate (FDR) <0.05 and the expression fold change (FC) >2 (log₂|FC| > 1). Firstly, the transcriptional differences between tolerant and susceptible families were focused on; therefore, the DEGs in comparison groups (Group I and II) under both normal and high salinity were identified. Secondly, the response of shrimp to salinity change at the transcriptional level was analyzed, and the DEGs for each family were analyzed when they were cultured under normal or high-salinity seawater. In order to reduce the false positives, the DEGs with similar expression patterns in two comparison groups or similar tolerance families were selected for further analysis.

All selected DEGs were used for Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Briefly, DEGs were mapped to GO terms in the GO database (http://geneontology.org) and pathways in the KEGG database (http://www.genome.jp/kegg/). Then, significantly enriched GO terms and KEGG pathways were recognized by a hypergeometric test-taking FDR ≤ 0.05 as a threshold.

The samples used for qPCR detection were the same as those used for transcriptome sequencing. Four DEGs were quantified with qPCR to validate the transcriptome analysis result. The Primer3plus (http://primer3plus.com/cgi-bin/dev/primer3plus.cgi) was used for the design of gene-specific primers based on their sequences. Primers used for qPCR detection are shown in Supplementary Table 1 , and 18S rRNA was used as the reference gene. The qPCR was conducted as follows: denaturation at 94°C for 2 min, 40 cycles of 94°C for 30 s and annealing at temperature 58°C for 20 s, and a final extension of 30 s at 72°C. The melting curve was used to check the specificity of primers, and the relative expression level was calculated with the 2^−ΔΔCt method (Livak and Schmittgen, 2001). The two-way ANOVA (factor: salinity and family) was utilized to determine the statistically significant differences between different factors, and a one-way ANOVA was used to reveal the significant differences among different samples. These statistical analyses were performed in R.

A total of 24 RNA-seq libraries were constructed for two tolerant families and two susceptible families cultured under two salinity conditions. All transcriptome sequencing data were uploaded to National Center for Biotechnology Information (NCBI) with the accession number PRJNA799461. The sequencing information was summarized in Supplementary Table 2 . After filtering, 99.68%–99.85% of raw reads were taken as clean reads, and 85.82%–89.70% of the clean reads were mapped to the reference genome of L. vannamei (Zhang et al., 2019). A total of 31,753 unigenes were annotated, and each library contained 60.74%–78.57% of total genes.

A total of 963 DEGs were identified in Group I between tFam1-25 and sFam1-25, while 1,242 DEGs were identified in Group II between tFam2-25 and sFam2-25. In order to reduce the false-positive results caused by the family itself, the DEGs with similar expression patterns within two comparison groups were chosen for further analysis. A total of 98 DEGs presented similar expression profiles in these two groups. Among them, 58 DEGs were upregulated and 40 DEGs were downregulated in two tolerant families ( Figure 1A ). GO enrichment analysis on these DEGs showed that they were enriched in metabolic process, single-organism process, cellular process and response to a stimulus for the biological process, membrane and membrane part for cellular component, and binding and catalytic activity for molecular function ( Figure 1B ). Further analysis of the GO terms for the biological process showed that DEGs enriched in response to stimulus were upregulated in tolerant families ( Figure 1C ). Nine genes encoding crustacyanins were found among these DEGs. In contrast, there were no significant enriched GO terms for DEGs, which were highly expressed in susceptible families ( Figure 1D ).

A total of 865 DEGs were identified in the comparison group tFam1-45 and sFam1-45 and 1,660 DEGs in the other comparison group tFam2 and sFam2-45. In order to reduce the false positive caused by the family itself, 56 DEGs with the same expression pattern in the two comparison groups were chosen for further analysis. Among these DEGs, 30 DEGs showed higher expression levels in tolerant families, and 26 DEGs presented higher expressions in susceptible families ( Figure 2A ). Further analysis of these DEGs showed that metabolic process, cellular process, single-organism process, and response to stimulus were the top four enriched GO terms in the category of the biological process. Binding and catalytic activity were the top two enriched GO terms in the category of molecular functions, which were very similar to those under normal salinity ( Figure 2B ). Statistical analysis showed that genes highly expressed in tolerant families were enriched significantly (Q value <0.05) in terms of response to a stimulus for the biological process. Two genes encoding crustacyanins were found in the enriched response to the stimulus GO term. As for the DEGs upregulated in susceptible families, no significant enrichment term was found ( Figures 2C, D ).

Further analysis of the DEGs between tolerant and susceptible families under both normal and high salinities showed that four DEGs had higher expression levels in tolerant families, and six DEGs presented higher expression levels in susceptible families ( Table 2 ). Four highly expressed genes in tolerant families were crustacyanin subunit C, 26S protease regulatory subunit 6A-B, actin, and toll-like receptor. Among the six DEGs highly expressed in susceptible families, four genes had annotations, including trypsin-1-like, compound eye opsin BCRH2-like, triose-phosphate isomerase, and mucin-2-like.

A total of 891, 409, 1,649, and 1,149 DEGs were identified from tFam1, sFam1, tFam2, and sFam2 families, respectively, after high-salinity treatment. Further analysis showed that 26 DEGs presented similar expression patterns in these four families. Among them, 21 DEGs showed upregulated expression levels and five DEGs showed downregulated expression levels after high-salinity treatment ( Figure 3A ). GO enrichment analysis on these DEGs showed that they were enriched in the items including single-organism process, metabolic process, cellular process for biological process, membrane and membrane part for cellular components, catalytic activity, and binding for molecular functions ( Figure 3B ). Further analysis showed that oxidation–reduction process was the most enriched item for the biological process ( Figure 3C ). Among the 21 upregulated DEGs, six unigenes were related to the oxidation–reduction process, including cytochrome P450 epoxidase-like protein, d-aspartate oxidase, alkylglycerol monooxygenase-like, spermine oxidase-like isoform X2, type I iodothyronine deiodinase-like isoform X1, and la-related protein 6-like.

In order to elucidate how the tolerant family adapts to high-salinity treatment, the DEGs for the two tolerant families under normal and high salinity were further analyzed. A total of 233 DEGs presented similar expression profiles in these two families ( Figure 4A ). Among these DEGs, 160 DEGs were upregulated and 73 DEGs were downregulated after high-salinity treatment in comparison with their controls at a salinity of 25 ppt ( Figure 4B ). GO enrichment analysis indicated that metabolic process was the most enriched GO term for the biological process, membrane was the most enriched one for cellular component, and catalytic activity was the most enriched GO term for molecular function ( Figure 4C ). Among them, serine-type endopeptidase activity, serine-type peptidase activity, and serine hydrolase activity were the top three enriched GO terms ( Figure 4D ). Amino acid metabolic process including tyrosine metabolism pathway and beta-alanine metabolism pathway were the top two KEGG enrichment pathways ( Figure 4E ). Genes encoding serine proteases, prophenoloxidase, and some osmoregulation genes, such as solute carriers and sodium/glucose cotransporter, showed significant upregulations in these two tolerant families after high-salinity treatment. These data indicated that the amino acid metabolism pathway and osmoregulation process might play a key role in shrimp to handle with high-salinity stress. Interestingly, five rhodopsin genes were found to be upregulated in tolerant families under high-salinity treatment, and they were also DEGs between tolerant families and susceptible families under high-salinity conditions.

In order to know whether the DEGs shared by two tolerant families showed different expression patterns in susceptible families after high-salinity treatment, the expressions of those DEGs in susceptible families were further analyzed. A total of 63 upregulated DEGs and 45 downregulated DEGs were found to be specifically recognized in tolerant families after high-salinity treatment ( Tables 3 , 4 ). Serine-type peptidase activity, serine hydrolase activity, and serine-type endopeptidase activity were the top three enriched GO terms in the DEGs specifically upregulated in the tolerant families. Four genes belonging to serine proteases families, including clip domain serine proteases, serine proteases, chymotrypsin-like protein, and prophenoloxide-activating enzyme III, were shared by the above three GO terms ( Figures 5A, B ), which showed higher expressions in tolerant families after high-salinity treatment. These data suggested that serine protease activity should play key roles in shrimp tolerance to high salinity. As for the DEGs specifically downregulated in tolerant families, no significant enriched GO term was found ( Figures 5C, D ).

The qPCR detection of four representative DEGs, including LVAN11695, LVAN07881, LVAN23509, and LVAN04611, are shown in Figure 6 . The expression tendency of these four genes was consistent with the transcriptomic sequencing results, suggesting that the transcriptome analysis data were reliable.