Sex-Dependent RNA Editing and N6-adenosine RNA Methylation Profiling in the Gonads of a Fish, the Olive Flounder (Paralichthys olivaceus)

Introduction

RNA transfers genetic information from DNA to protein and occupies a central position in all cellular processes. RNA modifications provide a rapid direct way to regulate gene expression. To date, more than 150 chemical modifications have been identified in eukaryotic mRNA (Helm and Motorin, 2017). Among these RNA modifications, adenosine-to-inosine (A-to-I) editing and N6-methyladenosines (m6A) are two of the most abundant. Both of these RNA modifications occur on A bases, but there is a negative correlation between them (Xiang et al., 2018). RNA editing is an enzyme-mediated post/cotranscriptional mechanism that alters the nucleotide composition of an RNA sequence without affecting the encoding DNA sequence (Farajollahi and Maas, 2010; Nishikura, 2016). This modification leads to differences between the transcript sequence and its template DNA, thereby recoding the open reading frame, affecting alternative splicing, and influencing RNA stability (Heale et al., 2010; Nishikura, 2016; Tajaddod et al., 2016). A-to-I editing is the predominant type of RNA editing and is catalyzed by Adenosine Deaminases Acting on RNA (ADARs) (Nishikura, 2016). ADAR enzymes (ADAR1, ADAR2, and ADAR3) bind to double-stranded RNAs (dsRNAs) and deaminate A to I. As the inosine (I) is generally read as guanosine (G) with the cellular machinery in the translation event, A-to-I substitution leads to A-to-G transition in the edited substrate (Nishikura, 2010). Thus, A-to-I editing is also named A-to-G editing (Bazak et al., 2014; Zhang et al., 2019). ADARs are essential for normal development, and knockout or homozygous null mutations in ADAR genes cause lethality in mice (Mus musculus) (Higuchi et al., 2000; Wang et al., 2000). RNA editing studies are being widely implemented in mammals by comparing matched RNA and DNA sequencing data (Tan et al., 2017), and most of these studies focus on the associations with diseases ranging from cancer to neurological diseases (Galeano et al., 2012; Slotkin and Nishikura, 2013; Wang et al., 2017). Some RNA editing events are tissue specific (Tao et al., 2012), but the available knowledge about RNA editing in gonads is very limited. In mouse testis, adar1 and adar2 genes are expressed in both Sertoli cells and the germline in a cell-type-dependent manner during germ cell development (Snyder et al., 2017), and 2012 genes exhibiting editing events in germ cells during spermatogenesis have been detected (Wang et al., 2019). However, the functions of RNA editing in gonads are still unknown.

N6-methyladenosine (m6A) was the first RNA modification shown to be reversible in eukaryotic messenger RNA. This modification is catalyzed by methyltransferases known as “writers” composed of a core complex consisting of Methyltransferase-like 3 (mettl3), mettl14, Wilms tumor 1-associating protein (wtap), and other subunits; can be reversed by “erasers” including fat mass and obesity-associated protein (fto) and AlkB homolog 5 (alkbh5); and can be recognized by “readers” such as YTH domain-containing proteins (ythdc1-2 and ythdf1-3) (Roundtree et al., 2017). It is an essential regulator of gene expression and can regulate the metabolism and processing of transcripts (Frye et al., 2018). The functions of m6A depend on its recognition by specific reader proteins that bind to methylated mRNAs. For example, in the nucleus, ythdc1 regulates splicing by binding to pre-RNAs, while in the cytoplasm, ythdf2 targets the transcripts for degradation, and ythdf1, ythdf3, and ythdc2 promote the translation of methylated mRNA (Patil et al., 2017; Liao et al., 2018). M6A plays critical roles in a variety of biological processes such as the maintenance of mRNA stability, miRNA maturation (Alarcón et al., 2015), stemness (Batista et al., 2014), circadian rhythms (Jean-Michel et al., 2013), and disease (Du et al., 2019). Increasing evidence has revealed crucial roles for m6A modification during sex determination, gonadal development, and gametogenesis in various taxa. In Drosophila, Ime4, the mettl3 homolog, can regulate sex by facilitating the alternative splicing of sex lethal (sxl), an important regulator of dosage compensation and the key sex determination factor (Haussmann et al., 2016; Lence et al., 2016). Other works have shown that m6A is involved in mammalian fertility and spermatogenesis (Hsu et al., 2017; Xu et al., 2017), amphibious oocyte meiotic maturation (Qi et al., 2016), and avian follicle selection (Fan et al., 2019). Deletion of mettl3, mettl14, fto, alkbh5, ythdf2, or ythdc2 leads to impaired fertility and defects in gametogenesis (Frye et al., 2018; Huang and Yin, 2018).

Fish are the most abundant vertebrates on earth, and their sex determination mechanisms exhibit extraordinary plasticity, diversity, and adaptability (Mei and Gui, 2015). To date, master sex-determination genes have been identified for only a few fish species, and fish sex determination and differentiation can be affected by the genetic system, exogenous hormones, and environmental factors (Capel, 2017). Epigenetic mechanisms provide fish with the ability to integrate genomic and environmental information for the development of a particular sex phenotype. Several steroidogenic enzymes and transcription factors are likely to be epigenetically regulated (Piferrer, 2018), including the important steroidogenic enzyme cytochrome P450 aromatase (encoded by cyp19a), but the understanding of the mechanism is limited to the DNA methylation of these genes’ promoters. As mentioned above, RNA modifications provide additional mechanisms of gene expression regulation without a DNA sequence change (Li et al., 2017) and could provide new insight into fish sex determination. Recently, m6A and RNA editing have been found to play important roles in the development of the gonads and other tissues. However, scant information is available for fish gonads barring one report addressing m6A in zebrafish (Danio rerio) showing that knockout of the m6A methyltransferase “writer” mettl3 led to a defect in sperm maturation and reduced sperm motility as well as a significant decrease in 11-ketotestosterone (11-KT) and 17β-estradiol (E2) levels (Xia et al., 2018). There has been no information about RNA editing in fish gonads reported to date. Therefore, these two RNA modifications remain largely unknown in fish gonads. The distribution patterns and functional relevance of these two modifications should be studied as an initial step to reveal the roles of m6A and RNA editing in gonads.

