M6A transcriptome-wide map of circRNAs identified in the testis of normal and AZ-treated Xenopus laevis

Background Evidence showed that N6-methyladenosine (m6A) is strongly associated with male germline development. However, the role of m6A methylation on circRNAs in amphibians remains unknown. In this study, we conducted m6A sequencing analysis to explore the m6A transcriptome-wide profile of circRNAs in testis tissues of Xenopus laevis (X. laevis) with and without treatment with 100 µg/L atrazine (AZ). Results The analysis showed that m6A modification of circRNAs enriched in sense overlapping in testes of X. laevis. We identified the differential m6A modification sites within circRNAs in testes of AZ-exposed X. laevis and compared that with animals from control group. The results showed that a total of 1507 methylated m6A sites was induced by AZ (760 up-methylated and 747 down-methylated). The cross-analysis exhibited a negative correlation of differentially methylated m6A peaks and circRNAs expression level. The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis indicated that 20 key pathways may be involved in the mechanism of testis damage of AZ-exposed X. laevis. Conclusions These findings indicated that differentially m6A-methylated circRNAs may play important roles in abnormal testis development of AZ-exposed X. laevis. This study is the first report about a map of m6A modification of circRNAs in male X. laevis and provides a basis for further studying on the function and mechanism of m6A methylation of circRNAs in the testis development of amphibian. Supplementary Information The online version contains supplementary material available at 10.1186/s41021-023-00279-0.


Introduction
N 6 -methyladenosine (m 6 A), as the most widespread post-transcriptional RNA modification, has attracted extensive attention in the field of epigenetics [1].Emerging evidence suggested that m 6 A modifications played an important role in RNA metabolism, and its dynamic regulation was significantly correlated with gene expression [2].Additionally, a previous study reported that knockout of RNA m 6 A regulators in the testis led to abnormal metabolism of the RNAs, which eventually caused spermatogenetic disorders and infertility [3].Recently, it is reported that the process of spermatogenesis was influenced by gene regulation at transcriptional, posttranscriptional and epigenetic levels [4].RNA m 6 A modification played an important role in spermatogenesis [5].
Circular RNAs (circRNAs) are a class of non-coding RNA with the length more than 200 nucleotides [6].Its closed-loop structures lacking a 5'-cap structure and a 3'-poly (a) tail, were generated by back-splicing of linear mRNAs [7].They were abundant, stable and highly conserved in evolution [8].CircRNAs have pivotally important role in regulating gene expression by sequestering specific miRNAs like a sponge or buffering inhibition of mRNA targets [9].At present, m 6 A modification has been identified in circRNAs, and most m 6 A circRNAs are expressed in a cell type-specific manner suggesting that m 6 A circRNAs play different biological functions in specific cell types [10].For example, Xu et al. found that METTL3/FTO/YTHDF1/2-mediated m 6 A modification can enhance the stability of circRNA-SORE and induce sorafenib resistance in hepatocellular carcinoma (HCC) [11].METTL3-mediated m 6 A methylation regulates the expression of circMETTL3 and promotes breast cancer cell proliferation and migration [12].M 6 A modification improves the expression of circCUX1 and makes hypopharyngeal squamous cell carcinoma radioresistant [13].Mettl14-mediated m 6 A methylation could also export circGFRa1 to cytoplasm from nucleus and promote female germline stem cells (FGSCs) self-renewal [14].Tang et al. found that m 6 A modification modulates the biogenesis of a subset of circRNAs in male germ cells and further regulates spermatogenesis [15].To our knowledge, the m 6 A transcriptome-wide map of circRNAs in amphibians remains to be explored.
Studies have found that environmental pollution could negatively affect amphibians' health status and fertility, ultimately leading to a significant decrease in their population [16][17][18].Atrazine (2-chloro-4-ethylamino-6-isopropylamino-s-triazine, AZ), an endocrine disrupting chemical (EDCs), is a highly effective triazine herbicide widely used in China, the USA and other countries [19].It is often detected in groundwater, surface and drinking water [20].It has been reported that AZ could cause metamorphosis of tadpoles and gonadal dysplasia of in several frog species [21][22][23], such as demasculinization and complete feminization of X. laevis [21].This may be one of the factors causing global amphibian declines [24,25].In our previous study, we examined the damage of X. laevis exposed to AZ (0.1, 10, or 100 µg/L) for 90 days in the water environment.The results showed that AZ caused the decrease of gonad weight and gonad somatic index and the histological damage of testis in X. laevis [26].In addition, we reported that m 6 A transcriptome-wide map of an amphibian species X. laevis.These findings provided clues to reveal the role of m 6 A which may affect amphibian testis development [27].However, the mechanism of reproductive toxicity in male X. laevis chronically exposed to AZ remains unclear.Therefore, it is necessary to investigate the potential changes of m 6 A modification of circRNAs.
Thus, this study aimed to address whether changes in m 6 A modification of circRNAs are associated with abnormal testicular development in X. laevis chronically exposed AZ.We comprehensively analyzed the m 6 A modification profiles of circRNAs in normal and 100 µg/L AZ-exposed X. laevis for 180 days, respectively.Then, we predicted the signal pathways in which dysregulated m 6 A methylation of circRNAs involved by KEGG pathway analysis.Meanwhile, our work will provide a basis for further studying on the function and mechanism of m 6 A methylation of circRNAs in the abnormal testis development of X. laevis and also provide novel insights into its function and biological significance in amphibians.

