Skip to main content

Circular RNA expression profiles and CircSnd1-miR-135b/c-foxl2 axis analysis in gonadal differentiation of protogynous hermaphroditic ricefield eel Monopterus albus



The expression and biological functions of circular RNAs (circRNAs) in reproductive organs have been extensively reported. However, it is still unclear whether circRNAs are involved in sex change. To this end, RNA sequencing (RNA-seq) was performed in gonads at 5 sexual stages (ovary, early intersexual stage gonad, middle intersexual stage gonad, late intersexual stage gonad, and testis) of ricefield eel, and the expression profiles and potential functions of circRNAs were studied.


Seven hundred twenty-one circRNAs were identified, and the expression levels of 10 circRNAs were verified by quantitative real-time PCR (qRT–PCR) and found to be in accordance with the RNA-seq data, suggesting that the RNA-seq data were reliable. Then, the sequence length, category, sequence composition and the relationship between the parent genes of the circRNAs were explored. A total of 147 circRNAs were differentially expressed in the sex change process, and GO and KEGG analyses revealed that some differentially expressed (such as novel_circ_0000659, novel_circ_0004005 and novel_circ_0005865) circRNAs were closely involved in sex change. Furthermore, expression pattern analysis demonstrated that both circSnd1 and foxl2 were downregulated in the process of sex change, which was contrary to mal-miR-135b. Finally, dual-luciferase reporter assay and RNA immunoprecipitation showed that circSnd1 and foxl2 can combine with mal-miR-135b and mal-miR-135c. These data revealed that circSnd1 regulates foxl2 expression in the sex change of ricefield eel by acting as a sponge of mal-miR-135b/c.


Our results are the first to demonstrate that circRNAs have potential effects on sex change in ricefield eel; and circSnd1 could regulate foxl2 expression in the sex change of ricefield eel by acting as a sponge of mal-miR-135b/c. These data will be useful for enhancing our understanding of sequential hermaphroditism and sex change in ricefield eel or other teleosts.

Peer Review reports


Circular RNAs (circRNAs) are a class of endogenous noncoding RNAs characterized by covalently closed continuous loops, single strands and no free 3' or 5' terminus [1, 2]. In the past few years, the use of high-throughput sequencing technologies has demonstrated the widespread existence of circRNAs in nematodes [3, 4], zebrafish [5], fruitflies [6], mice [4, 7], pigs [8] and humans [4, 9]. Although the function of most circRNAs is largely unclear, current studies have revealed that some circRNAs are important modulators that regulate gene expression at multiple levels [4, 10,11,12,13,14,15,16,17].

CircRNAs have important biological functions during the development of the reproductive system [18,19,20,21,22]. Expression pattern analysis of many circRNAs in germ cells has shown typical developmental stages or tissue-cell specific characteristics, and their biological functions may be related to follicular development, ovarian senescence and spermatogenesis [23,24,25,26]. For example, circINHA has been shown to promote pig granulosa cell proliferation and inhibit granulosa cell apoptosis via CTGF as a competing endogenous RNA (ceRNA) that directly binds miR-10a-5p [27]. During spermatogenesis in mice, a large number of circRNAs with open reading frames are synthesized and then translated into proteins to compensate for the shortage of proteins caused by uncoupling of transcription and translation [28]. In medaka (Oryzias latipes), circ880 may be related to gonadal development by combining with miR-375-3p [22]. Despite these sporadic reports, the functional role and molecular mechanism of circRNAs in the reproductive system is still in its infancy.

Foxl2, a female sex-associated gene, is principally expressed in granulosa cells and important for ovarian development and maintenance [29,30,31,32,33,34,35]. In mice, primary follicles are not formed and cannot complete maturation when foxl2 function is lost, which may be caused by granulosa cells being reprogrammed into testis-specific Sertoli-like cells in a cell autonomous and oocyte-independent manner [35, 36]. Mechanistically, foxl2 induces XX male sex reversal mainly by activating testis-specific genes (notably Dmrt1 and Sox9) and repressing ovary-specific genes (notably Cyp19) [35, 37,38,39]. Research in olive flounder indicated that the expression of foxl2 may be repressed via maintenance of DNA methylation [40]. In goats, foxl2 mRNA has been detected in testes at several developmental stages, but no FOXL2 protein has been detected, suggesting that foxl2 undergoes translational or posttranscriptional regulation [41]. Lines of evidence have shown that the expression of foxl2 is significantly downregulated by miRNAs in humans and mice [42,43,44,45]. All of these results suggested that foxl2 undergoes transcriptional or posttranscriptional regulation.

Ricefield eel is a protogynous hermaphroditic fish that changes from female to male through the intersex stage during its lifetime. Ricefield eel has been gradually recognized as a new model species for the study of biological development and evolution due to its sex change and small genome [46,47,48]. The expression level of foxl2 in the ovary of ricefield eel has been shown to be high before sex change but to decrease sharply after the gonad developed into the ovotestis and testis [49,50,51]. Downregulation of foxl2 has been thought to be responsible for gonadal differentiation of ricefield eel, but what and how foxl2 expression declines is still unclear [51]. As important modulators involved in the transcriptional and posttranscriptional regulation of gene expression and miRNA sponges, circRNAs may be involved in the decreased expression of foxl2 in the sex change of ricefield eel. However, there are no reports about circRNAs in ricefield eel to date. In the present study, we first performed transcriptome analysis in ricefield eel gonads to identify the potential circRNAs involved in ricefield eel sex change. Then, bioinformatics was used to analyse the function of differentially expressed circRNAs and their interaction with miRNAs. Finally, a dual luciferase reporter assay and RNA immunoprecipitation (RIP) were used to demonstrate the ceRNA relationship between circSnd1, foxl2, mal-miR-135b and mal-miR-135c. To our knowledge, this study is the first to analyse the expression profile of circRNAs in the gonads of fish with sex change, providing insights into the identification of novel targets for sex change in ricefield eel.


Preliminary analysis of circRNA sequencing

Fifteen circRNA libraries from the gonads of ricefield eel at different developmental stages were constructed and sequenced, and the overview data for each library are shown in Additional file 1: Table S1. Briefly, the most and least raw and clean reads were obtained from IM1 (IM: middle intersexual gonad) and OV3 (OV: ovary), respectively. Clean base data of all libraries were > 15Gb, with Q20 > 97%, Q30 > 94%, error rate < 0.04%, GC content between 46% and 48%, and the ratio of clean reads mapped to the unique genomes was 76.46%-85.26%.

After rigorous selection, 721 circRNAs were identified from the 15 libraries, and all the circRNAs were novel because this was the first study to identify circRNAs in ricefield eel. The sequences and junction sequences of all the circRNAs are listed in Additional file 2: Data Base. Size distribution analysis revealed that the length of the 719 circRNAs ranged from 40-800 nt, the average length was 307 nt, and the majority (65.37%) ranged from 201-400 nt (Fig. 1A). Specifically, the lengths of novel_circ_0002103 (1619 nt) and novel_circ_0004589 (1480 nt) were significantly longer than those of the others. Statistical analysis of the parental genes showed that the 721 circRNAs were generated from 499 precursor mRNAs by back-splicing, and 80.36% of the precursor mRNAs produced one circRNA (Fig. 1B). Further statistics showed that all the circRNAs were exonic, intronic and intergenic, of which the majority were exonic (Fig. 1C). The numbers of exons, introns and intragenic regions contained in the circRNAs are also shown in Fig. 1D-F.

Fig. 1
figure 1

Characterization and classification of circRNAs. A CircRNA length and the number of circRNAs in the corresponding length range. B The relationship between the paternal genes and the circRNAs. “1” means that one paternal gene corresponds to one circRNA, “2” means that one paternal gene corresponds to two circRNAs, and so on. C CircRNA category and the number and proportion of circRNAs in the corresponding category. One circRNA contains one or more exons (D), introns (E), and intergene regions (F)