The olive flounder, Paralichthys olivaceus, is a commercially important marine fish in China, Japan, and Korea. It exhibits an XX (female)/XY (male) sex determination mechanism (Yamamoto, 1999), but its sex determination and gonadal development are also influenced by environmental factors, especially temperature (Sun et al., 2010). Previous studies have proven that the promoter DNA methylation levels of some steroidogenic-related genes are different between the ovary and testis (Wen et al., 2014; Si et al., 2016; Fan et al., 2017), which indicates that sex determination and gonadal development could be regulated by epigenetic mechanisms in this species. Similar to most fish, there is no available information on RNA modifications of the flounder gonad. Gynogenesis in fish can rapidly produce offspring with increased genetic similarity and is an efficient way of breeding (Yamamoto, 1999; You et al., 2001). The genotype of gynogenetic flounders is XX (female), but the phenotypic sex can be reversed to male through changes in environmental factors. Thus, this species could provide an appropriate model for studying sex determination and maintenance of phenotypic sex in fish. In this study, we systematically characterized the RNA editome based on high-throughput RNA sequencing data and whole-genome sequencing data for the gonads and muscle of the wildtype and meiogynogenetic diploid flounders. Additionally, RNA methylation profiling was performed via methylated RNA immunoprecipitation sequencing (MeRIP-Seq) to examine the differences in m6A modifications between the ovary and testis. This RNA modification profiling of the male and female flounder gonads adequately elucidates the sexual molecular mechanism and may provide new insights into the molecular mechanisms underlying the gene regulatory network controlling sex determination and subsequent maintenance of phenotypic sex in fish.

Results

Developmental Stage of the Flounder Gonadal Samples

According to the HE results from the gonadal histological sections (Figure 1), the developmental stage of all the ovaries used in this study was stage II, which included a large number of oocytes of phase II. In oocytes, one or numerous large nucleoli were at the periphery of the nucleus, and the nucleus was evident. Follicular cells appeared outside of the cell membrane of a small number of oocytes. The sizes of oocytes were bigger than those of oogonia, and cytoplasm were bluish violet because of basophilia. The developmental stage of the testes used was IV-V. The primary spermatocytes, secondary spermatocytes, and spermatids were all observed, which showed typical characteristics of stage IV of the testis. In the meanwhile, more mature spermatozoa were also found and spermatogonia were almost invisible in the sections showing the testis were developing into stage V.

Genome-Wide Analysis of RNA Editing Sites in the Flounder Gonads

Identification of RNA Editing Sites

Candidate RNA editing sites at the transcriptome-wide level were detected by RNA sequencing in the gonads and matched DNA sequencing in the muscle from the WM, WF, GM, and GF samples. A summary of the sequencing results is provided in Table 1. Approximately 96.76% of the 520.34 million pass-filter reads obtained from DNA sequencing were successfully aligned to the flounder reference genome. After quality trimming, 49.79 million RNA reads were generated from each sample on average, with an average mapping rate of 81.98%. Possible RNA editing events were detected by comparing matched RNA-seq and DNA-seq data. A site was considered to be a potential editing site when it was detected in at least two individuals. Combining these approaches, a total of 742 editing sites were identified from the RNA-seq data of the gonads (Supplementary Table S1).

Characterization of the Flounder Gonadal RNA Editing Sites

The distribution of the genomic locations of A-to-I, C-to-U, and non-canonical editing sites is shown in Figure 2A. Among all editing sites, a large fraction (364, 49.06%) were located in the non-coding 3′UTRs, and 61 (8.22%) sites were resided in coding sequences, 21 of which resulted in an amino acid change (Figure 2A). All intergenic-associated sites were near annotated genes, and some might represent extended 5′UTRs or 3′UTRs. Among the 742 detected editing sites, 551 editing events were located within annotated genes and were used for further analysis. All 12 possible base substitutions resulting from RNA editing except for A-to-C substitution were detected in this study (Figure 2B). A-to-G substitution was found to be the most common type, and there were a total of 459 A-to-G editing sites detected in the flounder gonads, accounting for up to 83.30% of the identified RNA editing sites. C-to-U substitutions accounted for approximately 3% of the sites (12 sites). Furthermore, most of the edited genes were edited in only one type of genic regions, such as an intron, CDS or 3′UTR.

Different RNA Editing Distributions in the Flounder Ovary and Testis

The RNA editing sites varied among the samples. Overall, number of the RNA editing sites in the testis was higher than that in the ovary. Average numbers of the editing events in the gonads of the WM, WF, GM, and GF samples were 610, 307, 614, and 349, respectively. The sample with the highest number of the edited sites was WM3 (n = 626, Figure 3A), while the sample with the lowest number was WF3 (n = 295, Figure 3A). However, the average RNA editing ratios were similar across different groups; in the WF, WM, GF, and GM samples, the ratios were 0.50, 0.47, 0.46, and 0.46, respectively (Figure 3B). Most of the mean editing ratios between the samples were not significantly different (Supplementary Table S2).

A-to-I editing catalyzed with ADAR enzymes was the most abundant editing type in the flounder gonads. The expression levels of the ADAR enzyme genes (adar1, adar2, and adar3) in the gonads were also analyzed to investigate whether the differences in the RNA editing pattern were related to the expression of the adar genes. The results determined by RNA-seq and qPCR showed that the adars presented higher expression in the testis than in the ovary, which was in accordance with the higher number of the RNA editing sites in the testis (Figure 3C). As shown in Figure 3C, there was a relationship between the adar expression and RNA editing patterns in different samples. Overall, we found a trend of a positive correlation (Pearson correlation = 0.948, p < 0.01) between the expression of the adar1 and tissue-specific RNA editing.

To explain the differential or sex-specific RNA editing sites in the flounder gonads, the expression of the genes edited in the ovary and testis was investigated. Among the 255 edited genes, 97 and 93 genes were differentially expressed in the WF vs. WM and GF vs. GM comparisons, respectively. Most of the genes expressed in the testis showed higher expression than those expressed in the ovary (Figure 3D). Within the detected editing sites, the editing ratios of 299 sites were different between the ovary and testis, and most editing events occurred only in the testis. Some of these sex-differential editing genes were related to gonadal development and gametogenesis, such as Cullin-3 (cul3), SRSF protein kinase 2 (srpk2), and DAZ-associated protein 2 (dazap2). Among the 255 genes, 73 were mapped to KEGG pathways. These genes were mainly associated with the “ErbB signaling pathway,” “MAPK signaling pathway,” and “Salmonella infection” (Figure 3E). KEGG enrichment revealed that DEGs with different RNA editing levels were related to metabolic pathways, cell adhesion molecules (CAMs), and microRNAs in cancer.

Transcriptome-Wide N6-methyladenosine Methylome Profiling of the Flounder Gonads

Transcriptome-Wide Detection of m6A Modifications

The MeRIP-Seq technique was applied to identify the target regions modified by m6A methylation in the flounder gonads. After filtering out low-quality data, 56,340,494, 59,429,986, 57,077,566, and 58,165,726 clean reads were obtained from the immunoprecipitated (IP) samples of the WM, WF, GM, and GF gonads, respectively. More than 50% of the IP reads uniquely mapped to the reference genome. Comparison of the m6A-IP data with the input data revealed 19,404 potential m6A peaks among 11,017 expressed genes in WF, 23,087 m6A peaks among 12,502 expressed genes in WM, 19,349 m6A peaks among 11,140 expressed genes in GF, and 22,084 m6A peaks among 12,176 expressed genes in GM. The information indicated that the flounder transcriptomes presented an estimated 1.793, 1.746, 1.702, and 1.637 m6A peaks per actively expressed transcript in WF, WM, GF, and GM, respectively.

