Skip to main content

Comprehensive analysis on the regulation of differentially expressed of mRNA and ncRNA in different ovarian stages of ark shell Scapharca broughtonii

Abstract

Background

Ovarian development is an important prerequisite and basis for animal reproduction. In many vertebrates, it is regulated by multiple genes and influenced by sex steroid hormones and environmental factors. However, relative information is limited in shellfish. To explore the biological functions and molecular mechanisms of mRNA and non-coding RNA that regulate ovarian development in Scapharca broughtonii, we performed whole transcriptome sequencing analysis on ovaries at three developmental stages. Furthermore, the biological processes involved in the differential expression of mRNA and ncRNA were analyzed.

Results

A total of 11,342 mRNAs, 6897 lncRNAs, 135 circRNAs, and 275 miRNAs were differentially expressed. By mapping the differentially expressed RNAs from the three developmental stages of Venn diagram, multiple groups of shared mRNAs and lncRNAs were found to be associated with ovarian development, with some mRNA and ncRNA functions associated with steroid hormone. In addition, we constructed and visualized the lncRNA/circRNA-miRNA-mRNA network based on ceRNA targeting relationships.

Conclusions

These findings may facilitate our further understanding the mRNA and ncRNAs roles in the regulation of shellfish reproduction.

Peer Review reports

Background

Gonadal development is an important part of the individual organism development and is the basis for reproductive activity. The stable and sustainable development of shellfish aquaculture industry depends on the supply of high quality and sufficient seedlings [1]. To maintain the steady development of seedling production, an in-depth study on the ovarian development of shellfish is necessary. Ovarian development is regulated by combination of endogenous and exogenous factors, where the endogenous factors are genetic materials such as genes, while the exogenous factors are environmental factors such as temperature, pH and photoperiod, and usually these exogenous factors regulate gonadal development by affecting the expression of the relevant genes [2]. Steroids generally play important roles as hormones in a variety of physiological processes, such as biological sex differentiation, reproduction, development and maintenance of endostasis [3]. Sex steroid hormones play various important physiological roles in vertebrates, including the regulation of growth, development, reproduction, differentiation, metabolism and homeostasis [4]. They mainly include progesterone, testosterone and estradiol, which have traditionally been defined by their role in normal reproductive function [5]. Sex steroid hormones are synthesized under the catalysis action of regulatory proteins and hormone synthase, and the genes encoding these proteins or enzymes mainly include StAR (Steroidogenic acute regulatory protein), CYP (Cytochrome P450) and HSD (Hydroxysteroid dehydrogenases) gene family members. Many earlier studies suggested sex steroid hormones existed only in vertebrates, however, more and more evidences proved they were also found in lots of mollusk [6,7,8,9], and varied seasonally with the reproductive cycle [10,11,12]. Therefore, the sex steroid hormones are also believed to be important roles in shellfish reproductive regulation.

Non-coding RNAs (ncRNAs) are a class of RNAs that do not normally encode proteins but functionally regulate protein expression [13]. Based on the size and biogenesis pathways, they are subdivided into several families including long non-coding RNAs (lncRNAs, with a length of > 200 bp), circular RNAs (circRNAs, with a closed continuous loop), and microRNAs (miRNAs, with a length of about 22 nucleotides) [14]. There are increasing evidences that lncRNAs are potential regulatory molecules as functional roles in high-order chromosomal dynamics, telomere biology and subcellular structural organization [15,16,17]. Evidences have been accumulating that circRNAs have great potential to exert effects a variety of physiological processes, such as affect microRNA regulation, regulating parental gene transcription, cell proliferation, and RNA-binding proteins [18,19,20]. Competitive endogenous RNA (ceRNA) regulation hypothesis suggests that mRNA, lncRNA, circRNA and pseudogene transcripts regulate the stability or translation activity of target genes by competitive binding with miRNA, thus achieving post-transcriptional regulation of genes [21,22,23]. Tang et al. used whole transcript sequencing to discover a dual ceRNA molecular pathway and to resolve the epigenetic mechanism of ceRNA regulation of sex determination and differentiation in Chinese tongue sole (Cynoglossus semilaevis), demonstrating for the first time in fish the regulation of sexual development by ncRNA mediated ceRNA crosstalk and revealing a new regulatory mechanism of sex determination and differentiation [24]. However, the studies of ncRNA regulation on shellfish reproduction, especially for the sex steroid hormones, still remain very limited.

The ark shell, Scapharca broughtonii, which is widely distributed along the northwest Pacific coast, is a kind of bivalve with important economic value [25]. In recent years, it has become an aquaculture variety with general approval in north coast of China because of its good market prospects, delicious taste and high nutritional value [26]. Because the gonads of S. broughtonii are covered by foot muscles, their gonad development stage cannot be judged by naked eyes, which brings trouble to the artificial breeding of seedlings. Therefore, it is very important to study the regularity and regulation of gonadal development. However, little is known about the regulating mechanisms of ncRNA in reproduction by sex steroid hormones of the ark shell. In this study, we performed whole transcriptome analysis for the ovaries from three development periods. The ceRNA networks were mapped based on the correlation between ncRNAs associated with ovarian development. These findings may facilitate our further understanding of ncRNA roles related to ovarian development in the reproductive regulation of the ark shell and provide new information for future genetic breeding in aquaculture.

Results

Determination of the period of ovary development

The gonad development in marine shellfish is distinctly cyclical and can be divided into different periods with the season change. According to previous study [27], the ovary development stage of S. broughtonii was divided into five phases based on the observation of ovary sections containing stage I (Early active stage), stage II (Development stage), stage III (Ripe stage), stage IV (Spawning stage) and stage V (Spent stage). The histological sections at five stages were obtained (Fig. 1), and the ovarian development regularity of S. broughtonii was similar to most other shellfish species. The oogonia were located on the follicle wall of ovary, and the developing oocytes were connected to the follicle wall of ovary by short stem-like structures, and the mature oocytes wander in the follicle with irregular shape (Fig. 1f). In stage I, follicles were sparse and small in number, and oogonia on follicle wall increased continuously. The number of follicles increased in stage II, and most developing oocytes were connected to the follicle wall by short stem-like structures. Stage III follicles were filled with mature free mature oocytes. The number of mature oocytes in stage IV decreased. Stage V follicles rupture and disperse, leaving only a small number of oocytes, which are gradually absorbed by surrounding cells.

Fig. 1
figure 1

