Female Pacific white shrimp (L. vannamei) were cultured at Jinyang Biotechnology Co. Ltd., Maoming, Guangdong, China. The shrimp were maintained in artificial seawater [pH 8.2 and 30 parts per thousand (ppt)] at 28°C and fed every morning with a normal shrimp diet. Female shrimp grown until body weights of 47.2±1.5 g and body lengths of 18.1±0.3 cm were used for analysis of ovarian development. At 7 days prior to artificially induced maturation, the shrimp were fed fresh clam worms, squid and oysters for nutritional reinforcement (Figure 1A). The animal experiments were conducted in accordance with the guidelines and approval of the Ethics Committees of South China Sea Institute of Oceanology, Chinese Academy of Sciences.

Ovaries and hepatopancreas samples from female shrimp were collected from 1 to 6 days after unilateral eyestalk ablation for stages I, II, III and IV, respectively (Figure 1B). Ovarian maturation stages (stage I–IV) were assigned after histological processing of the anterior end of ovarian explants with the number 2–9 lobes (Liaaen-Jensen, 1997; Bae et al., 2017). Stage I ovaries are mainly comprised of oogonia and previtellogenic oocytes; stage II ovaries primarily comprise previtellogenic oocytes and endogenous vitellogenic oocytes; stage III ovaries process a large number of oocytes developing from endogenous vitellogenetic oocytes to exogenous vitellogenetic oocytes; and exogenous vitellogenic oocytes and mature oocytes predominate in stage IV ovaries (Figure 1D). The gonadosomatic index (GSI) of each shrimp was calculated as a percentage of ovary weight to body weight, and the GSIs for shrimp in stages I, II, III and IV were 0.5–1.0%, 1.0–2.5%, 2.5–5.0% and >5.0%, respectively (Figure 1C). All samples were frozen immediately in liquid nitrogen after dissection, and total RNA was extracted by using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA).

The RNA integrity value (RIN) was checked by an Agilent 2100 Bianalyzer (Agilent, USA), and RNA-seq libraries were prepared following the standard Illumina protocols by the Beijing Genomics Institute (Beijing Genomics Institute, Beijing, China). In brief, mRNA equally contributed from 3 individuals was enriched by oligo-dT selection. Double-stranded cDNAs were synthesized using the fragmented mRNA as templates with random hexamer primers. The short fragments were connected with adapters, and suitable fragments were selected for PCR amplification. Then, the library was sequenced using an Illumina HiSeq 6000 (Illumina Inc., San Diego, CA, USA) in BGI.

The raw sequencing data were first filtered to remove the noise reads that contained additional low-quality, adaptor-polluted and high content of unknown base (N) reads by fastp software. Transcriptome assembly was performed based on the sequence of the L. vannamei genome (GenBank GCA_003789085.1) as a reference (Zhang et al., 2019). Clean data (reads) were compared with the reference genome using TopHat2 (http://tophat.cbcb.umd.edu/) to obtain the mapped data (reads) for subsequent analysis.

Functional annotation of the transcriptome was performed with 6 functional databases, including the nonredundant protein sequence database (NR), Swiss-Prot, Pfam, EggNOG (Clusters of Orthologous Groups of proteins), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), as described previously (Jinxia et al., 2015).

Using RNA-Seq by Expectation-Maximization (RSEM), the read counts of each sample gene/transcript were obtained using the results of comparison to the genome as well as the genome annotation file. They were then converted to Fragments Per Kilobase of exon model per Million mapped fragments (FPKM), which in turn yielded normalized gene/transcript expression levels. After obtaining the number of read counts of genes/transcripts, differential expression analysis was performed on multiple samples (≥2) to identify differentially expressed genes/transcripts between samples and to investigate the function of differentially expressed genes/transcripts. Differentially expressed genes (DEGs) were detected with PossionDis based on the Poisson distribution, and genes with a fold change≥2.00 and a false discovery rate (FDR) ≤0.001 were considered significant DEGs. For pathway enrichment analysis, all DEGs were annotated using GO and KEGG databases.