Analysis of differentially expressed circRNAs

The expression patterns of 147 differentially expressed (DE) circRNAs during the development of gonads are shown in the heatmap (Additional file 3: Fig. S1). In all, 82.31% of the DE circRNAs were strongly downregulated in the testis. Compared with the TE (testis), the number of DE circRNAs in the OV was the lowest (34 DE circRNAs), gradually increased, reached a peak (99 DE circRNAs) in IM vs. TE, and then gradually decreased (Fig. 2A-D). Notably, the expression of 22 circRNAs differed significantly during the whole sex change of ricefield eel (Fig. 2E and Additional file 4: Table S2).

Fig. 2
figure 2

Differentially expressed circRNAs during gonadal development. A OV vs. TE. B IE vs. TE. C IM vs. TE. D IL vs. TE. E Venn diagram of differentially expressed circRNAs. OV: ovary, IE: early intersexual gonad, IM: middle intersexual gonad, IL: late intersexual gonad, TE: testis

Functional analyses of the parent genes of the DE circRNAs

The function of circRNAs could be reflected to some extent by their parent genes. Gene ontology (GO) term analysis showed that the parent genes extensively participated in many biological processes (Fig. 3A and Additional file 5: Fig. S2A-C). Biological process terms were abundant among the 4 groups, but cellular component terms were rarely involved and even absent in the OV vs. TE group (Additional file 5: Fig. S2A). Furthermore, the number of main functional groups of the parent genes increased gradually in the process of sex change. In addition, notably, the GO term “regulation of photosynthesis” was enriched in each group, but the number of genes enriched in this GO term was the largest in the IE vs. TE groups (Fig. 3A). PDZ domain and LIM domain 2 (pdlim2), the parent gene of novel_circ_0000659, was enriched in this term.

Fig. 3
figure 3

GO and KEGG enrichment of the parent genes of the DE circRNAs. A GO enrichment of the parent genes of the DE circRNAs in the IE vs. TE groups. IE: early intersexual gonad, TE: testis, BP: biological process, CC: cellular component, MF: molecular function. B KEGG enrichment of the parent genes of the DE circRNAs

KEGG enrichment analysis was performed for all the parent genes of DE circRNAs. The results indicated that the top 20 most highly enriched (P-adj<0.05) pathways were involved in reproduction, such as the Toll-like receptor signalling pathway, progesterone-mediated oocyte maturation, mTOR signalling pathway and FoxO signalling pathway (Fig. 3B). More detailed analysis showed that 10 parent genes of 13 DE circRNAs were involved in many pathways (Additional file 6: Table S3). For example, pik3r3, LOC109952263, aldh2 and akt1 were involved in 13, 9, 13 and 14 pathways, respectively. Notably, pik3r3 and akt1 seem to be sequentially differentially expressed during the sex change of ricefield eel, pik3r3 in OV/IM vs. TE, and akt1 only in IL vs. TE. In addition, spire2 and pc were enriched in all 4 groups, but gab1, aldh2 and akt1 were enriched in only one group. A comparative analysis of circRNAs in Additional file 4: Table S2 and Additional file 6: Table S3 shows that only novel_circ_0005856, novel_circ_0003251 and novel_circ_0003252 are common.

Target miRNAs of the DE circRNAs

The results showed that 96 DE circRNAs interacted with 91 miRNAs. Most of the DE circRNAs could act as miRNA sponges, and correspondingly, one miRNA could also be bound by several circRNAs. For example, novel_circ_0004943 can interact with 10 miRNAs, and mal-miR-16a can interact with 8 circRNAs (Fig. 4).

Fig. 4
figure 4

Schematic diagram of interactions between DE circRNAs, miRNAs and mRNAs. For presentation purposes, “novel_circ_000” and  "mal-miR-" were shortened to “circ” and "miR-", respectively in the diagram

Novel_circ_0001738 (named circSnd1) and novel_circ_0004005 interact with mal-miR-135b/c and mal-miR-196a/b, respectively. In addition to mal-miR-135b/c, circSnd1 also interacts with another 7 miRNAs (mal-miR-124, mal-miR-141, mal-miR-15a, mal-miR-16a/b, mal-miR-184 and mal-miR-200a). mal-miR-135b/c has been shown to be differentially expressed during sexual change of ricefield eel (Fig. 7B). To elucidate the mechanism of ceRNA, we further predicted target mRNAs for mal-miR-135b/c and mal-miR-196a/b. The results indicated that mal-miR-135b/c and mal-miR-196 may bind to the 3' UTRs of foxl2 (Additional file 7: Fig. S3A) and akt1, respectively (Additional file 7: Fig. S3C).

Confirmation of circRNAs by PCR and RT–qPCR

We experimentally tested 10 candidate circRNAs from the gonads of ricefield eel by PCR and Sanger sequencing. Five DE and 5 non-differentially expressed (NDE) circRNAs were randomly chosen, and all 10 circRNAs were validated as expertly described. Specifically, circRNA forms were amplified in cDNA samples transcribed by random hexamer primers but not in cDNA samples transcribed by Oligo(dT)18 primers, and RNase R resistance experiments showed that circRNAs could protect against the degradation of RNase R (Fig. 5A and Additional file 8: Fig. S4). Moreover, Sanger sequencing further confirmed the circRNA sequencing and back-splicing sites.

Fig. 5
figure 5

Characterization of the predicted circRNAs. A Amplification of the circRNAs using divergent primers with cDNA. Lanes 1, 3 and 4 refer to the cDNA samples transcribed by random hexamer primers, and lane 2 refers to the cDNA samples transcribed by oligo (dT)18. The cDNA used in lane 3 was synthesized by total RNA digested with RNase R, and lane 4 was a control. : represents back-splicing sites. B Verification of circRNA expression patterns. The circRNA expression level of RNA-seq was calculated as Log50(TPM), and the RT–qPCR data were computed as the mean ± S.E.M. TPM: transcripts per million

To further verify the accuracy of the RNA-seq data, the aforementioned 10 circRNAs were quantified. The results indicated that the expression patterns of 8 out of 10 circRNAs were highly symmetrical with those of RNA-seq (Fig. 5B), and only novel_circ_0000873 and novel_circ-0001675 were not validated (data not shown). Although some subtle differences may exist, we consider them objective and acceptable. In summary, the RNA-seq results were reliable and precise.

CircSnd1 as a miRNA sponge of mal-miR-135b/c

Among the 4 DE circRNAs (novel_circ_0000659, circSnd1, novel_circ_0004005 and novel_circ_0005865) mentioned above, novel_circ_0004005 has the lowest expression level (data not shown). Based on the binding sites of miRNA and circRNA, only miR-135b/c could specifically bind to the splice site of circSnd1. Thus, circSnd1 was selected for further analyses in order to ensure that we could effectively avoid the interference of mRNA.

RNA-seq indicated that circSnd1 was derived from 6 discontinuous exon sequences of Snd1 pre-mRNA by back splicing, and bioinformatics analysis revealed that both mal-miR-135b and mal-miR-135c specifically bound to the back-splicing site of circSnd1 (Fig. 6A). Cytoplasmic and nuclear separation of ricefield eel ovaries was performed to explore the location of circSnd1. U6 and 18S were demonstrated to be mainly expressed in the nucleus and cytoplasm, respectively, and were used as reference genes to test the expression level of isolated nucleoplasm RNA and confirm that circSnd1 in the cytoplasm was detected (Fig. 6B), which further proved that circSnd1 may act as a miRNA sponge of mal-miR-135b/c.

Fig. 6
figure 6

circSnd1 interacts with mal-miR-135b/c. A The genomic loci of circSnd1 and divergent primers were designed (F: forward primer, R: reverse primer), and the expected size of the product was amplified and verified by Sanger sequencing; the cDNA samples in lanes 1-4 were the same as those in Fig. 5A. B Quantitative results of circSnd1 in the nucleus and cytoplasm. C Schematic representation of the mal-miR-135b/c target sequence within circSnd1. D, E Dual-luciferase assays for validating the interaction of mal-miR-135b/c with circSnd1. Values are the mean ± S.E.M., and the experiment was repeated at least three times. The letters above the bar chart indicated P<0.05

