CircRNA Identification and CircRNA–miRNA–mRNA Network in Cynoglossus semilaevis Sexual Size Dimorphism

Simple Summary Chinese tongue sole (Cynoglossus semilaevis) typically displays female-biased sexual size dimorphism (SSD), but its epigenetic regulatory mechanisms are poorly understood, especially the role of circRNAs. To explore the function of circRNAs in Chinese tongue sole SSD, we firstly identified differentially expressed circular RNAs (DE circRNAs) in female, male, and pseudo-male C. semilaevis. Secondly, the ceRNA network containing DE circRNAs, miRNAs, and mRNAs in the three genders was constructed. Among the ceRNA network, several circRNAs such as novel_circ_004374 and novel_circ_014597 may regulate hipk2 expression by sponging miR-130-x. It is also worth exploring whether or how novel_circ_008696 regulates setd2 by binding to novel-m0387-3p. The present study provided the cirRNA and its ceRNA network that potentially regulate C. semileavis female-biased SSD for the first time. Abstract Sexual size dimorphism (SSD), which is the sexual differences in body size, has been widely reported in various species including fishes. For Chinese tongue sole (Cynoglossus semilaevis), a flatfish exhibiting typically female-biased SSD, little is known for its epigenetic regulation mechanism, especially the role of circRNAs. Here, we identified the differently expressed abundances of circRNAs in females, males, and pseudo-males to explore the potential functions of circRNAs in Chinese tongue sole SSD. In total, 14,745 novel circRNAs were screened, among which 1461 DE circRNAs were identified from the brain, gonad, liver, and muscle in female, male, and pseudo-male individuals. The ceRNA network was subsequently constructed, including 10 circRNAs, 26 mRNAs, and 11 miRNAs. These DE mRNAs were mainly related to the mRNA surveillance pathway, metabolic pathways, and cellular senescence. Importantly, the ceRNA network has revealed that several circRNAs such as novel_circ_004374 and novel_circ_014597 may regulate homeodomain interacting protein kinase 2 (hipk2) expression by sponging miR-130-x. It is also worth exploring whether or how novel_circ_008696 regulates SET Domain Containing 2, histone lysine methyltransferase (setd2), which in turn affects the epigenetic patterns of different sexual individuals. The present study not only enriches the knowledge on the potential roles of circRNA in the physiological process, but also provides new clues for the explanation of fish SSD. In future studies, the precise function and involvement of circRNAs in female-biased SSD will require more efforts.


Introduction
Many species have evolved to exhibit female and male body size differences, such as frogs, lizards, birds, and bats, and this sexual size dimorphism (SSD) may increase their survival advantage in the wild [1][2][3][4]. It has been reported that more than 600 kinds of fish species exhibit female or male-biased SSD [5], which would cause growth disadvantages for male or female individuals, restricting in this way the sustainable development of fish aquaculture [6]. For Chinese tongue sole (Cynoglossus semilaevis), a typically femalebiased SSD flatfish in aquaculture [7], its females grow 2-4 times faster than the males [8,9]. Moreover, the environment factors including high temperatures often cause female-tomale sex reversal, and the pseudo-males exhibit physiological characteristics and growth performances similar to the males [10]. Thus, a low proportion of females has affected the work of C. semilaevis farmers and increased the cost of aquaculture [11].
Therefore, the present study aimed to illustrate the circRNAs that might regulate C. semilaevis SSD by whole transcriptomic analysis of four important tissues including brain, liver, gonad, and muscle from female, male, and pseudo-male individuals. The competing endogenous RNA (ceRNA) network of DE circRNAs was constructed. Moreover, the enriched GO terms and pathways were analyzed to reveal the potential effect of circRNA on C. semilaevis SSD.