To understand the topological pattern of m6A methylation in the flounder gonadal transcriptome, the distribution profiles of m6A peaks in the entire transcriptomes of these gonads were also investigated. It was shown that m6A peaks were abundant in exons and 3′UTRs in both the ovary and testis, which showed higher abundance of these peaks than 5′UTRs (Figure 4A). Motifs that were enriched in regions surrounding m6A peaks were analyzed, and the results showed that the most frequent motif was GGACU (Figure 4B) in all groups, which was consistent with the classic consensus sequence of RRACH.

Different m6A Methylomes in the Flounder Ovary and Testis

The m6A marks in the transcripts of the ovary and testis were compared, and 4,741 and 5,313 m6A sites in 3,947 and 4,468 genes were found to be differentially methylated between the gonads of WF vs. WM and GF vs. GM, respectively. The different peaks were mostly present in CDs and 3′UTRs of the mRNAs but were also found in 5′UTRs (Figure 4C). Furthermore, the transcript expression levels of genes with differentially methylated m6A sites were analyzed. The association analysis between them was performed based on the RNA-Seq and MeIP-seq sequencing data. The results showed that among the 3,947 and 4,468 differentially methylated genes, the expression levels of 1,717 and 1,796 genes were significantly changed (Figure 4D). Among these DEGs that also exhibited differential m6A RNA methylation levels, 58% and 60% of the DEGs appeared to be negatively correlated with RNA methylation levels in the WF vs. WM and GF vs. GM comparisons, respectively.

GO and KEGG enrichment analyses were conducted to decipher the regulatory roles of the differential RNA methylation of genes with different expression levels between the ovary and testis. The top results of the GO enrichment analysis are shown in a bar graph in Figure 5A and Supplementary Table S3, and the gene functional enrichment analysis showed that these genes were significantly enriched in cellular functional categories relevant to gonadal development and gametogenesis, such as “ciliary basal body,” “cilium,” “sequence-specific DNA binding,” and “Cul3-RING ubiquitin ligase complex.” The results of the KEGG pathway analysis of differentially methylated transcripts between the ovary and testis showed that the genes identified in the gonads of the WF vs. WM and GF vs. GM comparisons mapped to 272 and 260 pathways, respectively (Figure 5B and Supplementary Table S4). Among these pathways, several well-documented pathways associated with gonadal development and maintenance were identified as m6A targets controlled by m6A modifiers, including the “TGF-beta signaling pathway,” “Wnt signaling pathway,” “ovarian steroidogenesis,” “estrogen signaling pathway,” and “AMPK signaling pathway.” Notably, the FoxO signaling pathway was one of the most enriched pathways identified in both the WF vs. WM and GF vs. GM comparisons, implying a potential important role in gonadal development and gametogenesis. Within this pathway, 48 genes exhibited both differential m6A RNA methylation and gene expression levels (Figure 5C). In addition, some pathways that have not been previously investigated in depth during gonadal development in fish, such as the PPAR and RNA degradation pathways, were identified. PPAR showed high dissimilarity in both RNA methylation and gene expression between the testis and ovary. The significant changes in this pathway were driven by alterations in the expression levels of 12 related genes.

To identify differentially expressed transcripts with differential RNA methylation in the male and female gonads, the ovaries or testes of the wild and gynogenetic fish were grouped together as the female and male gonads to be compared. In total, 728 genes showed upregulation of both RNA methylation and expression, while 960 genes showed upregulation of RNA methylation and downregulation of expression, 673 genes showed downregulation of RNA methylation and upregulation of expression, and 618 genes showed downregulation of both methylation and expression. The RNA methylation levels and expression patterns of sex-related genes in fish or other vertebrates were checked. These genes have been reported as playing important roles in sex determination/differentiation, sex maintenance or gametogenesis. In total, 36 sex-related candidate genes were selected (Table 2), among which 18 were significantly more highly expressed in the testis compared to the ovary, while the others were more highly expressed in the ovary. Genes that are essential for testis differentiation and maintenance of male-specified germ cells in fishes and other vertebrates, such as doublesex and mab-3-related transcription factor 1 (dmrt1), sry box-containing gene 8 (sox8), sry box-containing gene 9 (sox9), anti-Müllerian hormone (amh), Wilms tumor 1 (wt1), deleted in azoospermia like (dazl), growth arrest and DNA damage-inducible gamma (gadd45g), and SIX homeobox 4 (six4), were identified and shown a male-biased expression pattern in the flounder. Sperm-associated antigen 5 (spag5), sperm-associated antigen 8 (spag8), sperm-associated antigen 17 (spag17), and spermatogenesis-associated protein 7 (spata7) were significantly more highly expressed in the testis compared to the ovary. Among the female-biased genes, some genes that have been proven to be involved in the regulation of ovarian differentiation, including bone morphogenetic protein 15 (bmp15), catenin beta 1 (ctnnb1), foxo3, and LIM homeobox 9 (lhx9), and genes encoding proteins implicated in oogenesis, such as zona pellucida sperm-binding protein 3 (zp3), were significantly highly expressed in the ovary.