Ethic approval
All animal experiments were carried out according to relevant guidelines and care regulations.All procedures complied with the "Principles of Animal Care".The protocol was assessed and approved by the Committee on the Ethics of Animal Experiments of Shandong Academy of Occupational Health and Occupational Medicine [28].

Animal treatment
Three pairs of adult male and female X. laevis were purchased from the Chinese Academy of Sciences (Beijing, China).The offspring were generated by natural mating.Laboratory freshwater produced by UV treatment and carbon filtration was used for the acclimatization of frogs in the laboratory and for all subsequent exposures.The X. laevis were kept at an average water temperature of 22 ± 2 °C at pH 7.5, under 12 h light and 12 h dark cycle.Tadpoles were fed fairy shrimp (Artemia nauplii) eggs in a young age daily and pork liver three times per week ad libitum when the tadpoles completed metamorphosis [26].
At Nieuwkoop-Faber (NF) stage 47 (13 d post-hatch), mixed sex tadpoles from the same parents were discretionarily divided into two groups.One group was exposed to AZ at a dose of 100 µg/L (dissolved in solvent vehicle DMSO (0.01%)) for 180 days, and the other group was exposed to 0.01% DMSO only.Test solutions were renewed by 50% replacement every 48 h.The animals were monitored daily for health status and morphological changes.All of the females were excluded in this study.On day 180, the X. laevis in the control and AZ-exposed groups were sacrificed, respectively.The testes were isolated and weighed, and finally stored at − 80 °C for further analysis [8,27].

CircRNAs preparation
Each group was recommended to include at least three biological replicates [29].We randomly selected six testes (three from controls and three from 100 µg/L AZexposed groups) for circRNAs analysis.Total RNA from testes was extracted using Trizol reagent (Invitrogen Corporation, CA, USA) according to the manufacturer's instruction.The concentration and purity of RNAs were measured using a NanoDrop® ND-2000 spectrophotometer (Thermo, Waltham, MA, USA).RNA integrity was examined by denaturing gel electrophoresis experiments [8].RNA samples were further purified and converted to double-stranded cDNA for microarray analysis according to the Agilent ® Xenopus 4 × 44 k Gene Expression Microarray protocols.

CircRNAs m 6 A MeRIP sequencing
M 6 A of circRNAs were sequenced by MeRIP sequencing using Illumina HiSeq sequencer.In short, the fragmented RNA was incubated with anti-m 6 A polyclonal antibody (Synaptic Systems, Goettingen, Germany, 202,003) in IPP buffer for 2 h at 4 °C.Under the same conditions, the mixture was immune precipitated by incubation with protein-A beads (Thermo Fisher).Then, bound RNA was eluted and then extracted according to the manufacturer's instructions.Purified RNA was used for RNA-seq library generation with NEBNext® Ultra™ RNA Library Prep Kit (NEB).Input samples without immune precipitation and m 6 A IP samples were subjected to 150 bp paired-end sequencing on Illumina HiSeq sequencer.Paired-end reads were obtained, and quality control was performed by Q30.More specific operations have been elucidated in our previous study [26,27].