The five stages histological sections of S. broughtonii ovaries. Note: a Early active stage; b Development stage; c Ripe stage; d Spawning stages; e Spent stage; f the higher magnification of oogonia, developing oocyte and mature oocyte. CT: connective tissue; FW: follicle wall; OO: oogonia; DO: developing oocyte; MO: mature oocyte; RO: residual oocyte; bar: 50 µm

Sequencing and mapping to the reference genome

After filtration, 172.84 Gb clean reads were obtained from RNA libraries of the nine S. broughtonii samples and the Q30 base percentage of each sample was more than 92.06%. Approximately 56.18% ~ 68.89% of all reads were mapped to the S. broughtonii reference genome, while 42.64% ~ 54.72% and 8.82% ~ 22.16% were uniquely and multiply mapped, respectively. Additionally, 39.55% ~ 54.80% and 39.61% ~ 54.69 of all reads were mapped to sense strand and antisense strand, respectively. A total of 103.32 M clean reads was obtained after sequencing of the small RNA libraries, with each sample clean data not less than 10.01 M. In total, 45,453 mRNAs, 69,768 lncRNAs, 3810 circRNAs and 1239 miRNAs were detected, and 21,408 novel genes were discovered in this study.

Expressional differences of mRNAs, lncRNAs, circRNAs and miRNAs

To further explore the global transcript changes, the differentially expressed mRNAs (DEGs), lncRNAs (DELs), circRNAs (DECs) and miRNAs (DEMis) were detected. Compared to group A, a total of 8302 DEGs (4344 up-regulated and 3958 down-regulated), 3788 DELs (1564 up-regulated and 2224 down-regulated), 66 DECs (24 up-regulated and 42 down-regulated) and 164 DEMis (79 up-regulated and 85 down-regulated) was identified in group B, and 8680 DEGs (4599 up-regulated and 4081 down-regulated), 4090 DELs (1866 up-regulated and 2224 down-regulated), 76 DECs (43 up-regulated and 33 down-regulated) and 178 DEMis (76 up-regulated and 102 down-regulated) were identified in group C, respectively. Compared to group B, 1918 DEGs (926 up-regulated and 992 down-regulated), 1622 DELs (903 up-regulated and 719 down-regulated), 38 DECs (14 up-regulated and 24 down-regulated) and 46 DEMis (20 up-regulated and 26 down-regulated) were screened in group C. Additionally, the heat map of the hierarchical clustering analysis clustered genes with the same or similar expression patterns. The hierarchical clustering heat maps of DEGs, DELs, DECs and DEMis expression are shown in detail in Fig. 2 and Fig. S1. The three samples from each group were consistently clustered together in the heatmap, indicating that all DEGs, DELs, DECs and DEMis were highly reproducible. These results demonstrated the high reproducibility and reliability of the transcriptome analysis in this study.

Fig. 2
figure 2

Heat map for DEGs of nine samples

Functional enrichment analysis of DEGs

In this study, a total of 11,342 DEGs were identified, and the GO and KEGG enrichment analysis based on those DEGs were performed. It was found that 1243, 1055 and 315 GO terms were significantly enriched (P-value < 0.05) in the paired comparison of group A vs B, A vs C and B vs C, respectively. A part of significantly enriched terms for the Biological Process (BP), Cellular Components (CC) and Molecular Function (MF) were displayed in Fig. 3 and Fig. S2. In the GO enrichment analysis, the most significantly enriched GO terms in the three groups was binding and cellular process.

Fig. 3
figure 3

The GO and KEGG enrichment analysis of the DEGs between group A vs B. a GO enrichment analysis of DEGs in group A vs B. b KEGG enrichment analysis of DEGs in group A vs B

KEGG enrichment analysis was performed on all DEGs to further determine the metabolic process pathways. The results revealed that the DEGs were significantly enriched in 53 KEGG pathways. As shown in Fig. 3 and Fig. S2, the top 20 significant KEGG pathways for each group were chosen and displayed. Of those, the most significantly enriched pathway for DEGs in groups A vs B, A vs C and B vs C were “FoxO signaling pathway”, “Ribosome” and “Oxidative phosphorylation”, respectively. Twenty of the KEGG pathways significantly enriched in the A vs B and A vs C groups were identical, including “Ribosome”, “FoxO signaling pathway”, “Hippo signaling pathway” and “DNA replication”. The “Ribosomal” pathway facilitates oocyte development and may be associated with ovarian development. The primary biochemical and signal transduction pathways identified by the KEGG analysis will provide further insight into the future research direction of mRNA.

To investigate the potential functions of DEGs related to ovarian development, all DEGs from three groups were plotted on a Venn diagram. As illustrated in Fig. 4a, the DEGs shared by the three groups had 404 genes. Of these, some genes were annotated with functions related to gonad development and sex steroid hormones, such as CYP17A, StAR-related lipid transfer protein 10 and steroid dehydrogenase. The GO enrichment analysis of those DEGs was shown in Fig. 4b, and the top 10 significantly over-represented GO terms were related to “binding”, “cellular process” and “metabolic process”. As shown in Fig. 4c, the top 20 significant KEGG terms included “Ovarian steroidogenesis”, “GnRH secretion” and “Prolactin signaling pathway” which were involved in ovarian development and sex steroid hormone synthesis.

Fig. 4
figure 4

Venn diagram of DEGs among three groups and GO and KEGG enrichment analysis plots of the 404 DEGs shared among groups. a Venn diagram of all DEGs among three groups; b GO enrichment analysis of the 404 DEGs shared among three groups; c KEGG enrichment analysis of the 404 DEGs shared among three groups

Functional analysis of DELs and their target genes

lncRNAs with coding potential were filtered out using the software CPC2, CNCI, CPAT and Pfamscan, and a total of 69,768 lncRNAs were identified. Most of these were defined as intergenic lncRNA, accounting for 75.1%, and followed by intron-lncRNA accounting for 20.9%. Further, using the screening criteria of Fold Change > 1.5 and P-value < 0.05, 6897 DELs was obtained. To investigate the possible functions of lncRNA associated with ovarian development, a Venn plot of the three groups was drawn as Fig. 5a. It showed that the DELs shared among the three groups had 70 lncRNAs related to 21,682 target genes. The GO enrichment analysis showed that the top 10 significantly enriched GO terms were “binding”, “cellular process” and “metabolic process” (Fig. 5b). Seven KEGG pathways were significantly enriched (P < 0.05), and the top 20 significant KEGG terms included “Steroid biosynthesis” and “Ovarian steroidogenesis” (Fig. 5c). Of these, some target genes were associat with sex steroid hormones, such as 3 beta-hydroxysteroid dehydrogenase, 17β-HSD14, StAR -related lipid transfer protein 5, steroid 17-alpha-hydroxylase/17, 20 lyase-like isoform X2, Estrogen receptor and Steroid receptor RNA activator 1.