The genes encoding the methyltransferase ‘writer’ core complex (mettl3, mettl14, wtap, and other subunits, including Virilizer (kiaa1429), RNA-binding Motif Protein 15 (rbm15), and E3 ubiquitin ligase Casitas B-lineage lymphoma-transforming sequence-like protein 1 (cbll1, also known as HAKAI), were identified in the flounder transcriptome. In addition, the m6A “eraser” genes fto and alkbh5, the YTH domain-containing protein “reader” genes ythdf1, ythdf2, ythdf3, ythdc1, and ythdc2, and genes encoding other families of m6A readers, including eukaryotic initiation factor 3 (eif3), insulin-like growth factor 2 mRNA-binding protein 2 (igf2bp2), fragile X mental retardation 1 (fmr1), and Proline-rich coiled-coil 2A (prrc2a), were also detected in this study. The expression of mettl3, mettl14, wtap, fto, alkbh5, ythdf1, ythdf2, ythdf3, ythdc1, and ythdc2 was investigated in the flounder ovary and testis with qPCR. The detected expression patterns of these genes were similar to those obtained by RNA-seq, indicating the reliability of the transcriptome expression analysis. The results also showed that the m6A “writer” mettl14 was more highly abundant in the ovary than in the testis, while the “erasers” fto and alkbh5 (p < 0.05) and the “readers” ythdf1 and ythdc2 exhibited especially high levels in the testis (p < 0.05, Figure 6).

Comparison of Genes and Pathways Modified by A-to-I Editing and m6A Methylation

A total of 125 genes were modified by both m6A and A-to-I editing, and these genes were significantly enriched in signaling pathways regulating the pluripotency of stem cells, the Hippo signaling pathway, and the pancreatic cancer pathway. Thirty of these genes showed significantly different expression levels between the ovary and testis (Table 3). Both RNA modifications occurred at adenosines, but the RNA editing sites were mainly located in 3′UTRs, which did not coincide with the positions of the m6A modification. By comparing the pathways modified by A-to-I editing and m6A methylation, we found that the Wnt and MAPK signaling pathways were affected by both m6A and RNA editing modifications.

Intriguingly, several ubiquitin-related genes, such as the critical gene of the Cullin-3-based E3 ubiquitin ligase complex, cullin3, were regulated by RNA editing, while some GO terms related to ubiquitin, “cilium,” and the “Cul3-RING ubiquitin ligase complex” were enriched in the DEGs with different m6A modifications. The E3 ubiquitin ligase HAKAI, one of the core constituents of the plant m6A writer complex, was also identified and shown female-biased expression.

Comparison of Gene Expression Profiles in the Flounder Ovary and Testis

Multiple comparisons were conducted between the gonads of the WF vs. WM and GF vs. GM groups to evaluate the differential expression of genes between the ovary and testis. Overall, a total of 23,127 genes were identified from the 12 flounder gonadal transcriptomes in this study, and a large number of them showed sex-biased expression (fold change > 2 and adjusted p-value ≤ 0.05). The DEG numbers in different groups are shown in Figure 7A. There were 12,402 and 11,897 genes showing different expression levels between the WF vs. WM and GF vs. GM, respectively. Compared with WF, 6,838 upregulated and 5,564 downregulated unigenes were identified in WM, while 6,446 upregulated and 5,451 downregulated unigenes were identified in GM compared to GF (Figure 7A). Principal component analysis (PCA) plot and heatmap both showed that differences in gene expression were much more pronounced between the ovary and testis than between the wildtype and gynogenetic gonads (Figures 7B,C). Compared with the wildtype ovary, 300 upregulated and 927 downregulated unigenes were identified in the gynogenetic ovary, while compared with the wildtype testis, 288 upregulated and 212 downregulated unigenes were identified in the gynogenetic testis.

As expected, cyp19a1a [aromatase converting testosterone (T) into E2] expression was higher in the ovary. However, 11β-hydroxylase (cyp11b) and 11β-hydroxysteroid dehydrogenase 2 (hsd11b2), which convert T into 11-KT, were not detected in either the ovary or testis. In addition, some genes involved in sex steroid synthesis mostly showed higher expression levels in the testis than in the ovary, such as steroidogenic acute regulatory protein gene (star1), cytochrome P450c17 (cyp17a1), and 3β-hydroxy-Δ5-C27-steroid oxidoreductase (hsd3b7), while the 17β-hydroxysteroid dehydrogenase 1 gene (hsd17b1) showed female-biased expression.

The WF vs. WM and the GF vs. GM comparisons revealed 14 and 13 significantly different KEGG pathways, respectively (Figure 7D and Supplementary Table S5). The results showed that DEGs identified in the WF vs. WM comparison were related to the pathways of “Cellular senescence,” “Phenylalanine metabolism,” “Glycosylphosphatidylinositol (GPI)-anchor biosynthesis,” “Small cell lung cancer,” “Sphingolipid metabolism,” “Arachidonic acid metabolism,” “Ribosome biogenesis in eukaryotes,” “Terpenoid backbone biosynthesis,” “GnRH signaling pathway,” “Fatty acid biosynthesis,” “Protein processing in endoplasmic reticulum,” “Glycosaminoglycan biosynthesis – chondroitin sulfate/dermatan sulfate,” “Antifolate resistance,” and “Thyroid hormone signaling pathway,” while the DEGs identified in the GF vs. GM comparison were related to “Ribosome biogenesis in eukaryotes,” “Phenylalanine metabolism,” “Thyroid cancer,” “Primary bile acid biosynthesis,” “Sphingolipid metabolism,” “Arginine and proline metabolism,” “Terpenoid backbone biosynthesis,” “Arachidonic acid metabolism,” “Glycerophospholipid metabolism,” “Glycerolipid metabolism,” “Limonene and pinene degradation,” “Other types of O-glycan biosynthesis,” and the “PPAR signaling pathway.”

Discussion

RNA Editing in the Flounder Ovary and Testis

RNA editing can increase transcriptome diversity and flexibility without mutations occurring at the DNA level (Hwang et al., 2016). It can affect the structure and coding of RNA and generate proteomic diversity (Koslowsky, 2004). After the identification of the RNA editing in trypanosomes in 1986, the RNA editome was well characterized in mammals such as humans (Peng et al., 2012), mice (Danecek et al., 2012), bovines (Bakhtiarizadeh et al., 2018), and porcines (Zhang et al., 2019). However, no related reports have been published for fish species with the exception of the report by Guryev et al. (2006) describing the presence of C-to-U and A-to-I RNA editing events in zebrafish. In this study, we systematically detected and characterized the RNA editome in the gonads of the flounder based on RNA and corresponding DNA sequencing data. The results showed that the canonical A-to-I RNA editing event was the most common type of editing in the flounder, which was consistent with previous studies in mammalian species (Hwang et al., 2016; Bakhtiarizadeh et al., 2018). However, most of the editing sites in the flounder gonads targeted non-coding regions (3′UTRs), which differs from findings in humans (intronic) (Nishikura, 2016) and porcines (intergenic) (Zhang et al., 2019). Compared to mammals, the ratio of editing sites in coding regions is higher in the flounder. The large number of editing sites occurring in introns and 3′UTRs of the flounder gonads indicated that RNA editing in fish may play fundamental roles in the regulation of splicing and miRNA regulation, respectively (Hwang et al., 2016).

Notably, our results also showed that the number of editing sites was different between the ovary and testis. Compared to the testis, the lower number of editing sites identified in the ovary suggested lower RNA editing activity. Most of the genes harboring RNA editing sites were more highly expressed in the testis than in the ovary. Thus, RNA editing may regulate gene expression, but the regulatory mechanism will be subject to further study. Annotation of the edited genes in the flounder gonads revealed that some of the genes were associated with gonadal development and gametogenesis. For instance, cullin3, a component of the Cullin-3-based E3 ubiquitin ligase complex, is necessary for caspase activation. Mutations in this gene reduce the effect of caspase activation and block apoptosis during spermatid individualization (Richburg et al., 2014). Srpk2 is highly expressed in mouse testis and might participate in precursor mRNA splicing during mouse spermiogenesis (Yang et al., 2015). Dazap2 interacts with the germ-cell-specific RNA-binding proteins Daz and Dazl1 (Yen, 2004).

Similar to what is observed in mammals, three independent genes of the ADAR family were identified in the flounder transcriptome: adar (adar1), adarb1 (adar2), and adarb2 (adar3). According to studies in mice, ADAR enzymes are essential for normal life and development, and knockout or homozygous null mutations in ADAR genes lead to lethality (Higuchi et al., 2000; Wang et al., 2000). In the flounder, adar1 showed the highest expression level among these three genes, while adar3 showed a very low expression level in the gonads. Only adar1 and adar2 have been shown to be enzymatically active (Nishikura, 2016), while adar3 is exclusively expressed in the brain and can inhibit RNA editing by competitively binding double-stranded RNA (Picardi et al., 2015; Oakes et al., 2017). Our results also proved that the distribution of A-to-G sites was significantly and positively correlated with the expression of adar1, suggesting that editing enzyme expression may play an essential role in regulating tissue-specific editing levels.

M6A Profiles of the Flounder Gonads

M6A represents the most prevalent epigenetic modification of RNAs. Many m6A profiling studies have been performed in mammals, avians, and plants, but no related reports have been published for fish except for zebrafish (Zhang et al., 2017). A study in zebrafish revealed conserved features of the m6A methylome and preferential distribution of m6A peaks near stop codons with a consensus RRACH motif. By using the MeRIP-seq method, we produced the first transcriptome-wide m6A modification profile and detected m6A methylation sites in the flounder gonadal transcriptome. There were approximately 1.7 m6A peaks per transcript, which was similar to the situation in zebrafish (1.7) (Zhang et al., 2017) but slightly higher than the number found in mammals (1.5) (Dominissini et al., 2012). M6A modification can regulate gene function by controlling RNA splicing, nuclear export, and RNA degradation and translation (Frye et al., 2018), and is profoundly implicated in biological processes in animal gonads. Studies show that deficiency of the alkbh5 demethylase leads to aberrant spermatogenesis and apoptosis with impaired fertility in the testis (Zheng et al., 2013), while ythdf2 is required for the generation of the maternal transcriptome during oocyte maturation (Ivanova et al., 2017), and ythdc1 is essential for the production of spermatogonia in males and oocyte maturation in females (Kasowitz et al., 2018). Furthermore, the simultaneous deletion of mettl14 and mettl3 in advanced mouse germ cells results in a decrease in sperm motility and an increase in the abnormal sperm ratio (Lin et al., 2017). In fish, the understanding of the function of m6A is still limited. It has been proven that zebrafish mettl3 is essential for sperm maturation and motility and can regulate 11-KT and E2 levels (Xia et al., 2018). In our study, the expression of the “writers” mettl14 and cbll1 was found to be more highly abundant in the ovary than in the testis, while the “erasers” fto and alkbh5 and the “readers” ythdf1 and ythdc2, which increase the translation of methylated mRNA, were expressed at especially high levels in the testis.

Several genes that have been reported as playing important roles or participating in sex determination/differentiation or gametogenesis showed differential m6A RNA methylation in the flounder ovary and testis. For example, dmrt1, amh, dazl, and Sox protein genes (sox4, sox8, sox9, and sox11) are important for the regulation of testis development and spermatogenesis (Barrionuevo and Scherer, 2010; Pfennig et al., 2015; Webster et al., 2017). Spag5, -8, -17, and spata7 were significantly more highly expressed in the testis than in the ovary, and the expression difference might be associated with spermatogenesis (Yang et al., 2016). The zona pellucida (ZP) is an extracellular glycoproteinaceous matrix shell that is involved in oocyte and gamete development, and the expression levels of zp gene mRNAs were found to be significantly increased during oogenesis in fish, especially at the previtellogenic stage, when they showed a higher expression level than that at the undeveloped stage (Zeng and Gong, 2002). Zp3 showed higher expression levels in the flounder ovary compared to the testis, and previous studies showed that this gene represents a major class of female-specific molecule with a role in reproduction (Liu et al., 2006). Although differential m6A modification of these genes between the ovary and testis was detected, the functions of these genes and their regulation by m6A modification are still unknown and should be further studied.

The DEGs with differential m6A RNA methylation were enriched in several canonical signaling pathways. These included pathways involved in gonadal differentiation and development in fish and other vertebrates, such as the TGF-beta signaling pathway, Wnt signaling pathway, ovarian steroidogenesis pathway, and FoxO signaling pathway. Among these pathways, the FoxO signaling pathway was enriched in both the WF vs. WM and GF vs. GM comparisons. FoxO proteins of different subclasses are transcription factors that play important roles in cell cycle arrest, apoptosis, and stress responses in vitro (Taisuke et al., 2004). Biochemical and genetic studies in mice have implicated the Pten/Pi3k/Akt/FoxO3 pathway as a major signaling pathway involved in the regulation of dormancy and initial follicular activation in the ovary (Lee and Chang, 2019). In the current study, the mRNAs encoding pten, pi3k, and foxo3 were all modified by m6A methylation. Intriguingly, they showed opposite m6A distribution patterns between the ovary and testis. Our results also included pathways whose roles have not been previously investigated in depth in relation to fish sex differentiation and gonadal development, such as the PPAR and RNA degradation pathways. A previous study in mammals indicated that PPAR is capable of regulating the function of the reproductive organs by modulating the expression of related steroidogenic enzymes through Star (Kowalewski et al., 2009) and PPAR-γ can regulate progesterone production in ovarian granulosa cells by coactivation with sf1 and lhxr1 (Yazawa et al., 2010). PPAR and the key signaling molecule RXRG in this pathway can promote the differentiation of chicken embryonic stem cells into primordial germ cells (Cheng et al., 2018). In this study, the m6A RNA methylation levels of several critical genes in this pathway, including sirt1 and RXR, were found to be different in the ovary and testis, but its function remains unclear in the gonads. The GO analysis showed that the cilia-related genes were significantly enriched. Primary cilium was identified as a key coordinator of signaling during organogenesis, and it had been shown that the primary cilium participated in receiving signals from Notch, PDGFRα, FGF, Wnt, and Hedgehog pathways (Wainwright et al., 2014). Primary cilia-related genes are present in different cell types of gonads and may be important for the differentiation of gonad cells and sex differentiation of gonads in mammals (Piprek et al., 2019). However, the role of primary cilium in the sex determination and gonadal differentiation in fish is completely unknown. The possible direct and indirect effects mediated by m6A modification remain to be further elucidated.

The Relationship Between m6A and A-to-I Editing in the Flounder Gonads

Among the different types of RNA modifications, m6A and A-to-I editing are two of the most abundant. Although both modifications occur on A bases, their catalytic mechanisms are distinct. A-to-I conversion is catalyzed by ADARs that preferentially bind to double-stranded RNA substrates (Nishikura, 2016), which are dependent on the formation of RNA secondary structures (Bahn et al., 2015). In contrast, m6A RNA methylomes exhibit enrichment of the RRACH motif at m6A sites. Thus, the different sequences and structural features of the A-to-I and m6A modifications suggest that these two chemical modifications are unlikely to modify the same A bases. Xiang et al. (2018) demonstrated a global A-to-I difference between m6A-positive and m6A-negative RNA populations in mammals and indicated that there was a negative correlation between m6A and A-to-I modifications (Xiang et al., 2018). In this study, 125 genes were found to be modified by both the m6A and A-to-I editing, 30 of which showed significantly different expression levels between the ovary and testis. Most of the positions of the m6A modifications and RNA editing were not coincident, which was in accordance with previous research results (Xiang et al., 2018). However, the regulatory relationship between the A-to-I and m6A modifications is unclear in the flounder and requires further study.

The teleost gonadal development is a dynamic process with various phases that are similar in different species. Testicular stage is generally based on the relative proportions of spermatocytes, spermatids, and spermatozoa (Blazer, 2002). Histological analysis suggested the testes used in this study were at stage IV to V. At stage IV, the primary spermatocytes, secondary spermatocytes, and spermatids were arranged in tubule seminiferous. Some of the spermatids developed into spermatozoa. While at stage V, the mature spermatozoa predominated in the testis and the number of spermatocytes was obviously less than that at stage IV. Thus, the sampled testes were at spermiogenesis stage with an increasing proportion of spermatozoa. Numerous well documented sex-related genes involved in gonadal development and gametogenesis were differently modified between the ovary and testis in this study. Among the candidate genes, a subset of well-established genes involved in spermatogenesis (e.g., dazl, dazap2, amh, gadd45g, srpk2, Spag5, -8, -17, and spata7) were detected. For instance, dazl is a master translational regulator essential for spermatogenesis (Li et al., 2019). Amh is expressed in the Sertoli cells that surrounded the germ cells and involved in the maintenance of the flounder testis (Wang et al., 2020). Gadd45g is necessary for activation of the male sexual pathway in mice and its absence leads to sex reversal (Johnen et al., 2013). In contrast, the developmental stage of the ovary is classified based on ogenesis process, such as oocyte size, presence of follicular layer, and number of nucleoli. The ovaries in this study were at stage II with a large number of phase II oocytes which were bigger than oogonia. Correspondingly, several genes related to ovarian or oocyte development (e.g., bmp15, foxo3, and nr0b1) were identified. The bmp15, an oocyte-expressed signaling molecule, is required for maintenance of the female sexual fate in zebrafish (Dranow et al., 2016). The foxo3 is a critical transcriptional regulator of fish ovarian development by regulating the ovarian germ cells and follicular cells (Liu et al., 2016). Nr0b1, also called dax1, mainly locates in fish primary oocytes and previtellogenic oocyte cells, and its mutation in zebrafish causes female-to-male sex reversal (Chen et al., 2016). In this study, we obtained a global view of the RNA modification differences between male and female fish gonads. Nevertheless, the gonadal development is a dynamic process with different gene expression patterns at different developmental stages. As specific developmental stage samples were used, this study could not elucidate the dynamics of gene expression and RNA modifications over the course of gonadal development. Thus, the dynamic changes of the RNA modifications and gene expression at all gonadal development stages should be analyzed in the near future. Moreover, it should be noted that the RNA modification patterns of gonads may differ depending on the cell types. As the entire gonads were used in the present study, the RNA modifications in this study represent the combination of the different cell types. Further study focusing on single-cell analysis would increase our knowledge of gonadal development and sexual plasticity in fish.

RNA editing and m6A have been recognized as targets for drug discovery, and some drugs for inhibiting RNA editing enzymes or m6A modification regulators were discovered recently (Christofi and Zaravinos, 2019; Ianniello et al., 2019). Our study showed the levels of the RNA editing and m6A were different in the flounder ovary and testis, and some sex-related genes (e.g., dmrt1, amh, and wt1) may be regulated by RNA modifications. Therefore, these RNA modification regulators and sex-related genes provide new potential targets for therapeutic products in fish sex control. However, the mechanisms of sex determination in fish are extremely complex, and identifying and understanding biological functions of RNA modification in fish sex determination and differentiation are important for designing reasonable chemical/drug.

Despite the substantial results have obtained in this study, some restrictions due to sequencing limitations remained. MeRIP-seq is the most common method to elucidate global mRNA m6A sites, but this approach could not accurately reflect the positions of m6A residues (Liu et al., 2013). The methylated regions are detected as peaks of approximately 100–200 base pairs that contain multiple DRACH motifs in transcript relative to input RNA. Moreover, antibodies for m6A can also detect a second base modification, N6,2′-O-dimethyladenosine (m6Am), which has been found at a lower abundance than m6A (Mauer and Jaffrey, 2018). Another common method for transcriptome-wide studies of m6A was the m6A individual-nucleotide-resolution cross-linking and immunoprecipitation (miCLIP), which is not influenced by peak shapes and not restricted to DRACH motifs (Linder et al., 2015). However, MeRIP-seq is still more often used than miCLIP, as it follows a simpler protocol and requires less starting RNA (McIntyre et al., 2020). Next-generation sequencing (NGS) technologies enable genome-wide identification of RNA-editing sites, and RNA editing studies are widely implemented by comparing matched high-throughput RNA sequencing and DNA sequencing data currently (Wang et al., 2016; Zhang et al., 2019). Currently, many bioinformatics tools can identify RNA-editing sites from RNA-seq data, but the accurate identification of the RNA editing sites may suffer from methodical artifacts such as reverse transcription errors, sequencing errors, alignment errors, and single nucleotide polymorphisms (SNP) (Wang et al., 2016; Guo et al., 2017). Thus, individual-nucleotide-resolution of m6A modification in the flounder gonads and the verification of identified RNA editing sites should be carried out in the subsequent study.

Gynogenetic Flounder Is an Appropriate Model for Studying Epigenetic Factors Involved in Fish Sex Determination and Differentiation

Gynogenesis is a form of genetically all-female reproduction in which the male pronucleus does not fuse with the female pronucleus in gynogenetic zygotes (Neaves and Baumann, 2011). The flounder exhibits an XX (female)/XY (male) sex determination system, and all gynogenetic female flounders should exhibit an XX genotype. However, the sex determination and gonadal development of the flounder are determined by genetic factors but are also affected by environmental factors, and even gynogenetic individuals can produce a male phenotype. Thus, the gynogenetic flounder possess same genome is able to produce two distinct sexual phenotypes in response to environmental cues with mechanisms not involving a genomic change. Epigenetic modifications provide an interface to integrate genetic factors and environmental signals during sex determination. Epigenetic factors such as DNA modifications have been reported to participate in the flounder sex differentiation. Our previous studies showed that the cyp19a promoter exhibited dimorphic methylation levels between the flounder testis and ovary, which were inversely correlated with its expression levels (Wen et al., 2014; Fan et al., 2017). To date, the epigenetic studies in fish sex determination mainly focus on heritable modifications of DNA, histones, and chromatin structure (Ortega-Recalde et al., 2020). RNA modifications such as RNA editing and m6A RNA methylation may shed new light on fish sex determination and development. However, as epitranscriptomic biomarkers, the roles of RNA editing and m6A modifications in the regulation of sex differentiation remain largely unknown. XX female and male fish could provide an appropriate model for studying this topic. Our previous study demonstrated that gene expression in gynogenetic diploids was similar to that in common diploids (Fan et al., 2016), which is consistent with the results of this study. The present study also showed that the identified RNA modifications were similar in the wildtype and gynogenetic gonads. For instance, the RNA editing sites and levels were almost the same, and there was no significant difference in the position or level of the m6A peaks between the wildtype and gynogenetic gonads.

Up till now, the studies of RNA modifications in sex determination were carried out in only a few model animals, such as Drosophila and zebrafish (Haussmann et al., 2016; Lence et al., 2016; Xia et al., 2018). In Drosophila, m6A can regulate sex by facilitating the alternative splicing of the key sex determination factor sxl. In zebrafish, knockout of the m6A methyltransferase mettl3 led to a defect in fertility as well as sex differentiation. These studies all focused on a specific modification regulator, and used the mutants to study its functions. The relatively long generation interval of the flounder causes a great deal of discomfort on such function researches. Nevertheless, compared to the zebrafish with a complex polygenic sex determining system, the gynogenetic flounder exhibits an XX female genotype with different phenotypic sexes. Sex differentiation of the flounder may be regulated by mechanisms not only involving genomic changes, but also involving epigenetic factors. Thus, the gynogenetic flounder probably provides a better model for addressing the epigenetic modification on fish sex determination issues. Moreover, the flounder is an aquaculture fish that is yet to reach industrial scale production, and elucidation of it sex differentiation is also a key area of applied research.

In conclusion, we profiled the RNA editing and m6A methylomes of the flounder ovary and testis. The results indicated that most A-to-I sites were located in the 3′UTRs of the mRNAs, and the number of the editing sites in the testis was higher than that in the ovary. The flounder gonadal transcriptome was extensively methylated with the m6A modifications, and the different peaks were mostly located in CDs and 3′UTRs of the mRNAs. The RNA editing and m6A modifications of RNA may be involved in the regulation of gonadal development and gametogenesis in the flounder gonads. The high-throughput epitranscriptome data in this study provide a foundation for understanding these critical RNA modifications in fish and expand our knowledge of the function of RNA modifications. Further researches using gonads at all developmental stages and during gonadal differentiation period, even using single cell type of gonads, are needed to shed light on the roles of RNA editing and m6A modifications in sexual plasticity and gonadal development of fish.

Materials and Methods

Fish and Sample Collection

The meiogynogenetic flounders were obtained through artificial induction as previously reported (You et al., 2001). Briefly, sperm were irradiated under ultraviolet light and then used to be fertilized with eggs. Five minutes after fertilization, the fertilized eggs were treated under cold shock by being immersed in 0°C seawater for 45 min, after which the eggs were incubated in normal seawater, and fertilized eggs that were not subjected cold shock were used as the haploid control group to confirm meiogynognesis induction. The fertilized eggs obtained with non-irradiated sperm and ordinary eggs were reared under same conditions as the control individuals (wildtype). The sexes of the meiogynogenetic and control fishes were identified by morphological observation of the gonads at approximately 1.5 years of age, and the developmental stages of the gonads were detected by using the histological method described as Sun et al. (2009). Gonadal and muscle tissues of three males (wildtype males, WM, 29.82 ± 0.86 cm total length, TL; 322.96 ± 47.36 g body weight, BW) and three females (wildtype females, WF, 31.95 ± 1.51 cm TL; 399.17 ± 50.84 g BW) from the flounders in the control group, and three males (gynogenetic males, GM, 29.71 ± 0.30 cm TL; 337.83 ± 37.32 g BW) and three females (gynogenetic females, GF, 31.84 ± 0.71 cm TL; 417.20 ± 31.63 g BW) from the flounders in the induced meiogynogenetic diploid group were collected. Gonadal and muscle samples were retrieved from the fish after they were anaesthetized with MS-222 (Sigma, United States). Each gonadal sample was divided into two halves. One half was fixed in Davison’s fixative solution for histological analysis, and the other half was immediately stored in liquid nitrogen for RNA isolation. The muscle samples were frozen immediately and then stored at −80°C for DNA extraction.

Nucleic Acid Isolation

Genomic DNA was extracted from the muscle samples by using the standard phenol-chloroform protocol. Only DNA samples with OD 260/280 ratios of 1.8–2.0 and total contents greater than 1.5 μg were used in the subsequent steps. Total RNA was isolated from the gonads using TRIzol reagent (Invitrogen, United States) according to the manufacturer’s instructions. The quality and quantity of the RNA were analyzed with Nanodrop 2000 (Thermo Fisher Scientific, United States) and Agilent 2100 Bioanalyzer (Agilent Technologies, United States) instruments, respectively. Testis RNA samples with RNA integrity number (RIN) scores higher than seven were used in this study. As in other fishes, the flounder ovary samples showed a strong peak of low-molecular-weight RNA, which masked the 18S and 28S rRNA peaks used for calculating the RIN. Therefore, the RIN score could not be used to measure RNA integrity in the ovary.

Transcriptome Sequencing

Three micrograms of total RNA per gonad were used for the subsequent generation of RNA-seq libraries. In total, 12 sequencing libraries were generated using mRNA purified from total RNA using TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, United States) following the manufacturer’s protocol. Each library preparation was sequenced on an Illumina HiSeq platform (Shanghai Personal Biotechnology Co., Ltd., Shanghai, China), and 150 bp paired-end reads were generated. The sequencing data were filtered with Cutadapt v2.7¹ (Martin, 2011). RNA-seq reads were mapped to the olive flounder genome using HISAT2² (Kim et al., 2015) with the default mismatch being no more than 2. HTSeq³ (Anders et al., 2015) was used to compare the Read Count values on each gene as the original expression of the gene, and then reads per kilobase per million (FPKM) was used to standardize the expression. DESeq v1.39.0 (Anders and Huber, 2010) was used to analyze the genes of different expression. Differentially expressed genes (DEGs) were selected on the basis of a log2 (fold change) > 1 or log2 (fold change) < −1 and a false discovery rate (FDR) < 0.05. Bi-directional clustering analysis of all different genes of samples was performed by R language Pheatmap v1.0.12 software package. The topGO v2.40.0 was used to map all the genes to terms in the Gene Ontology database and calculate the numbers of differentially enriched genes in each term.