Data analysis
After sequencing, reads were trimmed with 3' adaptor and low-quality readings were removed to produce clean reads with Cutadapt software (v1.9.3).Clean reads from all libraries were aligned to genome by bowtie 2 software [30].Then circRNAs were detected and identified with find circ software [31], the circBase database [32] was also used to annotate the identified circRNAs, clean reads of all libraries were mapped to genome using hisat2 software (v2.04) [33].Methylated sites on circRNAs were detected by MACS software [34], and differentially methylated sites of circRNAs were identified by diffReps [35].In addition, the overlapping sites between the peaks, identified by both software based on the circRNAs exons were chosen for further analyzed.Pathway enrichment analysis of differentially methylated protein coding genes was performed through KEGG pathway database.

AZ induced physiological changes of X. laevis
Compared with froglets in the control, the mortality of froglets significantly increased, the weight and gonadosomatic index (GSI) of the testis were significantly reduced, the metamorphosis time of tadpole was prolonged in AZexposed X. laevis.Histopathological results showed large empty spaces and irregular shape of seminiferous lobules in the testicular tissue.The number of germ cells in testicular tissue was significantly reduced.These results have been reported in our previous paper [26].

M 6 A sites within circRNAs in the testes of control and AZ-exposed X. laevis
The m 6 A transcriptome-wide profile of circRNAs were performed in testes in three biological replicates from the controls (n = 3) and AZ-exposed X. laevis (n = 3).We found that there were 1680 m 6 A peaks were shared in controls and AZ-exposed X. laevis, while 276 m 6 A peaks were identified in controls but absent in AZ-exposed X. laevis, and 1357 m 6 A peaks were identified in AZexposed X. laevis but absent in controls (Fig. 1).

Distribution of m 6 A-methylated circRNAs in the testes of control and AZ-exposed X. laevis
To further illustrate the distribution profiles of m 6 A-methylated circRNAs, we analyzed the genomic position of total m 6 A-methylated circRNAs in the testes of control and AZ-exposed X. laevis, respectively.We found that m 6 A-methylated circRNAs in both groups were divided into 5 groups: antisense, exonic, intergenic, intronic, sense overlapping.Particularly, we found that the most of m 6 A-methylated circRNAs were concentrated in sense overlapping (52.9% in control and 53.1% in AZ-exposed groups) (Fig. 2a and b).

Differentially m 6 A modification sites of circRNAs in X. laevis exposed to 100 µg/L AZ
There were 1507 differentially methylated m 6 A sites were identified.Among them, 760 m 6 A methylated sites were significantly up-regulated, 747 m 6 A methylated sites were significantly down-regulated (Table S1 and Fig. 3).Tables 1 and 2 showed that top ten up-and down-methylated m 6 A sites of circRNAs with the highest fold change (FC) values.To further illustrate the distribution profiles of differentially methylated circRNAs, we found most significantly methylated circRNAs belonged to sense overlapping (Fig. 4a and b).

The enrichment pathways of genes related differentially m 6 A-methylated circRNAs by KEGG
To explore the effect of differential m 6 A-methylated cir-cRNAs in testis of AZ-exposed X. laevis.KEGG pathway analysis were performed for target genes related to the modified circRNAs.The result revealed that target genes related to up-methylated circRNAs were mainly enriched in ten pathways, such as Valine, Leucine and isoleucine degradation, Tryptophan metabolism, Retinol metabolism, MAPK signaling pathway, Linoleic acid metabolism, GnRH signaling pathway, Fatty acid elongation, Fanconi anemia pathway, Carbon metabolism and Calcium signaling pathway (Fig. 6a).Meanwhile, target genes related to down-methylated circRNAs were significantly involved in Tryptophan metabolism, RIG-I-Like receptor signaling pathway, Retinol metabolism, MAPK signaling pathway, Linoleic acid metabolism, Herpes simplex infection, GnRH signaling pathway, Cysteine and methionine metabolism, Butanoate metabolism and ABC transporters (Fig. 6b).

