In November, a total of 50 30-month-old H. schlegelii individuals were collected from Fuzhou Hongmen Reservoir of Jiangxi Province, China. The mussel gonads were immediately dissected. For each mussel, gonad tissues were sampled for RNA extraction (snap frozen in liquid nitrogen) and fixed in neutral poly-methanol stationary solution for histology. All tissues were dehydrated in an ethanol series and embedded in paraffin wax. The 5 μm sections of gonad tissue were stained with hematoxylin and eosin (H&E) according to standard protocols. Sex was identified by the presence of sperms or egg cells in the gonads. If sperms and egg cells both exist in the gonads, it is considered hermaphrodite.

The gonads of nine individuals (three males, three females, and three hermaphrodites) were used for short-read Illumina sequencing, and the gonads of two individuals (one female and one male) were used for long-read PacBio sequencing. According to the manufacturer’s instructions, total RNA was extracted using 3 ml of Trizol per sample (Invitrogen, Carlsbad, CA, United States). After treatment with RNase-free DNase I (Promega Madison, WI, United States), three male testes, three female ovaries, and three gonads from a double-sex individual were used for Illunima RNA-seq. Additionally, two samples, a male and a female, were adopted for PacBio Iso-Seq. High-throughput sequencing was performed at the Novogene Bioinformatics Institute (Novogene, Beijing, China).

The Iso-Seq method was used for sequencing long-reads of even proportions of tissue from gonads. The quality of extracted RNAs was evaluated using an Agilent 2100 Bioanalyzer (SA Pathology, Adelaide, SA, Australia). High-quality RNAs (RNA integrity number >6.0) were used for cDNA synthesis. The Iso-Seq library was prepared according to the Iso-Seq protocol using Oligo (dT) enriched mRNA. A Clontech SMARTer polymerase chain reaction (PCR) cDNA Synthesis Kit was used to reverse-transcribe mRNA into DNA. The BluePippin system was used to perform fragment screening of the full-length cDNA for damage repair, end repair, the connection of SMRT dumbbell joints, and exonuclease digestion. The Iso-Seq library was prepared as described by Pacific Biosciences (PN 100-092-800-03). Single-Molecule Real-Time (SMRT) sequencing was performed on the PacBio Sequel System.

Replicate testis and ovary RNA samples were used in mRNA paired-end library construction using a NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, United States) following the manufacturer’s recommendations. Briefly, mRNA was obtained via poly-A mRNA isolation with oligo-dT beads, and the fragments were used for cDNA synthesis with random hexamer primers (NEB). After end-repair, adenylation, adaptor ligation, cDNA purification, and PCR amplification, 3′ paired-end cDNA libraries were constructed, and their quality was evaluated on the Agilent Bioanalyzer 2100 system. After cluster generation, all libraries were loaded onto an Illumina Novaseq platform, and 150 paired-end reads were generated.

Raw Illumina reads were filtered by removal of low-quality reads and the adapters used for library construction. Reads containing undetermined nucleotides (N) and low-quality reads (Qphred ≥20) accounting for more than 50% of the entire read length were discarded. All downstream analyses were based on the high-quality clean data.