Fig. 5
figure 5

Venn diagram of DELs shared among three groups and GO and KEGG enrichment analysis plots of the target genes. a Venn diagram of all DEGs among three groups; b GO enrichment analysis on the target genes of 70 DELs shared among three groups; c KEGG enrichment analysis on the target genes of 70 DELs shared among three groups

Functional analysis on DECs

The length distribution of circRNAs were analyzed (Fig. 6a), and it was found that the sequence length of most circRNAs were greater than 3000 bp, and followed by the length range from 400 bp to 600 bp. There were 3810 circRNAs were identified using the CIRI software, of which 70%, 20% and 10% belonged to exon, intergenic region and intronic type, respectively (Fig. 6b). And a total of 135 DECs were detected among the three groups by using Fold Change > 1.5 and P-value < 0.05 as the screening criteria. We identified 66 DECs in group A vs B, 76 in groups A vs C and 38 in group B vs C. All host genes of DECs were functionally annotated, and the GO enrichment were shown in Fig. 6 and Fig. S3. The binding and cellular processes were the most significant among the three groups. KEGG pathway enrichment analysis was performed on DECs host genes, and 12 KEGG pathways were enriched including “Insulin signaling pathway”, “FoxO signaling pathway” and “Lysosome”.

Fig. 6
figure 6

Distribution of circRNAs and GO and KEGG enrichment analysis on DECs host gene in group A vs B. a Distribution diagram of circRNA length; b circRNA distribution in different regions of the genome; c GO enrichment analysis on the DECs host gene in group A vs B; d KEGG enrichment analysis on the DECs host gene in group A vs B

Construction of potential lncRNA/circRNA-miRNA-mRNA regulatory networks

To identify ceRNA interactions associated with ovarian development, we constructed ceRNA networks based on the interaction of differentially expressed mRNA and ncRNA, and the ceRNA networks were constructed and visualized using the cytoscape software. Among them, part of the network related to ovarian development was shown in Figs. 7 and 8. As shown in Fig. 7, the lncRNA-miRNA-mRNA network related to ovarian development consisted of 6 miRNAs, 36 lncRNAs and 60 genes, among which 4 genes were involved in ovarian related pathways, such as "Oocyte meiosis" and "GnRH signaling pathway". The circRNA-miRNA-mRNA network (Fig. 8) included 10 circRNAs, 4 miRNAs and 87 DEGs, including 6 ovarian development related genes such as CYP17, 17β-HSD12 and β-catenin. CYP17 encodes 17α-hydroxylase/17, 20-lyase, a key enzyme in the steroid synthesis pathway with significant roles in steroid synthesis, sex determination and differentiation and gametogenesis. In the ceRNA network, two CYP17 genes (EVM0010244, EVM0015185) both interacted with a miRNA (novel_miR_1238), and five circRNAs (Circ00168055, Circ18481716, Circ28051217, Circ08882033, Circ08909412) acted as sponges in this process. CYP17 is key rate-limiting enzyme in the sex steroid hormone synthesis pathway and plays an important role in steroid hormone synthesis of vertebrates. The novel_miR_1238 interacts with two CYP17 genes, suggesting that novel_miR_1238 may indirectly regulate CYP17 genes and thus participate in the regulation of sex steroid hormone synthesis. These results suggest that lncRNA and circRNA may function as ceRNAs, competing with miRNAs as a binding partner of mRNA, which may regulate the metabolism involved in the synthesis of sex steroid hormones.

Fig. 7
figure 7

Interaction diagram of lncRNA-miRNA-mRNA network associated with ovarian development. Red and green nodes represent up-regulated and down-regulated RNAs, respectively

Fig. 8
figure 8

Interaction diagram of circRNA-miRNA-mRNA network associated with ovarian development. Red and green nodes represent up-regulated and down-regulated RNAs, respectively

Quantitative real-time PCR validation

To validate the accuracy and reliability of transcriptome sequence data, a total of 24 differentially expressed RNAs were selected for qRT-PCR analysis, including three mRNAs, three lncRNAs, one circRNA and one miRNA in each groups. As shown in Fig. 9, the qRT-PCR results were coincided with the RNA sequence data, which confirmed the reliability of RNA sequence analysis.

Fig. 9
figure 9

Validation of the RNA-seq data by qRT-PCR. a RNA-seq and qRT-PCR data for group A vs B; b RNA-seq and qRT-PCR data for group A vs C; c RNA-seq and qRT-PCR data for group B vs C

Expression profile of nine genes related to sex steroid hormone

The expression pattern of nine genes associated with sex steroid hormones and reproduction in different stages of ovarian development were analyzed. During ovarian development, nine related genes exhibited two different expression patterns. The relative expression levels of 17BHSD14, CYP17A, CYP3A28, Steroid receptor, START5 and Beta-catenin showed an increasing trend during ovarian development from stage I to stage II, and reached the highest value in stage II, then there was downward trend during ovarian development from stage II to stage V (Fig. 10 a-f). The relative expression levels of 3BHSD7, CYP17 and Catenin beta-like increased during ovarian development from stage I to stage II, and reached the highest value in stage II, then decreased during ovarian development from stage II to stage IV, and slightly increased during ovarian development from stage IV to stage V (Fig. 10 g-i).

Fig. 10
figure 10

The expression pattern of nine related genes. a 17BHSD14, b CYP17A, c CYP3A28, d Steroid receptor, e START5, f Beta-catenin, g 3BHSD7, h CYP17, i Catenin beta-like. The different lowercase letters represent significant differences in expression of genes (P < 0.05)

Discussion

In this study, we explored the different expression patterns of mRNAs and ncRNA for the first time during the different development periods of ovary for screening the key factors influencing the breeding of ark shell at the whole-transcriptome level. In addition, the potential lncRNA/circRNA-miRNA-mRNA regulatory networks were constructed, which provided new insights into the molecular mechanisms underlying the effects of sex steroid hormones on shellfish ovarian development.