Discussion
Studies have shown that AZ could cause histological damage of the testis, reduced sperm count and generate abnormal sperm [36].In addition, X. laevis exposed to AZ also showed damaged germ cell and reduced testicular weight [37].These results were consistent with the previous results of our study [26].M 6 A methylation has been shown to be one of the most abundant internal modifications of RNA in higher eukaryotes [38,39].It was characterized by dynamic reversibility, wide spready and unique patterns [10].Emerging evidence have suggested that m 6 A methylation played important biological functions in RNA modification, tumorigenesis, fat metabolism, reproductive development, and so on [39,40].And the role of m 6 A circRNAs has also been gradually elucidated.Hence, we explored m 6 A transcriptome-wide map of circRNAs identified in the testis of normal and AZ-treated X. laevis, to predict the role of m 6 A modification of circRNAs in abnormal testicular development.In this study, we investigated the changes of m 6 A modification of circRNAs in the testis of X. laevis by exogenous environmental stimulation.Our findings showed that the number of m 6 A sites within circRNAs increased in the testis of AZ-exposed X. laevis compared with that of controls.The majority of m 6 A methylated circRNAs and differential m 6 A methylated ones enriched in sense overlapping.Interestingly, studies have found that the back splicing tends to occur predominantly in m 6 A-enriched sites in male germ cells.
The change allows the most circRNAs contain large open reading frames to ensure long-term and stable protein production for specific physiological processes in the absence of corresponding linear mRNA [15,41].And aberrant m 6 A methylation may affect the abnormal spermatogenesis [42].Therefore, we speculated that m 6 A modification enriched in sense overlapping may play an important role in abnormal development of testis in amphibian species.
To further explore the potential role of m 6 A modification in the testes of AZ-exposed X. laevis, we performed combined analysis of m 6 A-Seq and RNA-Seq.Results revealed that a negative correlation between differentially methylated peaks and circRNAs expression levels.It is reported that circRNAs were characterized by covalently closed structure and resistance to exonucleases, but m 6 A modification could induce the degradation process of cir-cRNAs by regulating their stability [29].Park et al. found that m 6 A-methylated circRNAs were cleaved by endoribonucleases via a YTHDF2-HRSP12-RNase P/MRP axis [43].This may affect the biological function of circRNAs [44], such as m 6 A modification of circNSUN2 could promote the liver metastasis of colorectal cancer by forming a circNSUN2/IGF2BP2/HMGA2 RNA-protein ternary complex [12].Additionally, it has been reported that the specific role of m 6 A modification on gene expression largely depended on the type of m 6 A 'readers' [45].Hence, we speculated that m 6 A modification in circRNAs may play a potential role in the abnormal We further analyzed the role of differentially m 6 A-methylated circRNAs related target genes in AZ-treated X. laevis for which 20 key pathways were obtained by KEGG pathway annotation method.In this study, the most GeneRatio term was "MAPK signaling Fig. 6 The annotated significant pathways targeted by the enrichment score of the differentially m 6 A-methylated circRNAs-related genes including upmethylated (a) and down-methylated (b) in testis of X. laevis exposed to 100 µg/L AZ.The horizontal axis is the -log10 (P-value) for the pathway and the vertical axis is the pathway category pathway".Studies have shown that some exogenous stimulants could activate MAPK signaling pathway, affected Sertoli cells (SCs) proliferation and blood-testis barrier (BTB) structure [46,47].Interestingly, In the process of spermatogenesis, SCs not only provided mechanical and nutritional support [48,49], but also protected germ cells by forming an immune protective environment through the BTB [50,51].Therefore, we assumed that MAPK signaling pathway regulated by circRNAs with m 6 A modifications might mediate the abnormal testis development of AZ-exposed X. laevis.
"GnRH signaling pathway" was also included in the result of KEGG.GnRH, as a local bioregulator, controlled the secretion of gonadotropin of the pituitary gland, including follicle stimulating hormone and luteinizing hormone, which were critical components affecting sperm development [52].In addition, it was reported that GnRH could also be synthesized in seminiferous tubules, and as a paracrine mediator of spermatogenesis [53].In GnRH deficient males, there would be incomplete testicular growth and maturation, cryptorchidism, androgen deficiency and other symptoms [54].Therefore, m 6 A-methylation of circRNAs involved in "GnRH signaling pathway" may play an important role in the growing development and functional activity of testis of AZ-exposed X. laevis.
In addition, it has been shown that Ca 2+ , which was produced by the Calcium Signaling Pathway, could directly control metabolism, secretion, fertilization, proliferation and other processes [55].And it also participated in the cAMP-induced steroidogenesis in Leydig cells [56].Therefore, we speculated that "Calcium Signaling Pathway" regulated by m 6 A-methylated circRNAs may involve testicular dysplasia of AZ-exposed X. laevis.Most of these pathways in the result were related to metabolism such as "Linoleic acid metabolism" and "Fatty acid elongation" signaling pathways.Fatty acid (FA) mainly included monounsaturated fatty acids (MUFA) and polyunsaturated fatty acids (PUFA) [57], while Linoleic acid were essential PUFA [58].Its synthesis was mediated by rate-limiting elongation and desaturation enzymes, which were regulated by fatty acid desaturase 2 (FADS2) and elongation of very-long-chain fatty acidslike 2 (ELOVL2) [58].It has been found that arrested spermatogenesis and lacked mature spermatozoa in the testis of mice with FADS2 and ELOVL2 gene knockout [59,60].Additionally, it was reported that the FA metabolic disorders could occur in infertile men [57].Hence, the role of m 6 A-methylated circRNAs involved in "Linoleic acid metabolism" and "Fatty acid elongation" signaling pathways deserve further study in male reproductive damage.Consequently, the results predicted that m 6 A methylation of circRNAs may regulate the differential expression of these target genes in abnormal testis development of X. laevis exposed to 100 µg/L AZ.