The Smrtlink (v5.1) program (http://www.pacb.com/support/software-downloads/) was used to process raw data produced by PacBio. We used the subreads.bam file in the offline data to implement the circular consensus sequence (CCS) algorithm (by ccs, parameters: ---minPasses = 2, minPredictedAccuracy = 0.8), to perform self-error correction to obtain CCS sequences (ccs.bam). These were classified by checking the polyA tail and the signals of 5 and 3′ primers to categorize the reads as FLNC (full-length non-chimera) and NFL (non-full-length non-chimera) reads by using pbtranscript classification (parameters: -min_seq_len = 200). We used the isoform-level clustering (ICE) algorithm to cluster FLNC reads and used Arrow to subsequently polish the NFL reads to create the final consensus transcriptome. The consensus Iso-Seq transcriptome was corrected with Illumina short reads using the software LoRDEC (v0.9) (Salmela and Rivals, 2014). Error-corrected consensus reads were clustered by using cd-hit-est from the CD-HIT (v4.6.8) (Fu et al., 2012) program (-c 0.95 -M 0 -T 0 -G 0 -aL 0.90 -AL 100 -aS 0.99 -AS 30) to obtain final transcripts for subsequent analysis.

To get comprehensive gene functional annotations of the unique sequences obtained after using the CD-HIT software to eliminate redundancy, we used blastn and blastx to align with NR (NCBI non-redundant protein sequences), NT (NCBI non-redundant nucleotide sequences), Swiss-Prot (a manually annotated and reviewed protein sequence database), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, setting the e-value as 1×e⁻¹⁰. Transcripts were aligned to the PFAM (Protein family) database to perform annotation using functional domains of proteins via Hmmscan. Finally, annotations were loaded into the Blast2GO (Conesa et al., 2005) program to obtain Gene Ontology (GO) terms. Transcription factors were performed using the animalTFDB 2.0 database by blastx. The Coding-Non-Coding-Index (CNCI) (Sun et al., 2013), Coding Potential Calculator (CPC) (Kong et al., 2007), Pfam-scan (Finn et al., 2016) and PLEK (Li et al., 2014) tools were used to predict the coding potential of transcripts. SSR of the transcriptome were identified using MISA (http://pgrc.ipk-gatersleben.de/misa/misa.html). It allows the identification and localization of perfect microsatellites as well as compound microsatellites which is interrupted by a certain number of bases.

Clean Illumina data were mapped to the full-length gonad transcriptomes using bowtie2 (Langmead and Salzberg, 2012), and read counts were obtained by RSEM (Li and Dewey, 2011) for each sample. For all comparisons, we considered the effects of sequence depth and gene length on fragments and normalized all read counts to the aligned RPKM (reads per kilobase per million mapped reads) to obtain the relative levels. FPKM >0.3 was defined as the threshold of significant gene expression.

The R package DESeq2 (Love et al., 2014) provides methods to identify DEGs between sex using negative binomial generalized linear models. The resulting p values were adjusted using a 5% false discovery rate to correct for multiple tests. Furthermore, a fold change |log2 (exp/ref)| > 1 was used to consider a gene as differentially expressed. Sets of DEGs and functions were obtained for the following comparisons: females versus males, hermaphrodites versus females, and hermaphrodites versus males.

To explore the potential functions of these DEGs in the three comparison groups, we performed GO and KEGG pathway enrichment analyses. DEGs were distributed into different GO terms based on their associations with cellular components, molecular functions, and biological processes. GO enrichment analysis of DEGs was executed using the clusterProfiler (Yu et al., 2012) R package (v3.14.3) based on the hypergeometric distribution. A Padj value <0.05 was assumed to indicate enriched functions. KOBAS software (3.0) (Bu et al., 2021) was used to find KEGG pathways in the statistical enrichment of DEGs.

H&E staining for visualization in histological sections was used to examine the development of the gonads. The female follicles could be observed to contain secondary oocytes and mature oocytes. Histology identified each stage of sperm cells in follicles of male gonads, including primary spermatocytes, secondary spermatocytes, spermatoblast, and sperm. Meanwhile, the morula structure during sperm development also could be observed. Both female follicle and male follicle were clearly detected in the hermaphrodite gonads (Figure 1).

The full-length transcriptome of H. schlegelii was generated using the PacBio Sequel platform, transcriptomes were obtained for both testis and ovary. PacBio single-molecule real-time sequencing is circular sequencing, and the effective inserts read produced in the single-molecule detecting during the sequencing process are called subreads. The resulting total of 12, 478, 014 subreads based on 36.1G was generated by two SMRT cells from the PacBio library. The average length of the subreads was 2,889. Standard Iso-Seq self-correction, classification, and clustering protocols (see Materials and Methods) were applied to the sequencing results. The 899,605 circular consensus sequences (CCS) were produced through conditional screening (full passes of 2 and quality of 0.80). And there were 682,624 full-length non-chimeric reads after filtering for full-length read classifications, which contained both poly-A tail and two primers. Finally, the consensus sequences for the gonads of H. schlegelii included 359,889 transcripts (Table 1). Meanwhile, more than 4.8 × 10⁸ reads numbers and a Q30 of over 90% the high-quality clean short reads were generated from nine gonads of H. schlegelii using the Illumina sequencing platform (Supplementary Table S1).

To build the full-length transcript with higher accuracy a lower error rate, all consensus sequences were corrected by short, clean reads as input data. After correcting the data, 359,889 transcripts, 3,779 N50, and 2,337 N90 were obtained (Supplementary Table S2). CD-HIT software removed redundant and similar sequences that resulted in 201,481 non-redundant transcripts with a mean length of 3,987 bp (Supplementary Table S3). About 99.2% of transcripts were longer than 1 kbp, and the main length distribution range from 2.5 to 5 kpb (Figure 2). These results show that the PacBio Iso-Seq is an efficient strategy for providing high-quality and full-length transcripts with Illumina short reads to correct.

For acquiring a comprehensive annotation of the H. schlegelii transcriptome, the 201,481 full-length transcripts were subjected to annotation analysis by aligning to public databases: NT NR, Swiss-Prot, KEGG, Pfam and GO. The Venn diagram shows that 78.31% (157,781) of the transcripts were successfully annotated in the least one database (Figure 3A). Optimal sequences alignments of sequence homology by NR database, 42.66% sequences had significant matches with Crassostrea gigas; 9.04% sequences were matched against Lottia gigantea, followed by Lingula anatina (5.18%), Aplysia california (5.09%), Octopus bimaculoide (4.58%). The remaining transcripts were homologous to those of other species (Figure 3B).

All the transcript annotated could be categorized into 55 GO terms (Figure 3C). Within the biological process category, the main GO terms were grouped in metabolic process, single-organism process and cellular component organization or biogenesis. Organelle was the most common annotation in cellular component category. Membrane and cell were the next most abundant GO terms. In the molecular function, the main proportion of GO terms referred to the binding and catalytic activity. Several putative hermaphroditism related to GO terms, including developmental process (995), reproduction (1092), and reproduction process (992), were significantly enriched.

A total of 61.65% transcripts from the gonad transcriptome of H. schlegelii were annotated with 46 subpathways from five pathways in the KEGG database (Figure 3D). The result indicated that signal transduction, endocrine system, and immune system were the top three pathways in abundance. Notably, some of these pathways may be involved in hermaphrodite, for example, oxytocin signaling pathway (1156), estrogen signaling pathway (965), and oocyte meiosis 734) (Supplementary Table S4).

A total of 87 candidate genes (1095) reportedly involved the gonad development, sex differences, and sex determination (Zhang et al., 2014; Yatsu et al., 2016; Li et al., 2020; Zhang et al., 2020; Broqurad et al., 2021) were present in the full-length gonad transcriptome of 30-month-old H. schlegelii (Supplementary Table S5). Four major sex determination-related genes in mollusks, dmrt-like, sox9, wnt4, and fem1 were present in the full-length gonad transcriptome. Eight sex-related receptor genes were present, including the four-hormone receptor esr1, gnrhr, lshr, and fshr, and other four are pgfra, spr, vgr, and zp3r. Nine spermatogenesis-related genes were significantly matched including atrx, boll, daz3, dnm3a, gtsf1, mycbp, piwl1, sycp1 and scmh1 from the full-length gonad transcriptome. Meanwhile, marf, vtg, hrom1, mta70 were associated with oogenesis, and the latter two were also related to spermatogenesis. Ten genes related to sex hormones synthesis and metabolism were identified. Four hydroxysteroid dehydrogenases hsd17b2(17β-hsd2), hsd3b1, hsd3b2, and hsd3b4. The following six genes belonged to the cyp450 family: cpy10, cyp17a1, cyp1a1, cyp1a2, cpy3a4, and cyp3a5. Eight gonad development-related genes were also identified from the full-length gonad transcriptome, including the female gonad development genes cot2, fst, vtg, vtg1, vtg6, and wnt4, and genes involved in male gonad development, including gata4 and sox9. Moreover, eight Sapta family genes and seven Spag family genes also were concerned. The remaining putatively hermaphrodite-related genes also appear in the Supplementary Table S5.

Transcription factors (TFs) play an important role in activating and inhibiting gene transcription, and they are also a current research focus. The animalTFDB2.0 database was used to identify the TFs of H. schlegelii, and the analysis predicted 61 TFs for gonads. The top 30 abundant terms are shown in (Figure 4A), of which the two most abundant were Zf-C2H2 and ZBTB. We used the unannotated transcripts to predict lncRNAs. The Venn diagram shows four methods for predicting the number of lncRNAs. Through CNCI, CPC, Pfam, and PLEK, 33,824 lncRNAs with high accuracy and reliability were identified (Figure 4B), these data aid in understanding how genes are expressed and regulated. More than 500 bp full-length transcripts were screened for SSR analysis. A total of 254,888 SSR sequences of six types (1–6 nucleotide repeats) were identified, and the distribution is shown in ( Figure 4C). Single nucleotide repeats (9-12) and dinucleotide repeats (5-8) were the most common, with counts of 76,387 and 49,402, respectively.

Differential level expression of transcripts at three shapes of H. schlegelii’s gonad was tested using the program DEseq2, and the screening criteria were Padj < 0.05 and log₂FoldChange >1. DEGs were detected between hermaphrodites and males (H vs. M), between hermaphrodites and females (H vs. F), and between males and females (M vs. F), with 2,289, 1,762, and 6,046 differentially expressed genes, respectively (Table 2). A total of 7,922 transcripts were expressed among all groups (Figure 5A). Hierarchical clustering was used to determine the distribution of DEGs among the three gonad forms (Figure 5B). Clustering showed that the male and hermaphrodite had similar expression levels, as they clustered into one branch.

For the hermaphrodite versus female group, the most enriched GO terms were “aminoglycan metabolic process” in the biological process category, “extracellular region” in the cellular component category, and “molecular function regulator” in the cellular component category. Most of the full-length transcripts were involved in the three top KEGG pathways, including “tuberculosis”, “glyoxylate and dicarboxylate metabolism” and “chemical carcinogenesis”. In the hermaphrodite versus male group, the most-enriched GO terms were “regulation of cell death”, “regulation of apoptotic process”, “regulation of programmed cell death” in the biological process category; “transferase activity, transferring alkylthio groups”, and “coenzyme-B sulfoethylthiotransferase activity” in the molecular function category. The KEGG enrichment analysis showed “tuberculosis”, “selenocompound metabolism”, and “phagosome” were the first three enriched pathways. For the male versus female group, the most-enriched GO terms were “organonitrogen compound metabolic process” in the biological process category, “extracellular region” in the cellular component category, “oxidoreductase activity” in the molecular function category. KEGG enrichment analysis concerned “glycine, serine and threonine metabolism”, “metabolism of xenobiotics by cytochrome P450”, and “porphyrin and chlorophyll metabolism” (Figure 6).

In further analysis 24 genes of the 87 early hermaphroditism-related genes were screened to have significant expression with three comparisons. Using FPKM as the expression level of 24 genes involved in hermaphroditism-related, make a gene expression heat map after normalization (Figure 7A). In total, ten genes (vtg, atrx, vtg6, wnt4, cpy17a1, nasp, cot1, vgr, cyp1a2, and marf1) were shown to be upregulated in the ovaries, while three genes (sox9, dazp1, and zp3r) exhibited higher expression in the testes. It is noteworthy of our attention that there are five genes (17β-hsd2, spg16, boll, hsp90, and spr) that are highly expressed in the hermaphrodite.

The genes 17β-hsd2, cpy1a2, and cpy17a1 were enriched in ovarian steroidogenesis, and each participates in sex hormone biosynthesis. Meanwhile, the full-length transcript is enriched to GnRH signaling pathway, oocyte meiosis, and steroid hormone biosynthesis. From these results, we may move forward in a single step to a hypothesized pathway for the formation of hermaphrodites in H. schlegelii (Figure 7B). FSH and LH are produced after being combined by GnRH and GnRHR. LH binds to LHCGR and activates the transcription of cyp17a1. FSH combinds FSHR to activate the transcription of 17β-hsd2 and cyp19a1. Under the action of cyp17a1, it catalyzes progesterone to produce androstenedione. However, there are two paths for androstenedione to generate 17β-estradiol. Androstenedione is catalyzed by 17β-hsd2 to testosterone, and then catalyzed by cyp19a1 to 17β-estradiol. Another way is that androstenedione is catalyzed by cyp19a1 to estrone, and then in 17β-estradiol was produced under the action of 17β-hsd2. Both 17β-hsd2 and cyp17a1 play important roles in the formation of sex hormones.