RNA Extraction and High-Throughput Sequencing
Total RNA was extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol and RNA quality was evaluated by Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and agarose gel electrophoresis. RNAs with RIN > 7.0 were submitted for library construction. Ribosomal RNAs (RNAs) were removed to enrich mRNAs and ncRNAs by Ribo-ZeroTM Magnetic Kit (Epicentre, Madison, WI, USA). The enriched ncRNA was fragmented into short fragments using fragmentation buffer and reverse transcribed into cDNA by using VAHTS Total RNA-seq (H/M/R) Library Prep Kit for Illumina (Vazyme, Nanjing, China). The segments' purification was conducted by QIAQuick PCR Extraction Kit (QIAGEN, Venlo, the Netherlands) and the second strand degradation was performed by Uracil-N-Glycosylase. The fragment was subsequently amplified with 15 cycles, and the cDNA quality and quantity were assessed on Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and Qubit 4.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). The cDNA library was sequenced by Illumina HiSeqTM 4000 in Gene Denovo Biotechnology Co. (Guangzhou, China). Raw reads were analyzed with the fastp software (version 0.20.1) [46] to filter low-quality data and obtain clean reads. The alignment tool bowtie2 (version 2.3.2) [47] was used to align the clean reads to remove the aligned ribosome reads without allowing mismatches, and the retained unmapped reads were used for subsequent analysis. HISAT2 (version 2.1.1) software [48] was used to align the reference C. semilaevis genome (NCBI GCA_000523025.1, Cse_v1.0) to construct a reference set of splice sites, and find_circ [20] was employed to identify circRNAs. DE circRNAs were assigned as |log2fold change| ≥ 1 and p value < 0.05 by edgeR (version 3.12.1) [49] package.

circRNA-miRNA-mRNA Association Analysis
A circRNA-miRNA-mRNA regulatory network was built with DE circRNAs, DE miRNAs, and DE mRNAs. The ceRNA network between DE lncRNAs, DE miRNAs, and DE mRNAs data was analyzed in our previous study [16]. The miRNA targets were predicted using TargetScan software (version 8.0). The correlation of DE mRNA-DE miRNA or DE circRNA-DE miRNA expression was evaluated by Spearman Rank correlation coefficient (SCC) < −0.7. Pearson correlation coefficient (PCC) was used to evaluate the correlation of DE circRNA-DE mRNA expression [50]. The circRNA-mRNA pairs with PCC > 0.9 were selected as co-expressed negatively miRNA pairs. The common miRNA sponges with p values < 0.05 were set as significant. The Sankey diagram and circRNA expression heatmap were made with OmicShare tool, a free online platform for data analysis (www.omicshare.com/tools, accessed on 6 March 2022). The Venn diagrams were created using the VennDiagram (version 1.7.3) package.

RT-PCR and qRT-PCR
Total RNA was extracted from the tissues by TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. Reverse transcription was performed to synthesize cDNA with random primers using PrimeScript TM RT Reagent Kit (Takara, Shiga, Japan). Real-time reverse transcription PCR (qRT-PCR) was performed using SYBR GREEN qPCR Super Mix (Invitrogen, Carlsbad, CA, USA) on 7500Fast Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). A 20 µL qRT-PCR reaction was conducted at 95 • C for 30 s, 40 cycles at 95 • C for 5 s, 60 • C for 33 s, with the next stage at 95 • C for 15 s, 60 • C for 60 s. According to a previous method [51], primers were designed to span the cleavage site to discriminate between linear and circular RNA transcripts. Sequences of all primers used in qRT-PCR are presented in Table 1. Detection of the length of PCR products by agarose gel electrophoresis and Sanger DNA sequencing of reverse transcription PCR (RT-PCR) products was used to demonstrate the presence of circRNAs. The qRT-PCR data were analyzed by the 2 −∆∆CT method [52] and the data obtained from qPCR were performed by two-way ANOVA test using GraphPad Prism 9 (version 9.3.0). Table 1. Primer sequences in this study.

Overview of circRNA Identification in Chinese Tongue Sole
In total, 14,745 novel circRNAs were identified from the brain, liver, muscle, and gonad tissues of male, female, and pseudo-male Chinese tongue soles. Over 70% of the reads were obtained by sequencing in the exon region ( Figure 1A), and the majority of cirRNA types were annot_exons (67%) and exon_intron (16%) ( Figure 1B). The length of most circRNAs (85%) ranged from 100 nt to 2500 nt and at 11% from 2501 nt to 10,000 nt ( Figure 1C). Most of the chromosomes, except for the W chromosome, produce more than 400 kinds of circRNAs ( Figure 1D).