Eight candidate genes from the apolipoprotein and carotenoid metabolic pathways were selected for fluorescent quantitative PCR validation. Total RNA was extracted from the ovaries and hepatopancreas of shrimp at four stages of ovarian development using TRIzol™ reagent (Invitrogen, Carlsbad, CA, USA), and RNA samples were reverse transcribed using Evo M-MLV RT Premix for RT−qPCR AG11706 (ACCURATE BIOTECHNOLOGY, Hunan, China), which allows for the amplification of cDNA from total RNA to amplify cDNA without isolating mRNA. Oligonucleotide primers (Table 1) specific to eight genes and the housekeeping gene (β-actin) were designed with Primer Premier 5.0 software and synthesized by Sangon Biotech (Shanghai) Co., Ltd. Each pair of primers was examined by PCR, which was conducted using SYBR Premix Ex Taq II (Takara, Kyoto, Japan) to ensure their availability for RT−qPCR. The PCR program for every target gene and β-actin was 95°C for 30 s, followed by 40 cycles of 95°C for 10 s and 60°C for 30 s. RT-qPCR was carried out using a SYBR® Green Premix Pro Taq HS RT-qPCR Kit AG11701 (ACCURATE BIOLOGY, China), performed with a Thermal® Cycler Dice Real Time System III (TaKaRa) in final volumes of 20 μL, each containing 1 μL of cDNA, 10 μL SYBR® Green Premix Pro Taq HS RT-qPCR Kit AG11701 (ACCURATE BIOLOGY, China), 0.5 μL each of forwards and reverse primers (10 μM), and 8 μL sterile water. Each sample was analyzed in triplicate along with the internal control gene. The target gene/β-actin mRNA ratio was calculated using the formula 2˗ΔCt, and the raw data were simply transformed into the selected gene-to-actin Ct ratio for statistical analysis.

For RT-qPCR, the expression levels of target genes were routinely normalized as a ratio of β-actin mRNA detected in the same sample. For RT-qPCR and GSI measurements, the data expressed as the mean±SE (standard error) were exported using TaKaRa DiceReal Time software and analyzed by using Student’s t test or one-way ANOVA followed by Fisher’s least significant difference (LSD) test.

Unilateral eyestalk ablation (UEA) and nutritional enrichment were applied to artificially induce ovarian maturation in female L. vannamei (Figure 1A). Based on our histological results (Figure 1D), the ovaries were found to develop to stage I at 0–3 days (Ov-I, Hp-I), to stage II at 1–4 days (Ov-II, Hp-II), to stage III at 3–6 days (Ov-III, Hp-III) and to stage IV at 4–6 days (Ov-IV, Hp-IV) and to spawn at 4–6 days after these administrations (Figure 1B). In this case, the mean gonadosomatic index (GSI) of the ovaries increased from 1.15 to 6.03 from day 0 to day 6 (Figure 1C), and the GSIs for stages I, II, III and IV were 0.82, 1.95, 3.29 and 7.11, respectively (Figure 1D). The GSI increasing pattern indicated a large amount of nutrient accumulation in ovaries during the exogenous vitellogenesis period from stage III to stage IV.

After filtering out low-quality reads using Illumina HiSeq 6000, fastp and SeqPrep software, we obtained data for shrimp ovaries and hepatopancreas in four stages. For the ovarian samples, 44,267,292, 44,175,442, 44,359,992 and 44,208,868 clean reads were generated for stages I, II, III and IV, respectively, and for the hepatopancreas samples, 44,630,016, 40,197,512, 44,540,914 and 44,043,748 clean reads were generated for stage I and for stages I, II, III and IV, respectively (Table 2). All sequence reads are deposited in the National Center for Biotechnology Information (NCBI) sequence read archive (SUB10677615). The Q20 values of all samples were higher than 95% by FastQC26 analysis, indicating high sequencing quality (Reference Gene Source: Penaeus_vannamei; GenBank GCA_003789085.1).