Whole-Genome Sequencing

Muscle DNA from the 12 fishes was used to construct 12 DNA libraries. These libraries were generated using a TruSeq Nano DNA HT Sample Preparation Kit (Illumina, San Diego, CA, United States) according to the manufacturer’s recommendations. Then, each library was sequenced on an Illumina HiSeq platform (Shanghai Personal Biotechnology Co., Ltd., Shanghai, China), and 150 bp paired-end reads were generated for further analysis (10 × depth). The quality control of the raw reads was conducted with FastQC software v0.11.2⁴. The raw data were trimmed and filtered with AdapterRemoval and trimmomatic v. 0.36 (Bolger et al., 2014) and then were aligned to the reference genome. The reference flounder genome was downloaded from Ensembl (GCF_001970005.1_Flounder_ref_guided_V1.0).

RNA Editing Detection

RNA editing sites were detected by direct comparison between the RNA-seq data from the flounder gonads and the corresponding genomic data using RES-Scanner. All thresholds used for the identification of RNA editing sites were the default parameters of RES-Scanner (Wang et al., 2016). The low frequency variant detection parameters were set to 5%, and the minimum coverage, minimum count, and minimum frequency were set to 10, 2, and 5%, respectively. A site was considered to be a potential editing site when it was detected in at least two individuals. The statistical significance of the differences between the samples in terms of the editing ratio was assessed with Student’s unpaired t-test was implemented in SPSS16.0.