To further confirm this relationship, wt and mu pSI-circSnd1 plasmids were constructed (Fig. 6C). As expected, there was a significant decrease in the luciferase activity mediated by wt-pSI-circSnd1 when mal-miR-135b and mal-miR-135c mimics were cotransfected (Fig. 6D-E).

CircSnd1 was a potential modulator of foxl2 by ceRNA

To further explore the potential biological function, the expression pattern of circSnd1 during the sex change process was quantified. circSnd1 was significantly downregulated in the IL and TE groups compared to the OV, IE and IM groups (Fig. 7A). Strongly increased expression of mal-miR-135b in the testes of ricefield eel compared to the ovaries and ovotestes (Fig. 7B). Obviously, the expression pattern of circSnd1 was opposite that of mal-miR-135b, and circSnd1 was a sponge of mal-miR-135b/c (Fig. 6). Based on these data, we speculated that there may be a target mRNA of mal-miR-135b/c in ricefield eel, whose expression is downregulated during the sex change process. Foxl2 is a well-suited candidate target mRNA, because its expression is dramatically reduced in testes during the sex change (Fig. 7C). Alignment results of mal-miR-135b/c sequences with the foxl2 3' UTR also indicated that mal-miR-135b/c binding sites were detected in the foxl2 3' UTR (Additional file 7: Fig. S3A).

Fig. 7
figure 7

CircSnd1 as a potential modulator of foxl2. A-C Expression patterns of circSnd1, mal-miR-135b/c and foxl2 during the sex change. D CircSnd1, mal-miR-135b/c and foxl2 were enriched in the RNA-induced silencing complex (RISC). E Schematic representation of the mal-miR-135b/c target sequence within the foxl2 3' UTR. F Dual-luciferase assays validating the -135b/c with foxl2. Values are the mean ± SEM, and the experiment was repeated at least three time. The letters above the bar chart indicated P<0.05

Based on the expression pattern, we performed an RIP test. As a component of the RISC enzyme in Caenorhabditis elegans, Drosophila and mammals [52], snd1 was chosen as a control. The results revealed that although the miR135b enrichment multiple was low, all 5 RNAs were enriched in RISC (Fig. 7D). Remarkably, foxl2 enrichment was abnormally higher than that of other genes. We then constructed wt and mu pSI-foxl2 plasmids targeting binding sites (Fig. 7E), and a dual-luciferase reporter assay revealed that mutation of the mal-miR-135b/c binding sites completely abolished the inhibition of luciferase reporter activities (Fig. 7F). Taken together, our results strongly supported that mal-miR-135b/c binds to the foxl2 3' UTR.


CircRNAs are important regulators of gene expression

CircRNAs express in multiple species that has the variety of the mechanism and molecular functions. Recent studies have detected 104388, 96675 and 82321 circRNAs in humans, macaques and mice, respectively [53]. A total of 70186 orthologous circRNAs exhibit conserved expression and splicing patterns across species, and numerous circRNAs may be involved in adaptive immunity, blood coagulation and neurotransmitters [53]. Current studies have shown that circRNAs can regulate gene expression extensively by circRNA-miRNA-mRNA regulatory network [54]. For example, the expression of miR-7 was inhibited by injecting zebrafish embryos with plasmid that can produce circular CDR1as, and then lead to significantly reduce of midbrain size [4]. CircPIKfyve acted as a ceRNA of miR-21-3p to reduce its inhibitory effect on mavs expression, which resulted in the activation of the NF-κB/IRF3 signaling pathway, and thereby enhancing the innate antiviral responses of Miiuy croaker [55]. CircRasGEF1B was also identified in Miiuy croaker with the same mechanism [56]. In addition to regulating gene expression, some circRNAs, such as circ-ZNF609, could be translated [28, 57]. Therefore, the research of the circRNAs function in sex determination and gonadal differentiation has great potential for understanding the intricate molecular events that underpin the extraordinary sex plasticity of fish.

CircRNAs were differentially expressed during sex reversal of ricefield eel. In other fish, 5052, 3762, 3868, 975 and 2727 circRNAs were identified in grass carp [58], half-smooth tongue sole [59], zebrafish [5], large yellow croaker [60] and medaka [22], respectively. In our study, the expression pattern of circRNAs in sex reversal of ricefield eel has been analyzed. A total of 721 circRNAs were obtained from the ovary, ovotestis and testis of protogynous hermaphroditic ricefield eel, and 144 DE circRNAs were identified from sex reversal. Then, 96 DE circRNAs interacted with 91 miRNAs. Bioinformatics analysis revealed that the DE circRNAs may play the important functions in the sex change process of ricefield eel, and the function of circSnd1 was preliminarily revealed, which initiated the dawn of the study of circRNAs in fish gonads and sex change in fish. Evidence from mammals has demonstrated that circRNAs act as ceRNAs to regulate sex determination or gonadal differentiation cues [10, 12, 27, 61]. The above results suggest that DE circRNAs may play a key role in sex reversal of ricefield eel.

DE circRNAs and their parent genes involved in regulation of sex change

GO and KEGG function enrichment of DE circRNAs were analyzed, the results revealed that some parent genes of circRNAs may be involved in sex change of ricefield eel. For example, the parent genes of novel_circ_0000659 (pdlim2), novel_circ_0004005 (akt1), and spire2 (novel_circ_0005865).

The parent gene pdlim2 may be a potential photosensitive gene in ricefield eel that closely involved in reproduction. In our study, pdlim2 was enriched in the GO term of “regulation of photosynthesis”. Pdlim2 is characterized by a PDZ domain (protein–protein interaction domains) in the N-terminus and an LIM domain in the C-terminus, and played a very important biological role in organ development by actin anchorage [62]. Plidm2 is strongly expressed in the corneal epithelium [63], high expression was detected in the inner and outer layers and lenses of zebrafish [64]. Opsins are the core of vision and are expressed primarily in ciliary photoreceptor cells, such as double cones and long and short single cones [65]. The consistent expression pattern of pdlim2 and opsins led us to make the bold assumption that Pdlim2 has a function similar to that of opsins. As we known, light is one of the important environmental factors in the reproductive process for seasonally breeding animals [66]. Light can affect gonadal development and sex hormone secretion in fish. In California grunion (Leuresthes tenuis, Ayres), a higher proportion of females were found at long photoperiod, while short photoperiod treatment produced the most male-biased ratios [67]. Furthermore, the irradiation with light of a specific wavelength can trigger female-to-male sex reversal of medaka [68]. Pdlim2 is an insulin-like growth factor-I-regulated protein [69], a study in juvenile rainbow trout (Oncorhynchus mykiss) revealed that a long photoperiod could upregulate IGF-I levels in plasma [70]. Therefore, it is suggested that pdlim2, as a potential photosensitive molecule, may mediate a photoperiodic pathway and activate the sex change process of ricefield eel. Moreover, the function of novel_circ_0000659 in the photoperiodic signalling pathway deserves attention.

Given the important role of the PI3K-Akt signalling pathway in ovarian development [71], many DE genes were enriched in the PI3K-Akt signalling pathway during ricefield eel gonadal development from ovary to ovotestis [72], and a recent study also showed that the PI3K-Akt signalling pathway was regulated by miR-430s, resulting in the influence of steroidogenesis and sex differentiation of ricefield eel [73]. We speculated that the PI3K-Akt signalling pathway plays an irreplaceable role in the sex change process of ricefield eel. Bioinformatics analysis showed that mal_miR_196a/b targeted both novel_circ_0004005 and the 3' UTR of akt1 (Fig. 4 and Additional file 7: Fig. S3B-C). Studies have demonstrated that circRNAs can alter the expression of their parent genes [13,14,15]. Taken together, our present study revealed that nove_circ_0004005 may act as a ceRNA to regulate the expression of akt1, in turn playing a potential regulatory role in the sex change process.