DEGs with known functions associated with ovary development and sex steroid hormones synthesis were discovered through comparative transcriptome analysis of overies at different development stages, such as CYP17A, StAR-related lipid transfer protein 10 and steroid dehydrogenase. CYP17A is a key rate-limiting enzyme in the biosynthesis of sex steroid hormones, progesterone and glucocorticoids [28]. Its activity and expression are important in regulating gonad development, germ cell differentiation and maintenance of secondary sexual characteristics [29]. The expression of CYP17A was high in the early stage of gonadal development and has sex dimorphism of Hyriopsis cumingii. After interfering with StAR3 gene, the expression level of CYP17A gene in gonad decreases, as did the content of sex steroid hormone, indicating that CYP17A gene plays important roles in gonadal development and participates in sex determination and gonadal development through its role in the synthesis pathway of sex steroid hormone [30]. This indicated that these differential expression CYP17 genes play important roles in the process of ovarian development and sex steroid hormones synthesis. The StAR transport of cholesterol from the outer mitochondrial membrane to the inner mitochondrial membrane is the first and rate-limiting step in sex steroid hormone synthesis [31]. Namely, StAR is an important gene in the rate-limiting step of the sex steroid hormone synthesis pathway and regulates the rate and level of sex steroid hormone synthesis. In Chlamys farreri, StAR3 gene was highly expressed in the gonadal tissue at the early stage of differentiation, and then decreased again in the ovaries at maturity after its expression decreased during growth, which may be because StAR3 provides an important enzyme for steroid hormone synthesis in the early stage of ovarian cell differentiation and ovarian maturation [32]. The StAR-related lipid transfer protein 10 may play an significant role in the transport of cholesterol during the early stages of ovarian development. The expression levels of CsDMRT3 and CsFOXL2 genes related to reproduction reached their peak in the mature stage of testes and ovaries development in the clam Cyclina sinensis, respectively. CsDMRT3 and CsFOXL2 expression levels were highest in the gonads when estradiol was treated with 5 µg/L, and oocyte diameter became larger with the increase of estradiol concentration. The results showed that, CsDMRT3 and CsFOXL2 play important roles in gonadal development and sex differentiation of C. sinensis [33, 34].

Non-coding RNAs, as important members of the intracellular gene regulatory network, influence various cellular life activities. Current studies suggest that non-coding RNAs also play important roles in the reproductive activities and participates in the regulation of spermatogenesis and oogenesis [35]. In recent years, researchers have discovered and identified a large number of non-coding RNA, but the studies on non-coding RNAs are very limited in aquatic animals. In this study, 69,768 lncRNAs, 3810 circRNAs and 1239 miRNAs were identified in nine samples, and 6897 DELs, 135 DECs and 164 DEMis were also detected. In mammals, ncRNAs have been shown to play a crucial role in the development of gonads [36, 37]. High-throughput sequencing of porcine ovaries revealed that circRNAs may play potential roles in regulating porcine litter size [38]. Currently, the main focus in aquatic animals is on the post-transcriptional regulation of non-coding RNA. Through miRNA sequencing analysis of Chinese mitten crab (Eriocheir sinensis) ovary, Song et al. found that miR-2 and miR-133 were differentially expressed during meiotic maturation in oocytes, indicating that miRNA may regulate cyclin B expression at the post-transcriptional level [39]. Transcriptome sequencing of the testes and ovaries of large yellow croaker (Larimichthys crocea) revealed multiple lncRNA targeting known genes related to sex differentiation and developmental. Among them, lncRNA MSTRG.24346 is extremely significant correlated with a sex determination candidate gene dmrt1, indicating that lncRNA may play an important role in sex differentiation [40]. In this study, we analyzed the expression changes of ncRNAs during ovarian development in ark shell, and screened the differentially expressed ncRNAs during different ovarian development stages by comparative analysis, and further performed target gene prediction and biological function annotation of differentially expressed ncRNAs with the help of bioinformatics tools. Functional annotation revealed that some lncRNAs were involved in reproductive signaling pathways such as steroid biosynthesis, ovarian steroidogenesis and estrogen signaling pathway involving target genes such as 3β-HSD, 17β-HSD14, StAR-associated lipid transfer protein 5, estrogen receptor and steroid receptor RNA activator 1, which are related to reproduction.

According to the ceRNA hypothesis, lncRNAs and circRNAs can function as miRNA's molecular sponges to restore the expression of mRNAs that miRNA has suppressed. The primary pathway of ncRNA-mRNA interactions is through ceRNA networks. Functional research of miRNA has found that ceRNA (most commonly lncRNA and circRNA) can regulate gene expression by competitively binding to miRNA. The ceRNA can affect miRNA-induced gene silencing through the miRNA via miRNA response elements (MREs) [41]. Sequencing of the longissimus dorsi muscle from two different pigs identified potential lncRNAs/circRNAs-miRNAs-mRNAs regulatory networks sharing genes involved in skeletal muscle muscular proliferation, differentiation/regeneration, and adipogenesis. These differentially expressed non-coding RNA as ceRNA for miRNAs may be involved in the molecular basis of muscle traits [42]. By mapping the known ovarian development related genes in ark shell, the ceRNA of the differentially expressed non-coding RNA and target genes were also analyzed. The lncRNA-miRNA-mRNA network consisted of 6 miRNAs, 36 lncRNAs and 60 genes, among which 4 genes were involved in ovarian related pathways. The circRNA-miRNA-mRNA network showed that a total of 101 RNAs interacted with each other, including 10 circRNAs, 4 miRNAs and 87 DEGs. And this circRNAs-miRNAs-mRNAs regulatory network includes six genes related to ovarian development. A total of 10 miRNAs associated with ovarian development were identified in this study, all were novel miRNAs. Furthermore, 36 lncRNAs and 10 circRNAs participated in the lncRNA/circRNA-miRNA-mRNA networks related to ovarian development, respectively, which may play potential roles in ovarian development. The lncRNA/circRNA-miRNA-mRNA network analysis revealed that miRNAs target and regulate genes were associated with ovarian development and sex steroid hormones. We found that the novel_miR_1238 regulates two CYP17 genes, which encode 17α-hydroxylase/17, 20-lyase that are key enzymes in the sex steroid hormone synthesis pathway. CYP17 gene controls the direction of steroid synthesis and is one of the key rate-limiting enzymes in sex steroid hormones biosynthesis pathway [43]. Therefore, novel_miR_1238 may indirectly regulate CYP17 gene to participate in sex steroid hormone synthesis.