High-Throughput m6A and Input RNA Sequencing

Poly (A) RNA was purified from 200 μg total RNA using Dynabeads Oligo (dT) 25-61005 (Thermo Fisher, CA, United States) with two rounds of purification. Following purification, the poly (A) RNA was fragmented into small pieces using Magnesium RNA Fragmentation Module (NEB, cat.e6150, United States) for 7 min under 86°C. Then the cleaved RNA fragments were incubated for 2 h at 4°C with m6A-specific antibody (No. 202003, Synaptic Systems, Germany) in IP buffer (50 mM Tris-HCl, 750 mM NaCl, and 0.5% Igepal CA-630). The IP RNA was reverse-transcribed to create the cDNA by SuperScript^TM II Reverse Transcriptase (Invitrogen, cat. 1896649, United States), which were next used to synthesis U-labeled second-stranded DNAs with E. coli DNA polymerase I (NEB, cat.m0209, United States), RNase H (NEB, cat.m0297, United States), and dUTP Solution (Thermo Fisher, cat.R0133, United States). An A-base was then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contained a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single- or dual-index adapters were ligated to the fragments, and size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme (NEB, cat.m0280, United States) treatment of the U-labeled second-stranded DNAs, the ligated products were amplified with PCR with the following conditions: initial denaturation at 95°C for 3 min; 8 cycles of denaturation at 98°C for 15 s, annealing at 60°C for 15 s, and extension at 72°C for 30 s; and then final extension at 72°C for 5 min. The average insert size for the paired-end libraries was ∼100 ± 50 bp. At last, we performed the 2 × 150 bp paired-end sequencing (PE150) on an illumina Novaseq^TM6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) following the vendor’s recommended protocol.

