Transcriptome analyses reveal the synergistic effects of feeding and eyestalk ablation on ovarian maturation in black tiger shrimp

To examine the effects of feeding and eyestalk ablation on ovarian maturation of domesticated female black tiger shrimp broodstock, shrimps were fed with either broodstock pellets or live polychaetes for four weeks before an eyestalk ablation was performed (Fig. 1). Both feeds could not induce ovarian maturation over the 4-week feeding period as evidenced by the pale white color of ovaries (Fig. 1A) and the insignificant changes in gonadosomatic index (GSI), a known indicator of ovarian maturation in the female black tiger shrimp (0.57 ± 0.18% to 0.54 ± 0.18% for the pellet-fed group and 0.57 ± 0.18% to 0.72 ± 0.21% for the polychaete-fed group; Fig. 1B, Supplementary Table S1).

After the eyestalk ablation, ovaries of the two groups matured differently. The ovaries of the pellet-fed group did not develop further than previtellogenic stage (stage I, see Materials and Methods for definition of ovarian stages) as they remained white-yellow in color and small in sizes with the GSIs of 0.61 ± 0.20%, 0.48 ± 0.08%, and 0.78 ± 0.31% for Day 1, 4, and 9 post ablation, respectively. On the other hand, the polychaete-fed shrimp were able to reach the reproductively mature stage after the eyestalk ablation as their ovaries continually became larger and eventually turned greenish by Day 9 post ablation. Their GSIs started from 0.72 ± 0.21% and significantly increased to 1.26 ± 0.57%, 3.61 ± 1.14%, and 5.88 ± 0.58% after eyestalk ablation at Day 1, 4, and 9, respectively (Fig. 1B, Supplementary Table S1).

To generate the reference ovarian transcriptome sequences for black tiger shrimp, we applied both short-read Illumina sequencing and long-read Pacific Biosciences (PacBio) isoform sequencing (Iso-seq). The depth coverage of Illumina sequencing allowed us to capture the sequences of rare transcripts (Supplementary Table S2) while the long PacBio reads enabled the identification of full-length transcripts (Supplementary Table S3). After processing raw read data from the PacBio RS II, we obtained a total of 546,473 reads of insert (ROIs) totaling 18.82 Gb from 8 SMRT cells (Table 1; Supplementary Table S3). Of the 546,473 ROIs, 26,609 were high-quality, full-length non-chimeric transcripts identified based on the presence of 5′ and 3′ cDNA primers and a polyA tail signal preceding the 3′-primer. The sizes of the ROIs ranged from 10 to 15,000 with an average read length of 2,087 and an N50 of 2,454 nucleotides. Short-read RNA sequencing was also carried out on the Illumina NovaSeq^TM 6000 and yielded 1,235,809,734 raw reads, totaling 187.92 Gb. After a de novo transcriptome assembly using Trinity v2.2.0¹⁷, we obtained 398,088 assembled transcripts ranging from 201 to 13,211 nucleotides. Of which, 204,819 sequences had a length of 300 nucleotides or longer (Table 1). To enhance the coverage of the reference transcriptome assembly, we combined high-quality, full-length transcripts from PacBio Iso-seq experiment with Trinity-assembled transcripts longer than 300 nucleotides. After clustering these sequences at 95% identity, we obtained a final set of 33,277 non-re dundant transcripts with an average contig length of 917 nucleotides and an N50 length of 1,728 nucleotides (Table 1). We investigated the coverage of the benchmarking universal single-copy orthologs (BUSCO)¹⁸ to assess the completeness of the gene space in the assembly. Even though we only examined transcripts deriving from ovary tissues of the black tiger shrimp, 92.7% of the assembly was identified as ‘complete’ (i.e., the length of the recovered genes are within two standard deviations of the BUSCO group mean length¹⁸), suggesting that the assembly is relatively complete with respect to the gene space. Gene ontology (GO) annotations describing biological processes, molecular functions, and cellular components were retrieved for 15,480 sequences (46.52% of the assembly; Supplementary Fig. S1). The most prevalent GO term associated with biological processes were oxidation-reduction process (239) followed by protein phosphorylation (157), translation (147), proteolysis (146), regulation of transcription (139), membrane transport (111) and signal transduction (103; Supplementary Fig. S1, Table S4). Among genes annotated to various molecular functions, the largest category was ATP binding (597), followed by metal ion binding (286) and nucleic acid binding (271), another ten subcategories with 100 or more terms (Supplementary Fig. S1, Table S4).