Conclusions

We compared the expressional characteristics of mRNAs and ncRNAs in the ovarian tissues of S. broughtonii. According to the results, ncRNAs were abundant in the ovarian tissues, and bioinformatics analyses indicated that ncRNAs were involved in ovarian development and sex steroid hormone. We identified 10 miRNAs, 36 lncRNAs and 10 circRNAs that may be involved in the regulation of ovarian development. This study provides new information about the effects of genetic and ncRNA on shellfish ovarian development and sex steroid hormone synthesis and provides new insights into future genetic breeding for aquaculture.

Methods

Sample

A total of 200 S. broughtonii individuals with an average shell length of 80.74 ± 5.08 mm at different stages of gonadal development were collected from Yantai in Shandong Province. Before sampling, all ark shells were cultured in tanks containing fully aerated filtered seawater at salinity of 30 ± 1, temperature of 18 ± 1 °C, pH of 7.75 ± 0.25 and dissolved oxygen of 7.56 ± 0.23 mg/L for one week. During the culture period, 100% of the seawater was replaced daily with fresh seawater, and fresh single-celled algae was provided every four hours. All sampled gonad tissues were orderly dissected, placed in liquid nitrogen, and then stored in a freezer at –80 ℃. Simultaneous, part of dissected gonads were placed in Bouin's fixative to fix for 24 h and then stored in 70% of the ethanol, used in paraffin sections to determine the gonad development period. Gonadal tissue sections were visualized using a Leica microscope to determine the specific periods of gonad development. In the following studies, nine ovarian samples were randomly selected from stages I (group A), stages III (group B), and stages IV (group C) for subsequent sequencing, with three replicates per group.

RNA extraction and quality monitoring

Nine ovarian samples from three developmental stages, with three replicates, were selected for RNA extraction. Total RNAs were extracted by using the miRNeasy Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. RNA concentration was measured using Nanodrop 2000 (Thermo Fisher Scientific, Wilmington, DE, USA), and RNA integrity was analyzed with agarose gel electrophoresis and the Agient 2100 Bioanalyzer (Agilent Technologies, Foster City, CA, USA). Samples with an RNA integrity number (RIN) value greater than seven were used for further analysis. The nine RNA samples that met the requirements were used for further constructing transcriptome library and sequencing.

RNA library construction and sequencing

Both the RNA libraries and the small RNA libraries were constructed using 3 µg RNA per sample. After removing ribosomal RNA (rRNA) using the Ribo-ZeroTM kit (Epicentre, Madison, WI, USA), the RNA libraries were built using the VAHTS Universal V6 RNA-seq Library Prep Kit for Illumina® (Vazyme, China) for mRNA, lncRNA and circRNA sequencing. The small RNA libraries were constructed using the VAHTSTM Small RNA Library Prep Kit for Illumina (Vazyme, China) after rRNA was removed. The concentration and insert size were detected using Qubit 3.0 and Agilent LabChip GX 2100, respectively, following the construction of the libraries. The libraries were sequenced on the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA).

Reads mapping and transcriptome assembly