The parent gene spire2 is one isoform of two SPIRE protein genes (spire1 and spire2) that play a key role in reproduction. SPIRE is one member of WH2-containing actin nucleators [74, 75]. The highest expression levels of spire1 and spire2 were detected in mouse oocytes [76]. Then, the high expression of spire1 and spire2 was specifically detected in rat Sertoli and germ cells and mouse spermatocytes, respectively, but only minor expression levels of spire1 were detected in mouse testes [77,78,79]. SPIRE proteins can initiate actin polymerization by binding actin monomers to four WH2 domains in the central part of the proteins [75]. In the mouse ovary, SPIRE1 and SPIRE2 can mediate asymmetric spindle positioning by assembling an action network that serves as a substrate for spindle movement and then drive polar body extrusion by promoting assembly of the cleavage furrow [76]. In rat testes, SPIRE1 is an important regulator of Sertoli cell actin and microtubule polymerization [79]. For example, spire 1 knockdown leads to gross disruption of F-actin and microtubule organization across the seminiferous epithelium, which prevents spermatids from crossing the blood- testis barrier [79]. In the present study, only the parent gene spire2 of novel_circ_0005865 was detected and enriched in all 4 groups. On this basis, we speculate that spire2 may be an important molecule for sex changes in ricefield eel.

CircSnd1 regulates foxl2 expression by functioning as a miR-135b/c sponge

CircSnd1 and foxl2 had similar expression patterns during the sex change of ricefield eel, and both of them could bind miR-135b/c through sequence complementary pairing. It was suggested that circSnd1 may regulate foxL2 expression through ceRNA mechanism to take part in the sex change of ricefield eel. Our present study demonstrated that circSnd1 was downregulated in the sex change process of ricefield eel, which was consistent with the foxl2 expression pattern but contrary to mal-miR-135b. Mechanistically, circSnd1 may act as mal-miR-135b/c sponge to downregulate foxl2 expression. Foxl2 has been found within rather well-developed seminiferous tubules, but it has never been strongly co-expressed with Sox9 in the same cell [80]. Another study showed that Sox9 immunoreactivity appears only after induction of foxl2 deletion [37]. A similar spatiotemporal pattern was observed in the gonadal development of Oryzias luzonensis [81]. Forced expression of foxl2 in XY transgenic mice induces seminiferous tubule disorganization and the development of ovotestis-like gonads, while downregulation of foxl2 leads to female-to male sex reversal [82]. In gibel carp, foxl2a-B deficiency leads to remarkably increased expression of some testis differentiation-related genes (dmrt1s and sox9bs) in the ovary but decreased expression of some oocyte-derived factors (cyp19a1as and gdf9s), which result in ovary development arrest and sex reversal [83]. Similar results have been found in Nile tilapia [39], zebrafish [84] and olive flounder [40]. All of the above evidence shows that foxl2 is mutually exclusive with male-specific genes and may induce XX male sex reversal by repressing testis-specific genes and activating ovary-specific genes.

The expression level of foxl2 has been shown to be dramatically decreased once the ovaries develop into the ovotestis and testis [51], and the results of the present study were consistent with this finding. In ricefield eel, a study demonstrated that disruption of foxl2 could significantly decrease the expression of cyp19a1a and the level of serum E2 but significantly increase the expression of dmrt1 and the levels of serum T and 11-KT [50]. Some scholars strongly believe that foxl2 is an important candidate gene responsible for sex change in ricefield eel [50, 51]. The expression of foxl2 can be suppressed by mal-miR-430a/b/c binding to the 3' UTR of foxl2 in ricefield eel [73]. In the present study, we demonstrated that mal-miR-135b/c binds to cirSnd1 and foxl2. The sharp downregulation of circSnd1 in the sex change process and upregulation of mal-miR-135b led to an increase in mal-miR-135b/c in the free state and binding to foxl2, ultimately resulting in the downregulation of foxl2. Here, we report for the first time that circRNAs regulate foxl2 expression in ricefield eel, even in a sex change fish model, providing clues for further research on the role of circRNAs in sex determination or gonadal development.


In conclusion, our present study identified 721 circRNAs in the gonads of ricefield eel for the first time, and 144 DE cricRNAs were revealed in the process of sex change of ricefield eel. GO and KEGG uncovered a lot of circRNAs, mRNAs and signalling pathways, which were closely involved in reproduction. CircSnd1 acts as a sponge of mal-miR-135b/c to reduce foxl2, and we strongly believe that the endogenous downregulation of foxl2 may be the main factor initiating the sex change of ricefield eel. Although our results enlightened the role of circRNAs in ricefield eel sex determination, more attention needs to be paid to identify more circRNAs related to gonadal development and finally clarify the sex change of ricefield eel, and even the teleost.


Experimental eels and sampling procedure

Ricefield eels (body weight = 65.20±4.59 g, body length = 41.26±7.13 cm) were purchased from a local market in Chengdu, Sichuan Province, China. The fish were temporarily maintained in the laboratory under a natural temperature and photoperiod for 24 h. The fish were randomly selected and anaesthetized with MS-222 (100 μg/mL, Syndel, Washington, USA) for 10 minutes. The gonads were collected after the ricefield eels were dissected and divided into three parts: one was fixed in Bouin’s solution for 24 h and then stored in 75% alcohol until used for determination the gonadal developmental stage. The remaining two fragments were immediately frozen in liquid nitrogen and then stored at -80°C until use; and one fragment was sent to Novogene (Novogene, Beijing, China) for RNA extraction for “library construction and RNA-seq”, and the other fragment was used for RNA extraction for “PCR amplification and Sanger sequencing”, “real-time quantitative PCR analysis”.

Histological sections that were 5 μm thick and stained with haematoxylin-eosin were used to determine the developmental stage of ricefield eel gonads based on a previous study [85, 86]. In the present study, the ovary (OV), early intersexual gonad (IE), middle intersexual gonad (IM), late intersexual gonad (IL) and testis (TE) were used (n = 3, Fig. 8).

Fig. 8
figure 8

Results of identification of the gonad development stages of ricefield eel. A Ovary (OV), B early intersexual gonad (IE), C middle intersexual gonad (IM), D late intersexual gonad (IL), E testis (TE). CAO: cortical alveolar oocyte, PGO: primary growth stage oocyte, GR: gonadal ridge, SC: spermatocyte, ST: spermatid, BV: blood vessel, EVO: early vitellogenesis oocyte

Library construction and RNA-seq

Total RNA was extracted from the OV, IE, IM, IL and TE samples of ricefield eel using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. RNA purity and contamination were assessed using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA), integrity was detected by 1% agarose gel electrophoresis, and concentration was measured using a Qubit® RNA Assay Kit in a Qubit® 2.0 Flurometer (Life Technologies, CA, USA). The RNA was qualified when A260/A280 was between 1.8 and 2.1, the brightness of the 28S band was approximately twice that of the 18S band, and there was no 5S band. The RNA samples of each stage were three biological replicates.

For circRNA enrichment, total RNA was treated with a Ribo-zero rRNA Removal Kit (Epicentre, Madison, USA) to remove ribosomal RNAs, and then 5 μg RNA was incubated with RNase R (Epicentre, USA) according to previously described procedures to remove linear RNAs [4]. Subsequently, the sequencing libraries were generated by the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s recommendations. Briefly, divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer were used to fragment the RNA. First-strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (RNase H). Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. In the reaction buffer, dNTPs with dTTP were replaced with dUTP. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of the 3' ends of DNA fragments, NEBNext adaptors with hairpin loop structures were ligated to prepare for hybridization. To preferentially select cDNA fragments of 150~200 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then, 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. Then, PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) Primer. Finally, the products were purified (AMPure XP system), and library quality was assessed on an Agilent Bioanalyzer 2100 system. After cluster generation, the libraries were sequenced on an Illumina HiSeq 4000 platform, and 150 bp paired-end reads were generated. Library construction and sequencing were carried out by Novogene (Novogene, Beijing, China).