To confirm the expression patterns of selected transcripts with differential expression levels in the RNA-seq result, real-time qPCR was performed for nine transcripts with significant expression levels in polychaete-fed group, especially after the eyestalk ablation (Supplementary Fig. S4). All examined transcripts were differentially expressed between pellet-fed and polychaete-fed groups after the 4-week feeding trial. As observed in the RNA-seq result, vitellogenin was expressed highly in polychaete-fed group after the 4-week trial. After the eyestalk ablation, all of these transcripts were expressed higher in the polychaete-fed group than the pellet-fed group, with cytochrome c oxidase and vitellogenin being the most obvious transcripts.

Nutritional requirements for shrimp reproduction are poorly defined. Generally, live feeds such as polychaetes are used to prepare broodstock before mating, but their nutritional value can vary with species, season of collection and life stage²¹. Their use is also labor-intensive, expensive and affects water quality²². The replacement of live feeds with dry feeds for shrimp reproduction will not be possible without understanding the molecular effects of how these live feeds prepare broodstock during ovarian maturation.

From the transcriptomic data, diets showed strong influences on gene regulation. The two examined types of feed regulated totally different sets of genes throughout the experiment as evidenced by an almost complete absence of co-expressed DEGs in both feed groups, suggesting that the pellet and polychaete utilized different molecular pathways. For the pellet-fed group, genes with altered expression levels belong to regular cell maintenance, and there was no significant alteration of expression levels for genes related to reproduction. Oppositely, the polychaete-fed group differentially expressed a large number of important genes whose functions are relevant to the preparation for ovarian maturation.

Polychaetes are high in polyunsaturated fatty acids (PUFAs), which have been reported to enhance ovarian maturation in many crustacean species such as the Pacific white shrimp Litopenaeus vannamei, the Chinese mitten crab Eriocheir sinensis, the mud crab Scylla serrata and P. monodon^7–9,23. This study reveals the benefits of polychaete at the molecular level, which are highlighted by a strong induction of an oogenesis gene marker, vitellogenin, just after four weeks of polychaete feeding and prior to the eyestalk ablation (Fig. 4: Cluster IIa). Vitellogenin, also known as apolipophorin, is a precursor for the synthesis of a major egg yolk protein vitellin through vitellogenesis process²⁴. This protein provides nutrition required for crustacean embryo development²⁵. Hepatopancreas and ovary are common sites for yolk biosynthesis in several penaeid shrimps such as Litopenaeus vannamei²⁶, Penaeus semisulcatus²⁷ and Penaeus monodon²⁸. While the level of vitellogenin transcripts was significantly induced in the polychaete-fed group within the first 24 hours after the removal of the eyestalk, the same induction was not observed until Day 9 for the pellet-fed group (Fig. 4: Cluster IIa). Coordination of reproduction and yolk biosynthesis is achieved via the endocrine system, which was partly located in the eyestalk^29,30. While it is not surprising that the vitellogenin transcript was induced upon removing GIH through eyestalk ablation, it is worth noting that its expression level was induced even before the ablation in the polychaete-fed group, providing the evidence to support that ovarian maturation not only depends on the release of hormone, but also on the readiness primed by specific nutrition obtained from the polychaetes. Moreover, the fact that the vitellogenin transcript was upregulated at Day 9 after the ablation in the pellet-fed group but did not result in ovarian maturation also suggests that the completion of the maturation process requires various factors to work together.

In addition to the induction of vitellogenin, polychaete feeding differentially induced heme-binding protein and collagen alpha-1 even before the eyestalk ablation. While both of these genes have not previously been studied in P. monodon, they have been linked to reproduction in other organisms. Heme-binding protein encodes a specific yolk protein in the insect Rhodnius prolixus, and its uptake in the ovary was modulated by eicosanoids³¹. Eicosanoids, hormone-like mediators of reproduction, have been reported to be involved in female reproduction of numerous crustaceans such as Daphnia magna³². Eicosanoids or prostanoids play crucial roles in regulation of female reproductive maturation^33,34. Injection of some PUFA derivatives such as prostaglandin F(2alpha) and prostaglandin E(2) significantly increased an ovarian index and oocyte diameter in fresh water crabs³⁵. PUFAs in polychaetes are the precursors of eicosanoids; thus, it is not surprising that heme-binding protein was induced after feeding polychaetes for four weeks. Similarly, collagen alpha-1 was induced before eyestalk ablation in the polychaete-fed group. Although the roles of collagen for ovarian maturation has not been studied in the penaied shrimp, collagens are a major component of the extracellular matrix, which regulates follicle maturation and a process of explosive follicle growth called folliculogenesis in several model organisms³⁶. Various types of collagens were present throughout the surface epithelium, follicular compartments, stroma, and granulosa cells of mouse ovary³⁷.