Clean reads were pooled and assembled into nonredundant transcripts with a reference genome, and a total of 52192 transcripts were generated with an average length of more than 1800 bp (Table 3).

Comparison of raw and clean data (reads) after quality control with the reference genome using HISAT2 revealed unique mappings of 29,617,891 (66.91%), 28,914,804 (65.45%), 29,535,553 (66.58%) and 29,529,228 (66.79%) for the four ovarian stages. The unique mapping of the four stages of the hepatopancreas reached 32,624,482 (73.1%), 28,488,888 (70.87%), 30,954,251 (69.5%) and 33,394,407 (75.82%), respectively (Table 4). The NR database annotated with the most functional information on transcripts.

All genes were annotated using BLASTx with the NCBI nonredundant (Nr), KEGG, SWISS-PROT, GO, COG and PFAM protein databases. Annotation information was retrieved from proteins with the highest sequence similarity. In this study, 36,025 (90.9%) genes were annotated, and among them, 19,214, 9,695, 12,927, 14,095, 14,082 and 14,077 genes were annotated in the Nr, KEGG, SWISS-PROT, GO, COG, and PFAM databases, respectively (Figure 2A). The Venn diagram shows the annotation of expressed genes across databases, with 7,530 genes annotated to multiple libraries simultaneously, in addition to 2,540 genes annotated to the NR database. (Figure 2B).

Based on the Gene Ontology (GO) database, the genes were annotated and assigned to three categories. Within the biological process category, “cellular process” (2,502 genes) and “metabolic process” (1,855 genes) were the dominant groups. Within the cellular component category, “membrane part” (2,705 genes) and “cell part” (2,375 genes) were the most abundant groups. Within the molecular function category, “binding” (3,396 genes) and “catalytic activity” (3,143 genes) were the dominant groups (Figure 2C). Based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, a total of 15,902 genes were mapped to six specific pathways, including cellular processes, environmental information processing, genetic information processing, metabolism, human diseases, and organism systems (Figure 2D). These annotated genes were further subdivided into 43 subcategory pathways. Among them, we focused on the most annotated Metabolism section in the First Category; most genes fell under Amino acid metabolism (289 genes), Carbohydrate metabolism (414 genes), lipid metabolism (422 genes), and Glycan biosynthesis and metabolism (289 genes).

The analysis of the relationship between samples from different stages of the hepatopancreas and ovary was demonstrated by a Venn diagram (Figure 3A). In this case, a total of 6730 genes were expressed in both tissues, of which 4079 and 1227 were genes specifically expressed in the hepatopancreas and ovary, respectively. Within the hepatopancreas, 8524 transcripts were commonly expressed at four stages, while at stages I, II, III and IV, the numbers of transcripts coexpressed in the hepatopancreas and ovary were 5541, 6098, 7508 and 6917, respectively (Figure 3A).

To identify DEGs involved in ovarian development, we used the RPKM value to compare the expression differences between different groups (Figure 3B). A large number of DEGs were screened with adjusted q < 0.05 and |log2Ratio|≥ 1. A total of 10288 differentially expressed genes (DEGs) were identified between ovarian and hepatopancreatic samples, including 3077 upregulated transcripts and 7211 downregulated transcripts. The ovarian and hepatopancreatic samples were divided into four stages for further comparison. Among them, 3267 and 5138 transcripts were upregulated and downregulated at stage II, and 2737 and 6337 DEGs were upregulated and downregulated at stage III in the hepatopancreas and ovary, respectively (Figure 3C).