Differentially Expressed circRNAs
Among all screened circRNAs, we identified 1461 differentially expressed circRNAs (DE circRNAs) through pairwise comparisons of four different tissues in female, male, and pseudo-male fish ( Figure 2). DE circRNAs were mainly identified in gonads, with 1187 differentially expressed circRNAs (Figure 2A). In the Venn diagram of the gonad DE circRNAs, there are 561 genes in the largest intersection ( Figure 2B), which are the differential genes shared by males and pseudo males among females, indicating that males and pseudo males have similar circRNA expression patterns.

Verification of Differentially Expressed circRNAs
To verify the reliability of circRNAs derived from RNA-seq, we selected 10 differentially expressed circRNAs (novel_circ_001329, novel_circ_001538, novel_circ_004002, novel_circ_004374, novel_circ_004559, novel_circ_004830, novel_circ_012476, novel_circ_012536, novel_circ_012863, novel_circ_013312), and amplified by RT-PCR primers ( Table 1). The results of agarose gel electrophoresis showed the length of the PCR products ( Figure 3A), as expected from transcriptome sequencing. The PCR products were sequenced by Sanger sequencing to confirm the junction sites, as well as upstream and downstream sequences around the splice junctions, shown in Figure 3B. These all support the presence of circRNAs.

Differentially Expressed circRNAs
Among all screened circRNAs, we identified 1,461 differentially expressed circRNAs (DE circRNAs) through pairwise comparisons of four different tissues in female, male, and pseudo-male fish ( Figure 2). DE circRNAs were mainly identified in gonads, with 1,187 differentially expressed circRNAs (Figure 2A). In the Venn diagram of the gonad DE circRNAs, there are 561 genes in the largest intersection ( Figure 2B), which are the differential genes shared by males and pseudo males among females, indicating that males and pseudo males have similar circRNA expression patterns.

Verification of Differentially Expressed circRNAs
To verify the reliability of circRNAs derived from RNA-seq, we selected 10 differentially expressed circRNAs (novel_circ_001329, novel_circ_001538, novel_circ_004002, novel_circ_004374, novel_circ_004559, novel_circ_004830, novel_circ_012476, novel_circ_012536, novel_circ_012863, novel_circ_013312), and amplified by RT-PCR primers ( Table 1). The results of agarose gel electrophoresis showed the length of the PCR products ( Figure 3A), as expected from transcriptome sequencing. The PCR products were sequenced by Sanger sequencing to confirm the junction sites, as well as upstream and downstream sequences around the splice junctions, shown in Figure 3B. These all support the presence of circRNAs.  Eight circRNAs including novel_circ_012476, novel_circ_004002, novel_circ_004374, novel_circ_013312, novel_circ_004688, novel_circ_012476, and novel_circ_013312 were selected for qRT-PCR validation. As a result, the expression trends of these cirRNAs by qRT-PCR corresponded well with the patterns of high-throughput sequencing (Figure 4). Eight circRNAs including novel_circ_012476, novel_circ_004002, novel_circ_004374, novel_circ_013312, novel_circ_004688, novel_circ_012476, and novel_circ_013312 were selected for qRT-PCR validation. As a result, the expression trends of these cirRNAs by qRT-PCR corresponded well with the patterns of high-throughput sequencing (Figure 4).

Differentially Expressed circRNAs ceRNA Network
CircRNAs may bind to sequence-complementary miRNAs and affect the post-transcriptional regulation of miRNAs on their target genes. To assess the underlying functions of DE circRNAs as ceRNAs, we constructed the ceRNA network for the differentially expressed circRNAs ( Figure 5A). A total of 10 circular RNAs competed with 26 mRNAs for binding to 11 miRNAs, of which 6 circular RNAs (novel_circ_004374, novel_circ_004559, novel_circ_007263, novel_circ_008696, novel_circ_009071, and novel_circ_014597) were DE circular RNAs in the gonad ( Figure 5). Novel_circ_009071 and novel_circ_013520 were highly expressed in pseudo-male brains ( Figure 5B). Novel_circ_004559 was highly expressed in female muscle. Novel_circ_007263 were was expressed in male liver.