circRNA identification

Raw data (raw reads) in FASTQ format were first processed through in-house prescripts. Clean reads were obtained by removing reads containing adapters, reads containing poly-N and low-quality reads from the raw data. At the same time, the Q20, Q30 and GC contents of the clean data were calculated. The reference genome (M_albus_1.0) was downloaded from NCBI (Monopterus albus - Nucleotide - NCBI) directly, and paired-end clean reads were aligned to the reference genome using Bowtie. Then, the circRNAs were detected and identified using find_circ [4] and CIRI2 [87].

RNA isolation and cDNA synthesize

Total RNA was extracted from the gonads of ricefield eel using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. gDNA was removed with gDNA Eraser (Trans Biotech, Beijing, China) before reverse transcription. The cDNAs were transcribed from the total RNA (1.0 μg) of OV, IE, IM, IL and TE using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, IN, USA, Cat#: K1622) according to the manufacturer’s instructions. Two kinds of cDNA were synthesized by random hexamer primers and oligo (dT)18 primers. In addition, cDNA was synthesized by random hexamer primers from the RNA, which was treated with RNase R (Epicentre, Biotechnologies Cat#: RNR07250) or ddH2O (control) as described in a previous study [4].

PCR amplification and Sanger sequencing

To validate the circRNAs, 10 circRNAs were randomly selected for PCR confirmation. Primer 5 was used to design divergent primers (both side sequences of circRNA junctions) for each candidate circRNA (Additional file 9: Table S4), and the cDNA synthesized in the previous section of this section was used as template. In general, divergent primers were expected to amplify circRNAs from cDNA synthesized by random hexamer primers but could not amplify circRNAs from cDNA synthesized by Oligo (dT)18 primers.

PCR was performed in a 10 μL final volume containing 0.5 μL of each 10 μM primer solution, 5.0 μL of 2×Taq Master Mix (Vazyme, Nanjing, China), 3 μL of ddH2O and 1 μL of gonadal cDNA template. The PCR conditions were as follows: an initial 3 min denaturation at 95 °C; 36 amplification cycles of 0.5 min at 95 °C, 0.5 min at melting temperature (Tm, °C), and 1.0 min at 72 °C; followed by a final extension for 15 min at 72 °C. PCR products were detected by 3% agarose gel electrophoresis, target products from cDNA (by random hexamer primer) were ligated into T-Vector pMDTM19 (TaKaRa, Dalian, China) for T-A cloning, and Sanger sequencing was performed by TsingKE Biological Technology Company Limited (Chengdu, Sichuan, China).

Real-time quantitative PCR analysis

The circRNAs selected above were further verified by quantitative real-time PCR (qRT–PCR) in gonads, and the primers were the same as those used to validate the circRNAs. Total RNA (100 ng) from gonad in OV, IE, IM, IL and TE stage (5 biological replicates for each stage) was reverse transcribed with random primers by using the RevertAid First-strand cDNA Synthesis Kit (Thermo Scientific, MA, USA) according to the manufacturer’s protocol, and the cDNA quality was verified by successful amplification of rpl17.

qRT–PCR was performed in a Bio–Rad CFX Connect system in a final reaction volume of 10 μl, comprising 5 μl of 2×SYBR Green Master Mix (TaKaRa, Dalian, China), 0.4 μL of each primer (10 μmol/l), 3.2 μL of nuclease-free water and 1 μL of cDNA template. The cycling parameters were 95 °C for 5 min, followed by 40 amplification cycles of 95 °C for 10 s, 59 °C for 15 s and 72 °C for 20 s. The specificity of qPCR amplification was confirmed by melting curve analysis, agarose gel electrophoresis, and sequencing of PCR products.

circRNA expression levels were quantified using a standard curve prepared with tenfold serial dilutions of the plasmid containing the corresponding DNA fragment, which ranged from 102 to 109 copies. The correlation coefficients and PCR efficiencies were not less than 0.990 and 90%, respectively. To minimize variation due to differences in cDNA loading, the expression levels of the target circRNAs were normalized to the geometric mean expression levels of rpl17 and ef1α. Target gene expression was calculated by \({\mathsf{C}}_{\mathsf{target\ circRNA}}/\sqrt{{\mathsf{C}}_{\mathsf{ef1}{\alpha}}\times {\mathsf{C}}_{\mathsf{rpl17}}}\).

Target miRNA and mRNA prediction and bioinformatic analysis

The potential miRNA binding sites of differentially expressed (DE) circRNAs were predicted by RNA22 v2 microRNA target detection, miRanda and BiBiServ2 RNAhybrid. The miRNAs of ricefield eel have not been included in miRBase, so 160 ricefield eel miRNAs that were used as libraries for target prediction. The data for mal-miR-135b/c in Fig. 7B were both cited from Gao Y et al [88], and the data for foxl2 in Fig. 7C were cited from our previous experiments [89]. Some sex-biased genes from our previous study or validated by other biologists were used as candidate mRNAs to predict miRNA binding sites. The circRNA-miRNA interaction networks were constructed by Cytoscape 3.3.0. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment [90] were performed on all host genes of DE circRNAs to determine the potential role of circRNAs.

Dual-luciferase reporter assay

The interactions between circSnd1/miR-135b/c and foxl2/miR-135b/c were measured using the Promega Dual-Luciferase System (Promega, Madison, USA). The circSnd1 sequences containing wild-type (wt) or mutant (mu, A-C and G-T mutagenesis for 4 nucleotides) miR-135b/c binding sites were synthesized, and the sequences of XhoI and NotI were synthesized upstream or downstream, respectively. For the generation of the luciferase reporter construct, the wt and mu sequences were ligated into the pSI-Check2 (cloud-seq biotech, Shanghai, China) vector, which was digested with XhoI and NotI (Invitrogen, USA) restriction enzymes. The recombinant plasmids were purified with the PureLink™ HiPure Plasmid Maxiprep Kit (Invitrogen, USA) following the manufacturer’s recommendations and then cotransfected with miR-135b/c mimics or control mimics using Lipofectamine 2000 (Invitrogen, USA) into 293T cells (Sangon Biotech, Shanghai, China). Dual-luciferase activities were tested with the Luciferase Assay Reagent (Promega, Madison, USA) 48 h later.

RNA immunoprecipitation