The fastp software⁵ was used to remove the reads that contained adaptor contamination, low quality bases, and undetermined bases with default parameter. Then sequence quality of IP and Input samples was also verified using fastp. We used HISAT2 to map reads to the genome of the flounder with default parameters. The mapped reads of IP and input libraries were analyzed with R package exomePeak⁶ (Meng et al., 2014), which identified m6A peaks with bed or bigwig format that were adapted for visualization on the IGV software⁷. MEME⁸ (Bailey et al., 2009) and HOMER⁹ (Heinz et al., 2010) were used for de novo motif identification, followed by localization of the motif with respect to peak summit. Called peaks were annotated by intersection with gene architecture using R package ChIPseeker¹⁰ (Yu et al., 2015). Then StringTie¹¹ (Pertea et al., 2015) was used to perform expression level for all mRNAs from input libraries by calculating FPKM. The differentially expressed mRNAs were selected with log2 (fold change) > 1 or log2 (fold change) < −1 and p-value < 0.05 by R package edgeR¹² (Robinson et al., 2010).

Validation of Transcriptomic Data

Real-time quantitative polymerase chain reaction (qPCR) was performed to verify the expression profiles obtained from the transcriptomic data. The RNA used for qPCR was the same RNA used for RNA-seq analysis. Ten m6A methyltransferase genes (mettl3, mettl14, wtap, fto, alkbh5, ythdc1, ythdc2, ythdf1, ythdf2, and ythdc3) and three adar genes (adar1, adar2, and adar3) were chosen for analysis. The primers used in this study were designed based on sequences from the transcriptome using Primer Premier 6.0 (PREMIER Biosoft International, United States) and are listed in Supplementary Table S6. Three replicates were performed for each gene in each individual. The relative quantification result for each target gene was expressed as the fold variation over the reference genes (ef1α and β-actin) (Zhong et al., 2008; Zheng and Sun, 2011) as calculated with the 2^–ΔΔCt comparative Ct method. PCR specificity was assessed through melting curve analysis. All data are expressed as the means ± SE (standard error). Post hoc Duncan’s multiple range tests were used to determine the differences between groups using SPSS 16.0. Significance was set at p < 0.05.