Differentially Expressed circRNAs ceRNA Network
CircRNAs may bind to sequence-complementary miRNAs and affect the posttranscriptional regulation of miRNAs on their target genes. To assess the underlying functions of DE circRNAs as ceRNAs, we constructed the ceRNA network for the differentially expressed circRNAs ( Figure 5A). A total of 10 circular RNAs competed with 26 mRNAs for binding to 11 miRNAs, of which 6 circular RNAs (novel_circ_004374, novel_circ_004559, novel_circ_007263, novel_circ_008696, novel_circ_009071, and novel_circ_014597) were DE circular RNAs in the gonad ( Figure 5). Novel_circ_009071 and novel_circ_013520 were highly expressed in pseudo-male brains ( Figure 5B). Novel_circ_004559 was highly expressed in female muscle. Novel_circ_007263 were was expressed in male liver.

Biological Functions of ceRNAs-Relevant mRNAs
In order to predict the potential functions of DE circRNAs, GO and KEGG pathway analysis for the ceRNAs-relevant mRNAs were performed. GO analysis of ceRNAs-associated mRNAs revealed that they were represented in metabolic, development, and others biological processes ( Figure S2A). KEGG pathway enrichment showed that the ceRNAsassociated mRNAs were involved in cellular senescence, mRNA surveillance pathway, metabolic pathways, and other pathways ( Figure S2B).

Discussion
Chinese tongue soles have a female heterogametic (WZ female, ZZ male) sex chromosome system and show obviously female-biased SSD [15]. As an important regulator for gene expression [53], circRNAs have been identified to be implicated in myogenesis, skeletal muscle and fat tissue growth, follicle development, and sperm motility [54][55][56][57]. Herein, we identified the different expression abundances of circRNAs in males, females, and pseudo males to examine the potential roles of circRNAs in C. semilaevis SSD.

Biological Functions of ceRNAs-Relevant mRNAs
In order to predict the potential functions of DE circRNAs, GO and KEGG pathway analysis for the ceRNAs-relevant mRNAs were performed. GO analysis of ceRNAsassociated mRNAs revealed that they were represented in metabolic, development, and others biological processes ( Figure S2A). KEGG pathway enrichment showed that the ceRNAs-associated mRNAs were involved in cellular senescence, mRNA surveillance pathway, metabolic pathways, and other pathways ( Figure S2B).