RIP was performed by using the Magna RIP RNA-Binding Protein Immunoprecipitation Kit (Millipore) according to the manufacturer’s instructions. Briefly, the gonads of ricefield eel were dispersed into single cells with a homogenizer and homogenized using RIP lysis buffer plus RNase inhibitor (Thermo Scientific) and proteinase inhibitor cocktail (Roche, Switzerland). After a soft freeze–thawing, 100 μL mixtures were then incubated with antibody targeting AGO2 (Invitrogen, Cat#: MA5-23515) or isotype IgG (Invitrogen, Cat#: 10400C) overnight at 4 °C, and 10 μL mixtures were taken as the input group. After washing six times with RIP wash buffer, the RNA–protein complex-associated beads were dissolved in proteinase K (Invitrogen, Cat#: 25530031) buffer, and then RNA was isolated with salt precipitation according to the manufacturer’s instructions, followed by quantitative reverse transcription PCR (RT–PCR) analysis. Primers for RT–PCR are shown in Additional file 10: Table S5. The enrichment levels were calculated by using POWER (0.5, CtAgo2 - CtIgG).

Separation of cytoplasm and nucleus

To determine the cytoplasmic and nuclear localization of circSnd1, the cytoplasmic and nuclear RNA of the fresh ovaries of ricefield eel were extracted according to the manufacturer’s recommendations for the RNA Purification Kit (Norgen Biotek, Canada) and then reversed transcribed with the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, IN, USA, Cat#: K1622) according to the manufacturer’s instructions. Small nuclear U6 and 18S were used as internal controls, and 3 biological duplicates and 5 technical duplicates were performed for each gene. Primers for RT–PCR are shown in Additional file 11: Table S6. The relative expression levels were normalized to the geometric mean expression levels of rpl17 and ef1α by \({\mathsf{C}}_{\mathsf{circSnd1}}/\sqrt{{\mathsf{C}}_{\mathsf{ef1}{{\alpha}}}\times {\mathsf{C}}_{\mathsf{rpl17}}}\).

Statistical analysis

GraphPad Prism 7 (GraphPad Software Inc., La Jolla, CA) was used to analyse the experimental data. The expression RNA levels were determined using one-way analysis of variance (ANOVA) followed by Duncan’s multiple comparisons test using SPSS 21.0 software (SPSS, Inc., Chicago, IL, USA). All the values are expressed as the mean ± S.E.M, and P < 0.05 was considered statistically significant.

Availability of data and materials

The datasets generated and analysed during the current study are available in the Sequence Read Archive ( of National Center for Biotechnology Information database with accession number SRR17761834 to SRR17761848. The link for reviewers is “”.The datasets analysed during this study are included in this published article and its supplementary information files. Please contact Zhi He ( if someone wants to request the data from this study.



Circular RNA


Competing endogenous RNAs


Monopterus albus miRNA

3' UTR:

3' untranslated region

5' UTR:

5' untranslated region


Ovary stage


Early intersexual stage


Middle intersexual stage


Late intersexual stage


Testis stage


RNA immunoprecipitation


RNA sequencing.


  1. Yang Q, Wu J, Zhao J, et al. Circular RNA expression profiles during the differentiation of mouse neural stem cells. BMC Syst Biol. 2018;12:8.

    Article  CAS  Google Scholar 

  2. Chen LL. The expanding regulatory mechanisms and cellular functions of circular RNAs. Nat Rev Mol Cell Biol. 2020;21:8.

    Article  CAS  Google Scholar 

  3. Ivanov A, Memczak S, Wyler E, et al. Analysis of Intron Sequences Reveals Hallmarks of Circular RNA Biogenesis in Animals. Cell Rep. 2015;10:2.

    Article  CAS  Google Scholar 

  4. Memczak S, Jens M, Elefsinioti A, et al. Circular RNAs are a large class of animal RNAs with regulatory potencyJ. Nature. 2013;495:7441.

    Article  CAS  Google Scholar 

  5. Shen Y, Guo X, Wang W. Identification and characterization of circular RNAs in zebrafish. FEBS Lett. 2017;591:1.

    Article  CAS  Google Scholar 

  6. Westholm JO, Miura P, Olson S, et al. Genome-wide analysis of drosophila circular RNAs reveals their structural and sequence properties and age-dependent neural accumulation. Cell Rep. 2014;9:5.

    Article  CAS  Google Scholar 

  7. Fan X, Zhang X, Wu X, et al. Single-cell RNA-seq transcriptome analysis of linear and circular RNAs in mouse preimplantation embryos. Genome Biol. 2015;16:1.

    Article  CAS  Google Scholar 

  8. Venø MT, Hansen TB, Venø ST, et al. Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development. Genome Biol. 2015;16:1.

    Article  CAS  Google Scholar 

  9. Salzman J, Gawad C, Wang PL, et al. Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types. PLoS One. 2012;7:2.

    Article  CAS  Google Scholar 

  10. Hansen TB, Jensen TI, Clausen BH, et al. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495:7441.

    Article  CAS  Google Scholar 

  11. Wu W, Wu Z, Xia Y, et al. Downregulation of circ_0132266 in chronic lymphocytic leukemia promoted cell viability through miR-337-3p PML axis. Aging. 2019;11:11.

    Google Scholar 

  12. Huang X, Wu B, Chen M, et al. Depletion of exosomal circLDLR in follicle fluid derepresses miR-1294 function and inhibits estradiol production via CYP19A1 in polycystic ovary syndrome. Aging. 2020;12:15.

    Google Scholar 

  13. Zhang Y, Zhang XO, Chen T, et al. Circular Intronic Long Noncoding RNAs. Mol Cell. 2013;51:6.

    Google Scholar 

  14. Li Z, Huang C, Bao C, et al. Exon-intron circular RNAs regulate transcription in the nucleus. Nat Struct Mol Biol. 2015;22:3.

    Google Scholar 

  15. Conn VM, Hugouvieux V, Nayak A, et al. A circRNA from SEPALLATA3 regulates splicing of its cognate mRNA through R-loop formation. Nature Plants. 2017;3:5.

    Article  CAS  Google Scholar 

  16. Zhang XO, Wang HB, Zhang Y, et al. Complementary Sequence-Mediated Exon Circularization. Cell. 2014;159:1.

    Article  CAS  Google Scholar 

  17. Ashwal FR, Meyer M, Pamudurti NR, et al. circRNA biogenesis competes with pre-mRNA splicing. Mol Cell. 2014;56:1.

    Article  CAS  Google Scholar 

  18. Talhouarne GJS, Gall JG. Lariat intronic RNAs in the cytoplasm of Xenopus tropicalis oocytes. RNA. 2014;20:9.

    Article  CAS  Google Scholar 

  19. Dong WW, Li HM, Qing XR, et al. Identification and characterization of human testis derived circular RNAs and their existence in seminal plasma. Sci Rep. 2016;6:2.

    Article  CAS  Google Scholar 

  20. Gao Y, MingliWu FY, et al. Identification and characterization of circular RNAs in Qinchuan cattle testis. Royal Society Open. Science. 2018;5(7).

  21. Tao H, Xiong Q, Zhang F, et al. Circular RNA profiling reveals chi_circ_0008219 function as microRNA sponges in pre-ovulatory ovarian follicles of goats (Capra hircus). Genomics. 2018;110:4.

    Article  CAS  Google Scholar 

  22. Zhang Y, Zhong Y, Guo S, et al. CircRNA profiling reveals circ880 functions as miR-375-3p sponge in medaka gonads. Comp Biochem Physiol Part D Genomics Proteomics. 2021;38(7).

  23. Lin X, Han M, Cheng L, et al. Expression dynamics, relationships, and transcriptional regulations of diverse transcripts in mouse spermatogenic cells. RNA Biol. 2016;13:10.

    Google Scholar 

  24. Shen M, Li T, Zhang G, et al. Dynamic expression and functional analysis of circRNA in granulosa cells during follicular development in chicken. BMC Genomics. 2019;20:1.

    Article  Google Scholar 

  25. Jia W, Xua B, Wu J. Circular RNA expression profiles of mouse ovaries during postnatal development and the function of circular RNA epidermal growth factor receptor in granulosa cells. Metabolism. 2018;185(4).

  26. Cheng J, Huang J, Yuan S, et al. Circular RNA expression profiling of human granulosa cells during maternal aging reveals novel transcripts associated with assisted reproductive technology outcomes. PLoS One. 2017;12:6.

    Google Scholar 

  27. Guoa T, Zhanga J, Yaoa W, et al. CircINHA resists granulosa cell apoptosis by upregulating CTGF as a ceRNA of miR-10a-5p in pig ovarian follicles. BBA - Gene Regul Mech. 2019;1862:10.

    Google Scholar 

  28. Tang C, Xie Y, Yu T, et al. m6A-denpendent biogenesis of circular RNAs in male germ cells. Cell Res. 2020;30:3.

    Article  Google Scholar 

  29. Cocquet J, Pailhoux E, Jaubert F, et al. Evolution and expression of FOXL2. J Med Genet. 2002;39:12.

    Article  Google Scholar 

  30. Wang DD, Zhang GR, Wei KJ, et al. Molecular identification and expression of the Foxl2 gene during gonadal sex differentiation in northern snakehead Channa argus. Fish Physiol Biochem. 2015;41:6.

    Google Scholar 

  31. Sridevi P, Senthilkumaran B. Cloning and differential expression of FOXL2 during ovarian development and recrudescence of the catfish, Clarias gariepinus. Gen Comp Endocrinol. 2011;174:3.

    Article  CAS  Google Scholar 

  32. Wang D, Kobayashi T, Zhou L, et al. Molecular cloning and gene expression of Foxl2 in the Nile tilapia, Oreochromis niloticus. Biochem Biophys Res Commun. 2004;320:1.

    Article  CAS  Google Scholar 

  33. Schmidt D, Ovitt CE, Anlag K, et al. The murine winged-helix transcription factor Foxl2 is required for granulosa cell differentiation and ovary maintenance. Developmental. 2004;131:4.

    Google Scholar 

  34. Boulanger L, Ml P, Gall L, et al. FOXL2 Is a Female Sex-Determining Gene in the Goat. Curr Biol. 2014;24(4).

  35. Uhlenhaut NH, Jakob S, Anlag K, et al. Somatic Sex Reprogramming of Adult Ovaries to Testes by FOXL2 Ablation. Cell. 2009;139:6.

  36. Uda M, Ottolenghi C, Crisponi L, et al. Foxl2 disruption causes mouse ovarian failure by pervasive blockage of follicle development. Hum Mol Genet. 2004;13:11.

    Article  CAS  Google Scholar 

  37. Veitia RA. FOXL2 versus SOX9: A lifelong “battle of the sexes”. Bioessays. 2010;32:5.

    Article  CAS  Google Scholar 

  38. Ottolenghi C, Omari S, Garcia-Ortiz JE, et al. Foxl2 is required for commitment to ovary differentiation. Hum Mol Genet. 2005;14:14.

    Article  CAS  Google Scholar 

  39. Zhang X, Li M, Ma H, et al. Mutation of foxl2 or cyp19a1a results in female to male sex reversal in XX Nile tilapia. Endocrinology. 2017;158:8.

    Google Scholar 

  40. Fan Z, Zou Y, Liang D, et al. Roles of forkhead box protein L2 (foxl2) during gonad differentiation and maintenance in a fish, the olive flounder (Paralichthys olivaceus). Reprod Fertil Dev. 2019;31:11.

    Google Scholar 

  41. Paihoux E, Vigier B, Chaffaux S, et al. A 11.7-kb deletion triggers intersexuality and polledness in goats. Nat Genet. 2001;163(4).

  42. Cho SH, An HJ, Kim KA, et al. Single nucleotide polymorphisms at miR-146a/196a2 and their primary ovarian insufficiencyrelated target gene regulation in granulosa cells. PLoS One. 2017;12:8.

    Google Scholar 

  43. Luo Y, Wu X, Ling Z, et al. microRNA133a targets Foxl2 and promotes differentiation of C2C12 into myogenic progenitor cells. DNA Cell Biol. 2015;34:1.

    CAS  Article  Google Scholar 

  44. Yu L, Chen J, Liu Y, et al. MicroRNA-937 inhibits cell proliferation and metastasis in gastric cancer cells by downregulating FOXL2. Cancer Biomark. 2018;21:1.

    CAS  Article  Google Scholar 

  45. Rosario R, Blenkiron C, Shelling AN. Comparative study of microRNA regulation on FOXL2 between adult-type and juvenile-type granulosa cell tumours in vitro. Gynecol Oncol. 2013;129:1.

    Article  CAS  Google Scholar 

  46. Zhao X, Luo M, Li Z, et al. Chromosome-scale assembly of the Monopterus genome. Gigascience. 2018;7:5.

    Google Scholar 

  47. Cheng H, Guo Y, Yu Q, et al. The rice field eel as a model system for vertebrate sexual development. Cytogenet Genome Res. 2003;101:3–4.

    Article  Google Scholar 

  48. Cheng H, He Y, Zhou R. Swamp eel (Monopterus albus). Trends Genet. 2021;37:12.

    Article  CAS  Google Scholar 

  49. He Z, Li Y, Wu Y, et al. Differentiation and morphogenesis of the ovary and expression of gonadal development-related genes in the protogynous hermaphroditic ricefield eel Monopterus albus. J Fish Biol. 2014;85:5.

    Article  CAS  Google Scholar 

  50. Feng K, Luo H, Li Y, et al. High efficient gene targeting in rice field eel Monopterus albus by transcription activator-like effector nucleases. Sci Bull. 2017;62:3.

    Article  CAS  Google Scholar 

  51. Hu Q, Guo W, Gao Y, et al. Molecular cloning and analysis of gonadal expression of Foxl2 in the ricefield eel Monopterus albus. Sci Rep. 2014;3:4.

    Google Scholar 

  52. Caudy AA, Ketting RF, Hammond SM, et al. A micrococcal nuclease homologue in RNAi effector complexes. Nature. 2003;425:6956.

    Article  CAS  Google Scholar 

  53. Ji P, Wu W, Chen S, et al. Expanded Expression Landscape and Prioritization of Circular RNAs in Mammals. Cell Rep. 2019;26:12.

    Google Scholar 

  54. Holdt LM, Kohlmaier A, Teupser D. Molecular roles and function of circular RNAs in eukaryotic cells. Cell Mol Life Sci. 2018;75:11.

    Article  CAS  Google Scholar 

  55. Su H, Chu Q, Zheng W, Chang R, et al. Circular RNA circPIKfyve acts as a sponge of miR-21-3p to enhance antiviral immunity through regulating MAVS in teleost fish. J Virol. 2021;95:8.

    Google Scholar 

  56. Chu Q, Zheng W, Su H, et al. A Highly Conserved Circular RNA circRasGEF1B Enhances Antiviral Immunity by Regulating miR-21-3p MITA Pathway in Lower Vertebrates. J Virol. 2021;95:7.

    Article  Google Scholar 

  57. Legnini I, Timoteo GD, Rossi F, et al. Circ-ZNF609 is a Circular RNA that Can be Translated and Functions in Myogenesis. Mol Cell. 2017;66:1.

    Article  CAS  Google Scholar 

  58. He L, Zhang A, Xiong L, et al. Deep Circular RNA Sequencing Provides Insights into the Mechanism Underlying Grass Carp Reovirus Infection. Int J Mol Sci. 2017;18:9.

    Google Scholar 

  59. Li J, Lv Y, Liu R, et al. Identification and characterization of a conservative W chromosome-linked circRNA in half-smooth tongue sole (Cynoglossus semilaevis) reveal its female-biased expression in immune organs. Fish Shellfish Immunol. 2018;82:3.

    Google Scholar 

  60. Xiao SX, Qiu C, Wang Z. Transcriptome-wide identification and functional investigation of circular RNA in the teleost large yellow croaker (Larimichthys crocea). Marine Genomics. 2017;32:1.

    Article  Google Scholar 

  61. Jia M, Li X, Jiang C, et al. Testis-enriched circular RNA circ-Bbs9 plays an important role in Leydig cell proliferation by regulating a CyclinD2-dependent pathway. Reprod Fertil Dev. 2019;32:4.

    Google Scholar 

  62. Velthuis AJW, Bagowski CP. PDZ and LIM Domain–Encoding Genes: Molecular Interactions and their Role in Development. Scientific World Journal. 2007;7:1.

    Article  CAS  Google Scholar 

  63. Torrado M, Senatorov VV, Trivedi R, et al. Pdlim2, a Novel PDZ–LIM Domain Protein, Interacts with α-Actinins and Filamin A. Invest Ophthalmol Vis Sci. 2004;45:11.

    Article  Google Scholar 

  64. Velthuis AJW, Ott EB, Marques IJ, et al. Gene expression patterns of the ALP family during zebrafish development. Gene Expression Patterns. 2007;7:3.

    Article  CAS  Google Scholar 

  65. Minamoto T, Shimizu I. Molecular cloning of cone opsin genes and their expression in the retina of a smelt, Ayu (Plecoglossus altivelis, Teleostei). Comp Biochem Physiol B Biochem Mol Biol. 2005;140(2).

  66. Bromage N, Porter M, Randall C. The environmental regulation of maturation in farmed finfish with special reference to the role of photoperiod and melatonin. Aquaculture. 2001;197:1.

    Article  Google Scholar 

  67. Brown EE, Baumann H, Conover DO. Temperature and photoperiod effects on sex determination in a fish. J Exp Mar Biol Ecol. 2014;461:3.

    Article  Google Scholar 

  68. Hayasaka O, Takeuchi Y, Shiozaki K, et al. Green light irradiation during sex differentiation induces female-to-male sex reversal in the medaka Oryzias latipes. Sci Rep. 2019;9(1).

  69. Loughran G, Healy NC, Kiely PA, et al. Mystique Is a New Insulin-like Growth Factor-I-regulated PDZ-LIM Domain Protein That Promotes Cell Attachment and Migration and Suppresses Anchorage-independent Growth. Mol Biol Cell. 2005;16:4.

    Article  Google Scholar 

  70. Taylor JF, Migaud H, Porter MJR, et al. Photoperiod infuences growth rate and plasma insulin-like growthf actor-I levels in juvenile rainbow trout, Oncorhynchus mykiss. Gen Comp Endocrinol. 2005;142:1–2.

    Article  CAS  Google Scholar 

  71. Makker A, Goel MM, Mahdi AA. PI3K/PTEN/Akt and TSC/mTOR signaling pathways, ovarian dysfunction, and infertility: an updat. J Mol Endocrinol. 2014;53:3.

    Article  CAS  Google Scholar 

  72. Chi W, Gao Y, Hu Q, et al. Genome-wide analysis of brain and gonad transcripts reveals changes of key sex reversal-related genes expression and signaling pathways in three stages of Monopterus albus. PLoS One. 2017;12:3.

    Google Scholar 

  73. Zhang L, Yang Q, Xu W, et al. Integrated Analysis of miR-430 on Steroidogenesis-Related Gene Expression of Larval Ricefield Eel Monopterus albus. Int J Mol Sci. 2021;22:13.

    Google Scholar 

  74. Quinlan ME, Heuser JE, Kerkhoff E, et al. Drosophila Spire is an actin nucleation factor. Nature. 2005;433:7024.

    Article  CAS  Google Scholar 

  75. Kerkhoff E. Cellular functions of the Spir actin-nucleation factors. Trends Cell Biol. 2006;16:9.

    Article  CAS  Google Scholar 

  76. Pfender S, Kuznetsov V, Pleiser S, et al. Spire-type actin nucleators cooperate with Formin-2 to drive asymmetric oocyte division. Curr Biol. 2011;21:11.

    Article  CAS  Google Scholar 

  77. Pleiser S, Rock R, Wellmann J, et al. Expression patterns of the mouse Spir-2 actin nucleator. Gene Expr Patterns. 2010;10:7–8.

    Article  CAS  Google Scholar 

  78. Schumacher N, Borawski JM, Leberfinger CB, et al. Overlapping expression pattern of the actin organizers Spir-1 and formin-2 in the developing mouse nervous system and the adult brain. Gene Expr Patterns. 2004;4:3.

    Article  CAS  Google Scholar 

  79. Wen Q, Li N, Xiao X, et al. Actin nucleator Spire 1 is a regulator of ectoplasmic specialization in the testis. Cell Death Dis. 2018;9:2.

    Article  CAS  Google Scholar 

  80. Hersmus R, Kalfa N, Leeuw B, et al. FOXL2 and SOX9 as parameters of female and male gonadal differentiation in patients with various forms of disorders of sex development (DSD). J Pathol. 2008;215(1).

  81. Nakamoto M, Muramatsu S, Yoshida S, et al. Gonadal sex differentiation and expression of Sox9a2, Dmrt1, and Foxl2 in Oryzias luzonensis. Genesis. 2009;47:5.

    Article  CAS  Google Scholar 

  82. Ottolenghi C, Pelosi E, Tran J, et al. Loss of Wnt4 and Foxl2 leads to female-to-male sex reversal extending to germ cells. Hum Mol Genet. 2007;16:23.

    Article  CAS  Google Scholar 

  83. Gan R-H, Wang Y, Li Z, et al. Functional Divergence of Multiple Duplicated Foxl2 Homeologs and Alleles in a Recurrent Polyploid Fish. Mol Biol Evol. 2021;38:5.

    Article  CAS  Google Scholar 

  84. Yang YJ, Wang Y, Li Z, et al. Sequential, Divergent, and Cooperative Requirements of Foxl2a and Foxl2b in Ovary Development and Maintenance of Zebrafish. Genetics. 2017;205:4.

    Google Scholar 

  85. He Z, Wu Y, Xie J, et al. Growth differentiation factor 9 (Gdf9) was localized in the female as well as male germ cells in a protogynous hermaphroditic teleost fish, ricefield eel Monopterus albus. Gen Comp Endocrinol. 2012;178:2.

    Article  CAS  Google Scholar 

  86. Chan STH, Phillip JG. The structure of the gonad during natural sex reversal in Monopterus albus. J Zool. 1967;151:1.

    Article  Google Scholar 

  87. Yuan G, Wang J, Zhao F. CIRI: An efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16:1.

    Article  Google Scholar 

  88. Gao Y, Guo W, Hu Q, et al. Characterization and Differential Expression Patterns of Conserved microRNAs and mRNAs in Three Genders of the Ricefield Eel (Monopterus albus). Sex Dev. 2014;8:6.

    CAS  Article  Google Scholar 

  89. He Z, Deng F, Yang D, et al. Crosstalk between sex-related genes and apoptosis signaling reveals molecular insights into sex change in a protogynous hermaphroditic teleost fish, ricefield eel Monopterus albus. Aquaculture. 2022;552:737918.

    CAS  Article  Google Scholar 

  90. Kanehisa M, Sato Y, Kawashima M, et al. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44:D1.

    Article  CAS  Google Scholar 

Download references


We thank Li Zheng, Yong Pu, Yuanyuan Jiao, Kuo Gao, Jinxing Xiong and Bolin Lai for their assistance with fish management and sample collection.


This research was supported by the National Natural Science Foundation of China [grant numbers 31972777, 2019; 31402286, 2015].

Author information

Authors and Affiliations



ZH Conceptualization, formal analysis, supervision, project administration. ZJM: Conceptualization, methodology, investigation, statistical analysis, writing draft. DYY: Conceptualization, methodology, supervision, writing, original and draft. QQC Investigation and statistical analysis. ZDH Resources and nvestigation. JXH: Investigation. FQD Methodology and investigation. QZ Investigation. JYH Investigation. LJY Resources. HJC Resources. LH Investigation and resources. XLH Conceptualization and review. WL Investigation. SYY Investigation. XBG Validation. MWZ Conceptualization and investigation. TMY Project administration, validation and editing. The authors read and approved the final manuscript.

Corresponding author

Correspondence to Taiming Yan.

Ethics declarations

Ethics approval and consent to participate

All experimental protocols involved in fishes in this study were conducted in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the Sichuan Agricultural University and ARRIVE guidelines ( All procedures and investigations were reviewed and approved by the Animal Research and Ethics Committees of Sichuan Agricultural University and performed in accordance with the guidelines of the committee (Approval No.20190031).

Consent for publication

Not applicable

Competing interests

The authors declare that they have 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 The Creative Commons Public Domain Dedication waiver ( 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

Verify currency and authenticity via CrossMark

Cite this article

He, Z., Ma, Z., Yang, D. et al. Circular RNA expression profiles and CircSnd1-miR-135b/c-foxl2 axis analysis in gonadal differentiation of protogynous hermaphroditic ricefield eel Monopterus albus. BMC Genomics 23, 552 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • circRNA
  • Monopterus albus
  • Foxl2
  • Gonadal differentiation
  • Reproductive
  • Expression pattern