The most interesting DEGs were those induced in the polychaete-fed group after the eyestalk ablation but not in the pellet-fed group. We hypothesized that those transcripts were up- or downregulated as a result of synergistic effects from polychaete feeding and eyestalk ablation. While most of those genes were involved in fatty acid regulation, energy production, and hormone-mediated oocyte maturation pathways, there were several genes whose functions have previously been linked to various biological processes worth examining in detail.

The shrimp culture conditions were followed a previously published work¹⁶. Female domesticated P. monodon broodstock (9-month-old) and live sand polychaetes Perinereis nuntia (6-month-old) were cultured at the Shrimp Genetic Improvement Center (SGIC, Surat Thani, Thailand). Brooders were maintained at a biosecure station containing seawater pumped from the Gulf of Thailand at a salinity of 30 ppt with water temperature at around 27 °C, and acclimated in tanks with 2 m in diameter for 2 weeks. The brooders were separated into two groups (n = 46 each): pellet-fed (Pe) and polychaete-fed (Po) groups. Pellet (Pe) and polychaete (Po) groups were fed five times/day (8:00, 11:00, 14:00, 16:30, and 22:00 h). The amount of feed/shrimp weight/day was 2% broodstock pellets (Charoen Pokphand Foods PCL) for Pe group and 7% polychaetes for Po group.

Ovaries of the brooders were collected at the beginning of the feeding trial (T0), 4 weeks after feeding trial before unilateral eyestalk ablation (T1), and 1, 4, and 9 days after unilateral eyestalk ablation (T2, T3 and T4, respectively; Fig. 1A). All ovary samples were quickly frozen in liquid nitrogen and stored at −80 °C until use. The gonadosomatic index (GSI) of each shrimp was calculated as a percentage of ovary weight to body weight, classifying different ovarian developmental stages as the following: previtellogenic (stage I, GSI < 2), vitellogenic (stage 2, GSI = 2–4), early cortical rod stage (stage 3, GSI =<4–6) and mature stage (stage 4, GSI > 6)⁵³. Differences among GSI from time points of the same feed group were considered significant Duncan; ANOVA (95% confidence).

RNA extraction was performed according a published protocol⁵⁴. Briefly, frozen ovary tissues were pulverized in liquid nitrogen for RNA extraction using TRI-REAGENT (Molecular Research Center, USA). Contaminating genomic DNA was removed by treatment with DNase I at 0.15 U/µg total RNA at 37 °C for 30 min. The quality and quantity of the RNA were assessed using gel electrophoresis and a NanoDrop (ND-8000) spectrophotometer before subsequent experiments.

To perform Iso-seq sequencing, we followed the protocol reported in Pootakham et al.⁵⁵. Equal amounts of RNA samples from all treatments were pooled (Supplementary Table S1 indicates the samples that were used for Iso-seq), and poly(A) mRNAs were enriched using the Dynabeads mRNA Direct Purification Kit (Thermo Fisher Scientific, Waltham, USA). The RNA integrity was assessed using the Fragment Analyzer (Agilent Technologies, Santa Clara, USA) prior to the construction of the Iso-seq library. We followed the Iso-seq protocol from Pacific Biosciences, using the SMARTer PCR cDNA Synthesis Kit (Clontech, Mountain View, USA) to construct the library and the BluePippin Size Selection System (Sage Science, Beverly, USA) to partition the library into 1–2 kb, 2–3 kb and 3–6 kb ranges. The sequencing was performed on a PacBio RS II instrument using P4-C6 polymerase and chemistry with 360 min movie times (Pacific Biosciences, Menlo Park, USA).

To obtain short-read RNA sequences, Illumina sequencing was performed at Macrogen (South Korea). Each library was constructed from 1 μg total RNA from the ovaries of 2–4 individuals per treatment (Supplementary Tableindicates the samples that were used for RNAseq) using the TruSeq stranded Library Preparation kit (Illumina, San Diego, USA) and sequenced using Illumina NovaSeq™ 6000 instrument with 150 bp pair end read. Sequence data were deposited in the NCBI GenBank under the BioProject ID PRJNA542029. S1