Discussion
Chinese tongue soles have a female heterogametic (WZ female, ZZ male) sex chromosome system and show obviously female-biased SSD [15]. As an important regulator for gene expression [53], circRNAs have been identified to be implicated in myogenesis, skeletal muscle and fat tissue growth, follicle development, and sperm motility [54][55][56][57]. Herein, we identified the different expression abundances of circRNAs in males, females, and pseudo males to examine the potential roles of circRNAs in C. semilaevis SSD.
In our study, a total of 14,745 novel circRNAs were identified, with most circRNAs derived from exon or exon_intron, with the exception of only 4% intergenic circRNAs ( Figure 1B). From the perspective of chromosome distribution, the number of circRNAs in the W chromosome is fewer than in the Z and autosomal chromosomes ( Figure 1D). Given that most circRNAs are generally produced from pre-mRNA back-splicing [24], the reason for there being few circRNAs in the W chromosome may be attributed to the fact that the W chromosome of C. semilaevis has the least number of functional genes (317), accounting for 17-34% of other chromosomes [9]. As in previous studies [43,[57][58][59], the C. semilaevis circR-NAs exhibited obvious tissue and sex-specific expression patterns, with gonads showing the largest number of differentially expressed circRNAs in different genders.
The lengths of most circRNAs (85%) ranged from 100 nt to 2500 nt in the present study, which is similar to grass carp circRNAs [60], while 78% of the circRNAs identified in the Japanese flounder were longer than 5000 nt [37]. In mammals, the length of the majority of circRNAs is less than 3500 nt in rats [58] and longer than 3000 nt in pigs [55]. These data suggest that the lengths of circRNAs might be affected by factors such as species, tissues, and the length of the gene itself.
CircRNAs have been involved in cellular and physiological regulation by decoying miRNAs [27][28][29]61]. The ceRNA network showed that 10 circRNAs competed with 26 mR-NAs for binding to 11 miRNAs. The enrichment of ceRNA-relevant mRNAs showed that cellular senescence, mRNA surveillance pathway, and metabolic pathways were involved. Noticeably, cellular senescence was also identified from the brown module, the negative SSD-related module in previous analyses based on the same branch of samples [15].
Homeodomain interacting protein kinase 2 (hipk2), a gene derived from a cellular senescence pathway, was highly expressed in male and pseudo-male gonads. Hipk2 is involved in transcription control, cellular apoptosis, and cell cycle regulation [62,63]. Impor-tantly, studies have reported that hipk2 may inhibit cell growth [64]. In the present study, novel_circ_004374 and novel_circ_014597 may modulate hipk2 by sponging to the miR-130. Interestingly, more than 30 lncRNAs were also predicted to be involved in hipk2 regulation by sponging miR-130 [16]. The detailed molecular mechanism needs further investigation.
Novel_circ_008696 was found to compete with SET Domain Containing 2, histone lysine methyltransferase (setd2), for binding to novel-m0387-3p. Setd2 is the main enzyme generating H3K36me3 (Histone 3 lysine 36 trimethylation), which provides docking sites for many chromatin regulators including DNA (cytosine-5)-methyltransferase 3A/B (dnmt3a/b) [65,66]. Interestingly, the present study and our previous transcriptome data have revealed that setd2 and dnmt3a/b exhibited male and pseudo-male biased expression levels in the gonad [15]. In the present study, novel_circ_008696, a competitor for setd2, was also expressed in higher levels in male and pseudo-male gonads. The overexpression of setd2 inhibited osteosarcoma cell growth by suppressing Wnt/β-catenin signaling [67]. Its role in growth regulation remains unknown. Thus, it is interesting to study whether or how novel_circ_008696 regulates setd2, which in turn affects the epigenetic levels in different sexual individuals. This information would be valuable for understanding the roles of non-coding RNAs in female-biased SSD.
In mammals, birds, and teleosts, it was reported that circRNA plays important roles in muscle growth regulation [41,[68][69][70]. Moreover, there are also studies addressing the involvement of circRNAs in the ovary development and sex determination for birds and teleosts [43,44,71]. Therefore, the present study not only enriches the potential roles of circRNA in physiological processes, but it also provides new clues for the explanation of fish SSD. Nevertheless, the circRNA-miRNA-mRNA regulatory networks require more research to elucidate the ceRNA mechanisms in Chinese tongue sole SSD.

Conclusions
The present study proved that cirRNA and its ceRNA network potentially regulate C. semileavis female-biased SSD for the first time. A total of 1461 DE circRNAs have been identified from the brain, gonad, liver, and muscle in female, male, and pseudo-male individuals. Importantly, the ceRNA network has revealed that several circRNAs such as novel_circ_004374 and novel_circ_014597 may regulate hipk2 expression by sponging miR-130-x. In addition, whether or how novel_circ_008696 regulates setd2 by binding to novel-m0387-3p was addressed. In future studies, the examination of the precise function and involvement of circRNAs in female-biased SSD requires more efforts. Only then will we be able to utilize a circRNA-or miRNA-based genome editing tool for sustainable aquaculture production.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/biology11101451/s1, Figure S1. The sex identification for male and pseudo male individuals; Figure S2  Informed Consent Statement: Not applicable.

Data Availability Statement:
The datasets presented in this study can be found in the NCBI SRA repository with accession number PRJNA743138 (https://www.ncbi.nlm.nih.gov/bioproject/ PRJNA743138/) (Accessed on 10 September 2021).

Conflicts of Interest:
The authors declare no conflict of interest.