Data Availability Statement

The datasets generated for this study can be found in the NCBI Sequence Read Archive BioProject PRJNA639001, PRJNA638896.

Ethics Statement

The study was conducted in accordance with the guidelines and regulations established under the Chinese Government Principles for the Utilization and Care of Animals  Used in Testing, Research, and Training. And the experiments were conducted strictly with the ethical standards of the Guidelines of Regulations for the Administration of Affairs Concerning Experimental Animals documented by the State Science and Technology Commission of Shandong Province. All the researchers who performed the experiments have been appropriately trained. The animal work and animal protocols were approved by the Institute of Oceanology, Chinese Academy of Sciences.

Author Contributions

FY designed the experiments and supervised the study. LW and ZW designed and participated the sample collection. LW, CZ, and SL conducted the experiments. LW, CZ, YZ, and YL performed the data analysis. LW and FY wrote and edited the manuscript. All the authors read and approved the final manuscript.

Funding

This study was supported by the National Key R&D Program of China (2018YFD0900202), the National Natural Science Foundation of China (31772834 and 31702337), the Youth Research Fund of Marine Biology and Biotechnology Laboratory, Pilot National Laboratory for Marine Science and Technology (Qingdao) (YQ2018NO01), the Key R&D (public welfare) program of Shandong Province (2019GHY112005), and the National Infrastructure of Fishery Germplasm Resources (2017DKA30470).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.00751/full#supplementary-material

TABLE S1 | List of the identified RNA editing sites in the flounder gonads.

TABLE S2 | Statistical significance analysis for comparing the editing ratio in different gonadal tissues.

TABLE S3 | GO enrichment analysis of DGEs with differential RNA methylation level.

TABLE S4 | Pathway enrichment analysis of DGEs with differential RNA methylation level.

TABLE S5 | Pathway enrichment analysis of DGEs between the flounder ovary and testis.

TABLE S6 | Primers used in real-time quantitative polymerase chain reactions with their melting temperature.

Footnotes

1. ^ http://code.google.com/p/cutadapt/
2. ^ http://daehwankimlab.github.io/hisat2
3. ^ https://pypi.org/project/HTSeq/
4. ^ https://www.bioinformatics.babraham.ac.uk/projects/fastqc
5. ^ https://github.com/OpenGene/fastp
6. ^ https://bioconductor.org/packages/exomePeak
7. ^ http://www.igv.org
8. ^ http://meme-suite.org
9. ^ http://homer.ucsd.edu/homer/motif
10. ^ https://bioconductor.org/packages/ChIPseeker
11. ^ https://ccb.jhu.edu/software/stringtie
12. ^ https://bioconductor.org/packages/edgeR

References