Transcriptome assembly based on the Iso-seq data was performed according to the method reported in Pootakham et al.⁵⁵. The PacBio SMRT Analysis Package (version 2.3; default parameters) was used to process raw reads from the PacBio RS II into error-corrected reads of insert. Adapter sequences, poly-A tails, artificial concatemers and 3′ truncated transcript sequences were removed using the Iso-seq protocol with default parameters. An additional round of error correction was carried out on the full-length non-chimeric transcripts using the PacBio ICE software without the Quiver step⁵⁶. If the consensus accuracy was ≥ 0.99 (default parameter; https://github.com/ben-lerch/IsoSeq-3.0/blob/master/README.md↗), the polished isoforms were classified as “high-quality.”

Illumina short reads were quality-filtered and adapter-trimmed using TrimGalore (https://github.com/FelixKrueger/TrimGalore↗). FastQC v0.11.8 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/↗) was used to check data quality before and after trimming. After the removal of low-quality reads, an Illumina-based de novo transcriptome assembly was performed using Trinity v2.2.0 using the default parameter setting¹⁷. In order to obtain a complete ovarian transcriptome, we combined high-quality, full-length transcripts from the Iso-seq experiment with Trinity-assembled transcripts. To remove redundant sequences and obtain reference transcripts from the combined set of transcripts, we clustered these transcript sequences at 95% identity using UCLUST⁵⁷. This reference transcriptome assembly was assigned function using BLASTX against custom NCBI non-redundant eukaryote protein database with an E-value cutoff of 1 × 10⁻¹⁰. BLASTX results were imported into Blast2Go Pro (version 5.2.4) for mapping and retrieving gene ontology (GO) terms associated with predicted sequences. The software also assigned the enzyme commission (EC) number and Kyoto encyclopedia of genes and genomes (KEGG) pathway annotations to the query sequences. The completeness of the transcriptome was evaluated using the coverage of the benchmarking universal single-copy orthologs (BUSCO)¹⁸.

The pair-end Illumina reads were aligned to the reference transcriptome using Bowtie 2 version 2.2.9⁵⁸. The total mapped reads in each bam file were counted by HTSeq package through the HThtseq-count union function [14].

The R package DESeq. 2 [15] provides a test for differential expression using negative binomial generalized linear models, will operate to identify the differently expressed genes between shrimp fed polychaetes or pellet at different time points after eyestalk ablation. After read count normalization, log2 fold changes between shrimp fed polychaetes or pellet after eyestalk ablation will be calculated. Significant differently expressed genes will be identified with Bonferroni adjusted p-value of <0.05. Genes displayed in Fig. 4 were selected based on their expression profiles. Our ovarian transcriptome assembly contained a total of 11,694 non-redundant features. We first selected genes that were significantly (p-value < 0.05) up- or down-regulated during at least one time point (1,559 features). We subsequently clustered them based on their expression patterns and selected groups of genes that were affected by feed effect and eyestalk ablation using Cluster 3.0 software⁵⁹. A total of 442 features exhibited changes in their expression level > 1.5 folds (Supplementary Table S5), and selected features with known functions were plotted in groups according to their similar expression patterns based on the results obtained from Cluster 3.0 software. Treeview version 1.6 was used to visualize gene expression data in Fig. 4⁵⁹. KEGG Automatic Annotation Server (KAAS) program was used to associate nucleotide sequences from 442 features to KEGG (Kyoto Encyclopaedia of Genes and Genomes) Orthology database (http://www.genome.jp/kegg/kaas/↗). Subsequently, the KEGGREST in R package was used to map KEGG Orthology (KO) to the pathways (Tenenbaum D (2019). KEGGREST: Client-side REST access to KEGG. R package version 1.26.1).

Among differentially expressed genes, nine genes were selected for qPCR validation. Primers were designed using Primer Premier Program (Table 2). We followed the previously published protocol¹⁶. Briefly, 15 μg of total RNA were treated with RQ1-RNase-free DNase (Promega, USA) to eliminate contamination of genomic DNA. The first strand cDNA were synthesized from 1.5 µg of DNase-treated total RNA using an ImProm- II^TM Reverse Transcription System Kit (Promega, USA) according to manufacturer’s instruction. Concentration and quality of the first strand cDNA were examined by Nanodrop UV (OD₂₆₀/OD₂₈₀). The cDNA of individual ovaries of black tiger shrimp from different time points were used as template. Each reaction was performed in a total volume of 20 µL containing 2x iQ™SYBR® Green Supermix (Bio-RAD, USA), 0.2 µg of cDNA template, 0.2 µM of primer pair. Cycling parameters were 95 °C for 2 min 30 s, followed by 40 cycles of 95 °C for 30 s, 58 °C for 20 s and 72 °C for 30 s. To ensure specificity of PCR products, melting curve analysis was performed from 55 °C to 95 °C with continuous fluorescence reading every 0.5 °C increment. The housekeeping gene Elongation factor 1 alpha (EF-1α) was used as an internal control⁶⁰. All qPCRs were performed in three biological replications (n = 3).

For construction of the standard curve for each transcript, a plasmid containing the transcript was constructed by cloning the PCR product of the transcript into a pGEM-T easy vector (Promega, USA). The resulting vector was transformed into E. coli JM109. The plasmid was extracted and used as the template for construction of the standard curve by 10-fold serial dilutions (10³–10⁸ copy numbers).

Expression levels between feed groups were statistically tested by t-test (p < 0.05), and those among different time points post ablation of the same feed group were statistically tested by ANOVA followed by Duncan method (p < 0.05).