Conclusion
We detected the m 6 A transcriptome-wide profile of cir-cRNAs in the testes of control and AZ-exposed X. laevis.The results showed that AZ could alter expression profile in 1507 m 6 A methylated peaks within circRNAs in which 760 were significantly up-methylated and 747 significantly down-methylated, and they mainly enriched in sense overlapping.Conjoint analysis indicated that a negative correlation of differentially methylated m 6 A peaks and circRNAs expression level, suggesting a regulatory role of m 6 A modification of circRNAs in amphibious gene expression.KEGG pathway analysis revealed that m 6 A related to "MAPK signaling pathway", "GnRH signaling pathway", "Calcium signaling pathway" and pathways related to metabolism such as "Linoleic acid metabolism" and "Fatty acid elongation" maybe play a pivotal role in the abnormal testis development of AZexposed X. laevis.Our study provided a basis for further studies on the function and mechanism of m 6 A methylation of circRNAs in the abnormal testis development of X. laevis.Meanwhile, this study presented the first m 6 A transcriptome-wide map of circRNAs in amphibian species X. laevis.This may help to further understand the role of m 6 A methylation in testis development and spermatogenesis in amphibians.

Fig. 1
Fig. 1 Venn diagram showing the overlap of m 6 A peaks within circRNAs in the testes of control and AZ-exposed X. laevis

Fig. 3 Fig. 2
Fig. 3 Volcano plots showing -log10 (P_value) versus log2FC in m 6 A methylated sites on circRNAs in the testes of control and AZ-exposed X. laevis.Red circles denote significantly up-regulated m 6 A methylated sites, whereas blue circles denote significantly down-regulated m 6 A methylated sites (p < 0.05 and fold change ≥ 4)

Fig. 4
Fig. 4 Overview the distribution of differentially m 6 A modification sites of circRNAs.a: Pie charts showing the percentage of up-methylated m 6 A in five segments.b: Pie charts showing the percentage of down-methylated m 6 A peaks in five segments

Fig. 5 a
Fig. 5 a: Dot plot of log2FC (circRNAs expression) against log2FC (differential m 6 A methylation) showing a negative correlation between overall m 6 A methylation and mRNA expression level (P = 0.02; Pearson R = -0.25).b: Four quadrant plots showing gene expression with differentially methylated m 6 A peaks

Table 1
The top 10 up-methylated m6A sites of circRNAs

Table 2
The top 10 down-methylated m6A sites of circRNAs txStart/txEnd: Start/end position of the differentially methylated RNA sites