In order to ensure the accuracy of information analysis, quality control of raw data was required to obtain high-quality reads order, which were clean reads. The clean data was obtained by removing low-quality reads containing over 10% of poly (N) and adapters from the raw data. And then the FastQC (http://www.bioinformatics. babraham.ac.uk/projects/fastqc/) was used to verify the sequence quality. The clean data were mapped to the S. broughtonii reference genome [44] sequence by the HISAT2 (version 2.0.4). Mapped reads were assembled into transcripts using StringTie (version 1.3.1). The assembled transcripts were annotated using the gffcompare program (http://ccb.Jhu.edu/software/stringtie/gffcompare.shtml).

Identification of lncRNA, circRNA, and miRNA

The identification of lncRNA were carried out by two steps including basic screening and potential coding capability screening. The transcripts were selected firstly by exon number ≥ 2, length ≥ 200 bp and the FPKM ≥ 0.1. And then, the coding potential calculator 2 (CPC2) [45], coding-non-coding index (CNCI, v2) [46], Coding Potential Assessment Tool (CPAT) [47] and Pfamscan (version 1.3) [48] were used to filter out the transcripts with coding potential. Transcripts that pass through all those above procedures were considered lncRNAs. CircRNAs were identified with the CIRI software [49]. The CIRI used the BWA software [50] to mapping with the reference gene sequence to generate SAM files. The CIGAR values in the SAM file were then analyzed, and the PCC (Paired Chiastic Clipping) signals was scanned from the SAM file. We mapped the reads to the reference genome to the mature sequences of the known miRNA in the miRBase (v22) database and their upstream and downstream for known miRNAs identification. For the biometric features of miRNAs, the new miRNAs were predicted using miRDeep2 software (version 2.0.5) [51]. The possible precursor sequences were obtained by reads alignment to the location information on the genome. Based on the distribution information of reads on the precursor sequence and the precursor structure energy information, the Bayesian model was scored to finally realize the prediction of the new miRNAs.

Target gene prediction and annotation

Neighjacent genes within 100 kb upstream and downstream of lncRNA were used as cis target genes using a Perl script. Trans target genes of lncRNA were predicted with correlation coefficient value > 0.9 and significance P-value < 0.01. Target gene of miRNA was predicted by using miRanda [52] and targetscan [53] based on the gene sequence information of the known miRNA and the newly predicted miRNA with the corresponding species. Since circRNAs comprise several miRNA binding sites, it is possible to identify circRNAs that bind to miRNAs using miRNA target gene prediction, which may then be used to determine the function of this subset of circRNAs based on the functional annotation of miRNA target genes. The predicted target gene sequences were aligned to the databases of NR, Swiss-Prot, GO, COG, KEGG, KOG and Pfam to obtain the annotation information of the target genes.

Differentially expressed RNA, GO and KEGG enrichment analyses

The DESeq software package (http://bioconductor.org/packages/release/bioc/html/ DESeq.html) was used to detect differentially expressed genes, differential expressed lncRNAs, differential expressed circRNAs and differential expressed miRNAs by satisfying different value of log2FoldChange|and FDR. The ClusterProfiler software was used to perform Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) [54] pathway enrichment analysis on DEGs, DELs, DECs and DEMIs, respectively. The KEGG enrichment analysis of differential genes used hyper-geometric distribution testing with ClusterProfiler software. GO terms and KEGG pathways with Benjamin-corrected P-value < 0.05 were considered to be significantly enriched.

Construction of the lncRNA/circRNA-miRNA-mRNA network

It has been proved that the lncRNAs/circRNAs could regulate the mRNA expression by competitively binding to the MREs. To explore the possible interactions between miRNA and mRNA, lncRNA and circRNA, respectively, the miRanda software was used to analyze the combinative ability between the S. broughtonii miRNA and the DEMs, DELs and DECs. Based on the ceRNA targeting relationship, the lncRNA/circRNA-miRNA-mRNA network was constructed and visualized using the Cytoscape software V3.9.0 (San Diego, CA, USA). Since lncRNA, circRNA, miRNA and mRNA can form endogenous competitive mechanisms, special attention was paid to target genes negatively correlated with miRNA expression levels in the analysis. That only follow up—down—up "or" down—up—down "trend of the expression of selected for further research.

Quantitative real-time PCR validation

The accuracy of omics data was validated by the relative expression levels of selected mRNAs, lncRNAs, circRNAs and miRNAs using quantitative real-time PCR (qRT-PCR). In the qRT-PCR experiment, the β-actin [55] was set as internal reference genes for mRNA, lncRNA, and circRNA and mir-125-x [56] for miRNA. The RNA used in the qRT-PCR was the same as that used in sequencing. The cDNA for mRNA, lncRNA and circRNA amplification was converted by HiScript® III RT SuperMix for qPCR Kit, and then the qRT-PCR reactions were performed using the ChamQ SYBR Color qPCR Master Mix Kit according the manufacturer’s instructions. Meanwhile, for miRNA, miRNA 1st Strand cDNA Synthesis Kit and miRNA Universal SYBR qPCR Master Mix were respectively used for cDNA synthesis and qRT-PCR. Each qRT-PCR experiment was carried out in triplicate. Relative gene expression levels were calculated using the 2−∆∆Ct method. The primers designed by Primer 3.0 were listed in Table 1. All data were analyzed using one-way ANOVA in SPSS 22.0. Quantitative data were presented as means ± standard deviation (SD), and the significance level for the statistical analysis was set as P < 0.05.

Table 1 Primers used in qRT-PCR experiment

Expression changes of nine genes during ovarian development

Three ovarian samples were randomly selected from each orary development stage, and from the identified DEGs, nine genes associated with sex steroid hormones and reproduction were selected for the qRT-PCR experiment, including CYP17A, CYP17, CYP3A28, 3BHSD7, 17BHSD14, START5, Steroid receptor, Beta-catenin and Catenin beta-like, then their expressionwere analyzed as aboved mentioned. The primers designed by Primer 3.0 were listed in Table 2.

Table 2 Primers used in qRT-PCR experiment of related genes

Availability of data and materials

All sequencing data was deposited in the NCBI Short Read Archive (SRA) database under the BioProject ID: PRJNA948553. Relevant supporting data can be found within the article and additional files.

Abbreviations

ncRNAs:

non-coding RNAs

lncRNAs:

long non-coding RNAs

circRNAs:

circular RNAs

miRNAs:

microRNAs

ceRNA:

competitive endogenous RNA

BP:

Biological Process

CC:

Cellular Components

MF:

Molecular Function

StAR:

Steroidogenic Acute Regulatory Protein

MREs:

miRNA response elements

RIN:

RNA integrity number

rRNA:

ribosomal RNA

PCC:

Paired Chiastic Clipping

DEGs:

differentially expressed genes

DELs:

differentially expressed lncRNAs

DECs:

differentially expressed circRNAs

DEMis:

differentially expressed miRNAs

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

qRT-PCR:

quantitative real-time PCR

SD:

standard deviation

References

  1. Qing W. Preliminary study on gonadal development and mechanism of oocyte maturation of hard clam Mercenaria mercenaria. Qingdao: Institute of Oceanology Chinese Academy of Sciences; 2007.

  2. Yamaguchi A, Lee KH, Fujimoto H, Kadomura K, Yasumoto S, Matsuyama M. Expression of the DMRT gene and its roles in early gonadal development of the Japanese pufferfish Takifugu rubripes. Comp Biochem Physiol D Genomics Proteomics. 2006;1(1):59–68.

    Article  PubMed  Google Scholar 

  3. Henley DV, Lindzey J, Korach K. Steroid hormones. In: Shlomo Melmed, P Michael Conn. Endocrinology. Humana Press; 2005.

  4. Edwards DP. Regulation of signal transduction pathways by estrogen and progesterone. Annu Rev Physiol. 2005;67(1):335–376. https://doi.org/10.1146/annurev.physiol.67.040403.120151.

    Article  CAS  PubMed  Google Scholar 

  5. Wierman ME. Sex steroid effects at target tissues: mechanisms of action. Adv Physiol Educ. 2007;31(1):26–33. https://doi.org/10.1152/advan.00086.2006.

    Article  PubMed  Google Scholar 

  6. Siah A, Pellerin J, Benosman A, Gagné JP, Amiard JC. Seasonal gonad progesterone pattern in the soft-shell clam Mya arenaria. Comp Biochem Physiol A Mol Integr Physiol. 2002;132(2):499–511.

    Article  CAS  PubMed  Google Scholar 

  7. Reis-Henriques MA, Guellec DL, Remy-Martin JP, Adessi GL. Studies of endogenous steroids from the marine mollusc Mytilus edulis L. by gas chromatography and mass spectrometry. Comp Biochem Physiol. 1990;95(2):303–9.

    Google Scholar 

  8. Matsumoto T, Osada M, Osawa Y, Mori K. Gonadal estrogen profile and immunohistochemical localization of steroidogenic enzymes in the oyster and scallop during sexual maturation. Comp Biochem Physiol Part B. 1997;118(4):811–7.

    Article  Google Scholar 

  9. Le Guellec D, Thiard MC, Remy-Martin JP, Deray A, Gomot L, Adessi GL. In vitro metabolism of androstenedione and identification of endogenous steroids in Helix aspersa. Gen Comp Endocrinol. 1987;66(3):425–33.

    Article  PubMed  Google Scholar 

  10. Ketata I, Guermazi F, Rebai T, Hamza-Chaffai A. Variation of steroid concentrations during the reproductive cycle of the clam Ruditapes decussatus: a one year study in the gulf of Gabès area. Comp Biochem Physiol A Mol Integr Physiol. 2007;147(2):424–31. https://doi.org/10.1016/j.cbpa.2007.01.017.

    Article  CAS  PubMed  Google Scholar 

  11. Liu J, Zhang Z, Zhang L, Liu X, Yang D, Ma X. Variations of estradiol-17β and testosterone levels correlated with gametogenesis in the gonad of Zhikong scallop (Chlamys farreri) during annual reproductive cycle. Can J Zool. 2014;92(3):195–204. https://doi.org/10.1139/cjz-2013-0202.

    Article  CAS  Google Scholar 

  12. Yan H, Li Q, Liu W, Ke Q, Yu R, Kong L. Seasonal changes of oestradiol-17β and testosterone concentrations in the gonad of the razor clam Sinonovacula constricta (Lamarck, 1818). J Molluscan Stud. 2011;77(2):116–22.

    Article  Google Scholar 

  13. Mattick JS, Makunin IV. Non-coding RNA. Hum Mol Genet. 2006;15 Spec No 1(suppl 1):R17.

    Article  PubMed  Google Scholar 

  14. Cech TR, Steitz JA. The noncoding RNA revolution—trashing old rules to forge new ones. Cell. 2014;157(1):77–94.

    Article  CAS  PubMed  Google Scholar 

  15. Mercer TR, Dinger ME, Mattick JS. Long non-coding RNAs: insights into functions. Nat Rev Genet. 2009;10(3):155–9. https://doi.org/10.1038/nrg2521.

    Article  CAS  PubMed  Google Scholar 

  16. Ulitsky I, Shkumatava A, Jan CH, Sive H, Bartel DP. Conserved function of lincRNAs in vertebrate embryonic development despite rapid sequence evolution. Cell. 2011;147(7):1537–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Batista PJ, Chang HY. Long noncoding RNAs: cellular address codes in development and disease. Cell. 2013;152(6):1298–307.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Chen I, Chen CY, Chuang TJ. Biogenesis, identification, and function of exonic circular RNAs. Wiley Interdiscip Rev RNA. 2015;6(5):563–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Lasda E, Parker R. Circular RNAs: diversity of form and function. RNA. 2014;20(12):1829–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Xuan L, Qu L, Zhou H, Peng W, Sun Y. Circular RNA: a novel biomarker for progressive laryngeal cancer. Am J Transl Res. 2016;8(2):932–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  21. Denzler R, McGeary SE, Title AC, Agarwal V, Bartel DP, Stoffel M. Impact of MicroRNA levels, target-site complementarity, and cooperativity on competing endogenous RNA-regulated gene expression. Mol Cell. 2016;64(3):565–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Xia T, Liao QI, Jiang X, Shao Y, Xiao B. Long noncoding RNA associated-competing endogenous RNAs in gastric cancer. Sci Rep. 2014;4:6088.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Peng L, Chen G, Zhu Z, Shen Z, Du C, Zang R, Su Y, Xie H, Li H, Xu X. Circular RNA ZNF609 functions as a competitive endogenous RNA to regulate AKT3 expression by sponging mir-150-5p in Hirschsprung’s disease. Oncotarget. 2017;8(1):808–18.

    Article  PubMed  Google Scholar 

  24. Tang L, Huang F, You W, Poetsch A, Nóbrega RH, Power DM, Zhu T, Liu K, Wang H, Wang Q, et al. ceRNA crosstalk mediated by ncRNAs is a novel regulatory mechanism in fish sex determination and differentiation. Genome Res. 2022;32(8):1502–15.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Zhao Q, Wu B, Liu Z, Sun X, Zhou L, Yang A, Zhang G. Molecular cloning, expression and biochemical characterization of hemoglobin gene from ark shell Scapharca broughtonii. Fish Shellfish Immun. 2018;78:60–8.

    Article  CAS  Google Scholar 

  26. Tanaka T, Aranishi F. Comparative genetic characterization of ark shell Scapharca broughtonii in Northeast Asia. J Shellfish Res. 2016;35(2):421–7. https://doi.org/10.2983/035.035.0216.

    Article  Google Scholar 

  27. Fei TR. Studies on reproductive biology and artificial breeding technology of Scapharca broughtonii. Qingdao: Ocean University of China; 2014.

    Google Scholar 

  28. Dongjie S, Zheng T, Qiang Z, Shihong W, Dong W, Wentong L, Saisai W, Yansheng S, Jufeng J. Expression of cytochrome CYP17a1 in different developmental phases of gonad in koi carp(Cyprinus carpio). J Henan Agric Sci. 2022;51(2):139–46.

  29. Su T, Ijiri S, Kanbara H, Hagihara S, Wang DS, Adachi S. Characterization and expression of cDNAs encoding P450c17-II (cyp17a2) in Japanese eel during induced ovarian development. Gen Comp Endocrinol. 2015;221:134–43.

  30. Shangguan XZ, Mao YR, Wang XQ, Liu ML, Wang YY, Wang GL, Li JL. Cyp17a effected by endocrine disruptors and its function in gonadal development of Hyriopsis cumingii. Gen Comp Endocrinol. 2022;323–4:114028.

    Article  Google Scholar 

  31. Clark BJ, Wells J, King SR, Stocco DM. The purification, cloning, and expression of a novel luteinizing hormone-induced mitochondrial protein in MA-10 mouse Leydig tumor cells. Characterization of the steroidogenic acute regulatory protein (StAR). J Biol Chem. 1994;269(45):28314–22.

    Article  CAS  PubMed  Google Scholar 

  32. Thongbuakaew T, Suwansa-Ard S, Chaiyamoon A, Cummins SF, Sobhon P. Sex steroids and steroidogenesis-related genes in the sea cucumber, Holothuria scabra and their potential role in gonad maturation. Sci Rep. 2021;11(1):2194–2194.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Yan S, Xu M, Xie J, Liao X, Liu M, Wang S, Fan S, Dong Z. Characterization, expression profiling, and estradiol response analysis of DMRT3 and FOXL2 in clam Cyclina sinensis. FRONT MAR SCI. 2023;10:1120015.

    Article  Google Scholar 

  34. Liu M, Ni H, Rong Z, Wang Z, Yan S, Liao X, Dong Z. Gonad transcriptome analysis reveals the differences in gene expression related to sex-biased and reproduction of clam Cyclina sinensis. Front Mar Sci. 2023;9:1110587.

    Article  Google Scholar 

  35. He C, Wang K, Gao Y, Wang C, Li L, Liao Y, Hu K, Liang M. Roles of noncoding RNA in reproduction. Front Renet. 2021;12:12. https://doi.org/10.3389/fgene.2021.777510.

    Article  CAS  Google Scholar 

  36. Mcbride D, Carré W, Sontakke SD, Hogg CO, Clinton M, McBride D, Sontakke SD, Hogg CO, Law A, Donadeu FX. Identification of miRNAs associated with the follicular-luteal transition in the ruminant ovary. Reproduction. 2012;144(2):221–33. https://doi.org/10.1530/REP-12-0025.

    Article  CAS  PubMed  Google Scholar 

  37. Miao X, Luo Q, Zhao H, Qin X. Genome-wide analysis of miRNAs in the ovaries of Jining Grey and Laiwu Black goats to explore the regulation of fecundity. Sci Rep-UK. 2016;6:37983.

    Article  CAS  Google Scholar 

  38. Xu G, Zhang H, Li X, Hu J, Sun S, Xu G, Zhang H, Li X, Hu J, Yang G, Sun S. Genome-wide differential expression profiling of ovarian circRNAs associated with litter size in Pigs. Front Genet. 2019;10:1010.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Song YN, Shi LL, Liu ZQ, Qiu GF. Global analysis of the ovarian microRNA transcriptome: implication for miR-2 and miR-133 regulation of oocyte meiosis in the Chinese mitten crab, Eriocheir sinensis (Crustacea: Decapoda). BMC Genomics. 2014;15(1):547.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Weijia W, Zhaofang H, Wanbo L, Zhiyong W. The identification and analysis of long noncoding RNA in testes and ovaries of the large yellow croaker (Larimichthys crocea). J Fish Sci China. 2019;26(05):852–60.

    Google Scholar 

  41. Pasquinelli AE. MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship. Nat Rev Genet. 2012;13(4):271–82. https://doi.org/10.1038/nrg3162.

    Article  CAS  PubMed  Google Scholar 

  42. Wang J, Ren Q, Hua L, Chen J, Zhang J, Bai H, Li H, Xu B, Shi Z, Cao H, Wang J, Ren Q, Hua L, Chen J, Zhang J, Bai H, Li H, Xu B, Shi Z, Cao H, Xing B, Bai X. Comprehensive analysis of differentially expressed mRNA, lncRNA and circRNA and their ceRNA networks in the longissimus dorsi muscle of two different pig breeds. Int J Mol Sci. 2019;20(5):1107.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Conley AJ. The role of cytochrome P450 17 alpha-hydroxylase and 3 beta-hydroxysteroid dehydrogenase in the integration of gonadal and adrenal steroidogenesis via the delta 5 and delta 4 pathways of steroidogenesis in mammals. Biol Reprod. 1997;56(4):789.

    Article  CAS  PubMed  Google Scholar 

  44. Bai C, Xin L, Rosani U, Wu B, Wang Q, Duan X, Liu Z, Wang C. Chromosomal-level assembly of the blood clam, Scapharca (Anadara) broughtonii, using long sequence reads and Hi-C. Gigascience. 2019;8(7):giz067.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Kong L, Yong Z, Ye ZQ, Liu XQ, Ge G. ACPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35(Web Server issue):W345–9.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, Liu Y, Chen R, Zhao Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;17:e166.

    Article  Google Scholar 

  47. Wang L, Park HJ, Dasari S, Wang S, Kocher JeP, Li W. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 2013;41:e74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Mistry J, Bateman A, Finn RD. Predicting active site residue annotations in the Pfam database. BMC Bioinformatics. 2007;8:298.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Yuan G, Wang J, Zhao F. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16(1):4.

    Article  Google Scholar 

  50. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. 2013.

  51. Friedlnder MR, Mackowiak SD, Li N, Chen W, Nikolaus R. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;1:37–52.

    Article  Google Scholar 

  52. Betel D, Wilson M, Gabow A, Marks DS, Sander C. The microRNA.org resource: targets and expression. Nucleic Acids Res. 2008;36(Database issue):D149.

    CAS  PubMed  Google Scholar 

  53. Lewis BP, Shih I, Jones-Rhoades MW, Bartel DP, Burge CB. Prediction of mammalian microRNA targets. Cell. 2003;115(7):787–98.

    Article  CAS  PubMed  Google Scholar 

  54. Minoru K, Miho F, Yoko S, Masayuki K, Mari I. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2022;51:D587–92.

  55. Li M, Zhu L, Zhou C, Sun S, Fan Y, Zhuang Z. Molecular characterization and expression of a novel big defensin (Sb-BDef1) from ark shell, Scapharca broughtonii. Fish Shellfish Immun. 2012;33(5):1167–73.

    Article  CAS  Google Scholar 

  56. Bai C, Zhang X, Venier P, Gu L, Li Y, Wang C, Xin L, Rosani U. Paired miRNA and RNA sequencing provides a first insight into molecular defense mechanisms of Scapharca broughtonii during ostreid herpesvirus-1 infection. Fish Shellfish Immun. 2022;122:225–33.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We are grateful to the Biomarker technologies Co., Ltd. for assisting in the whole-transcriptome sequencing.

Funding

This research was funded by Natural Science Foundation of Shandong Province (ZR2022MC095) and the Project of Yellow River Fisheries Resources and Environment Investigation from the MARA, P. R. China.

Author information

Authors and Affiliations

Authors

Contributions

WW, BW and ZL conceived the project and designed the experiment. WW and BW performed the experiment. WW, XS, LZ, TY and YZ analyzed the data. WW, WX and SZ made substantial contributions to the interpretation of data. WW wrote the draft manuscript. All authors contributed to the production of the revised manuscript. All authors have read and approved the final manuscript for publication, and agree to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Corresponding author

Correspondence to Biao Wu.

Ethics declarations

Ethics approval and consent to participate

The animal study protocol was approved by the Institutional animal care and use committee of Yellow Sea Fisheries Research Institute, CAFS (protocol code YSFRI-2023002).

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, W., Wu, B., Liu, Z. et al. Comprehensive analysis on the regulation of differentially expressed of mRNA and ncRNA in different ovarian stages of ark shell Scapharca broughtonii. BMC Genomics 24, 563 (2023). https://doi.org/10.1186/s12864-023-09648-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09648-z

Keywords