The three most enriched pathways of metabolism were selected from the KEGG enriched pathways amino acid metabolism, carbohydrate metabolism, and lipid metabolism, and the gene sets of the three pathways were established separately. Pathway enrichment analysis was performed using the same principle as GO functional enrichment analysis, and the KEGG pathway function was considered to be significantly enriched when the corrected P value (Padjust) was <0.05. The enrichment results (Figure 4A) showed that most of the metabolism-related transcripts were upregulated in the hepatopancreas, while fewer were upregulated in the ovary. Three pathways related to lipid transport were identified from the lipid metabolism gene set: fatty acid biosynthesis (map00061), fatty acid degradation (map00062), and biosynthesis of unsaturated fatty acids (map01040), and the heatmap showed that the majority of transcripts associated with lipid transport were upregulated in the hepatopancreas and less in the ovary (Figure 4B).

Apolipoproteins are important carriers for lipid transportation, providing the required nutrients during the shrimp ovarian development cycle. In this study, the expression levels of apolipoproteins and their receptors were selected for analysis and illustrated by a heatmap (Figure 5A). The apolipoprotein included vitellogenin (Vg, LOC11380422), apolipoprotein D (ApoD, LOC113810220 and LOC113815027), beta-1,3-guucan-binding protein (BGBP, LOC113813397 and LOC113807239), hemolymph protein (BGBP, LOC113813397 and LOC113807239) and hemolymph clottable protein (CP, LOC113819742 and LOC113807427), and the apolipoprotein receptor included vitellogenin receptor (VgR, LOC113811544 and LOC113811783), very low-density lipoprotein receptor (VLDLR, LOC113811553 and LOC113802323) and low-density lipoprotein receptor-related protein 1 (LDLR-RP1, LOC113824560) (Jia et al., 2016; Hoeger and Schenk, 2020; Li et al., 2021). In this case, Vg, ApoD, BGBP, BP, VLDLR and LDLR-RP1 transcripts were predominantly expressed in the hepatopancreas but not in the ovary, while one of the VgR transcripts appeared to have low ovarian expression and high hepatopancreas expression. The RT-qPCR results for BGBP, VgR and LDLR were consistent with the RNA-seq results (Figure 5B). However, the expression of Vg genes appeared in two different expression patterns for the two isoforms, which will be covered later in the section.

By RT-qPCR, two transcripts of Vg genes showed completely opposite expression sites (Figure 6A), in which LOC113804022 had high expression in the hepatopancreas and was defined as hepatopancreatic vitellogenin (Hp-Vg), and LOC113826723 had high expression in the ovary and was defined as ovarian vitellogenin (Ov-Vg). Compared with Ov-Vg, HP-Vg transcriptional levels are higher than ovarian fragments. The nucleotide sequence similarity between the two Vg transcripts was 71.6%, indicating that the two genes do not encode the same protein. As shown in Figure 6B, prediction of structural domains indicated the presence of DUF1943, Pfam : DUF1081 and VWD domains in both Vg proteins. Ov-Vg contained the LPD_N domain, while it was replaced by the Pfam : Vitellogenin_N domain in Hp-Vg (Figure 6B). In addition, gene organizations showed that both types of Hp-Vg and Ov-Vg genes contain 15 exons separated by 14 introns, and all introns start with GT and end with AG, conforming to the GT/AG rule for intron splicing (Figure 6B).

The primary pathway associated with carotenoid metabolism is metabolism of cofactors and vitamins, and the heatmap shows that transcripts of this pathway were more highly expressed in the hepatopancreas than in the ovary (Figure 7A). The results for the secondary pathway retinol metabolism were similar to those for the primary pathway (Figure 7B), with most transcripts highly expressed in the hepatopancreas. Four typical genes on the pathway were screened for heatmap and quantitative analysis (Figures 7C, D), and the enzyme BCOM1 (LOC113809949), associated with carotenoid cleavage, was found to be highly expressed in the hepatopancreas and almost absent from the ovary, while other representative genes cytochrome P450 3A1-like (3A1), transcript variant X1 cytochrome P450 3A13-like (3A13), and cytochrome P450 9e2-like, transcript variant X1 (9e2), which are CYP family member genes related to carotenoid processing, were also highly expressed only in the hepatopancreas, which was consistent with the RT-qPCR results.