Alterations of the m6A Methylation Induced by TGF-β2 in ARPE-19 Cells

Background: N ${}^{6}$ -methyladenosine (m ${}^{6}$ A) participates in diverse physiological processes and contributes to many pathological conditions. Epithelial-mesenchymal transition (EMT) of retinal pigmental epithelial (RPE) cells plays an essential role in retinal-related diseases, and transforming growth factor $\beta{}$ 2 (TGF- $\beta{}$ 2) is known to induce EMT in vitro. However, the effect of TGF- $\beta{}$ 2 on m ${}^{6}$ A methylation in RPE cells is not yet known. Methods: RNA-seq and MeRIP-seq were performed to analyze changes at the mRNA and m ${}^{6}$ A levels after TGF- $\beta{}$ 2 treatment of human ARPE-19 cells. mRNA levels and total m ${}^{6}$ A levels were subsequently validated. Results: Sequencing revealed 929 differentially expressed genes and 7328 differentially methylated genes after TGF- $\beta{}$ 2 treatment. Conjoint analysis identified 290 genes related to microtubule cytoskeleton, focal adhesion, ECM-receptor interaction, cell division, cell cycle, AGE-RAGE, PI3K-Akt and cGMP-PKG pathways. Further analysis revealed that 12 EMT-related genes were altered at the mRNA and m ${}^{6}$ A levels after TGF- $\beta{}$ 2 treatment (CALD1, CDH2, FN1, MMP2, SPARC, KRT7, CLDN3, ELF3, FGF1, LOXL2, SHROOM3 and TGFBI). Moreover, the total m ${}^{6}$ A level was also reduced. Conclusions: This study revealed the transcriptional profiling of m ${}^{6}$ A modification induced by TGF- $\beta{}$ 2 in RPE cells. Novel connections were discovered between m ${}^{6}$ A modification and TGF- $\beta{}$ 2-induced EMT, suggesting that m ${}^{6}$ A may play crucial roles in the EMT process.

Keywords

m⁶A methylation

RPE cells

EMT

TGF-β2

transcriptomics

MeRIP-seq

1. Introduction

Retinal pigment epithelium (RPE) cells play an important role in the development of proliferative vitreoretinopathy (PVR). This isan abnormal repair process caused by failed retinal surgery for retinal detachment or serious ocular trauma [1]. During PVR, the surface of the nerve retina forms a fibrous membrane due to inflammation and wound healing, resulting in structural deformation of the retina. The fibrous membranes mainly include RPE cells, macrophages, myofibroblast-like cells and glial cells [2]. The fibroblast-like cells are derived from RPE cells via epithelial to mesenchymal transition (EMT) and are the underlying cause of PVR [3, 4, 5]. Mature RPE cells are mitotically quiescent under normal physiological conditions. However, various growth factors and cytokines are released from the serum when the blood-retinal barrier is damaged, thereby inducing RPE cells to enter the EMT process [6]. The essential growth factors and cytokines are primarily composed of transforming growth factor $\beta{}$ (TGF- $\beta{}$ ), fibroblast growth factors (FGF), insulin-like growth factors (IGF), platelet-derived growth factors (PDGF), epidermal growth factor (EGF), tumor necrosis factor alpha (TNF- $\alpha{}$ ), interleukin-1 $\beta{}$ (IL-1 $\beta{}$ ), IL-6 and IL-8 [7]. These factors cause the cells to acquire mesenchymal phenotypes with enhanced proliferation, migration, and invasion capabilities. TGF- $\beta{}$ 2 is the major regulator of EMT in vitro, with numerous investigators confirming that it induces EMT in RPE cells [8, 9, 10]. Moreover, ARPE-19 cells exhibit enhanced migratory capacity and increased expression of $\alpha{}$ -smooth muscle actin ( $\alpha{}$ -SMA) in response to TGF- $\beta{}$ 2 treatment [11, 12]. Although molecular targets of the EMT process have been described, other key factors in the EMT of RPE cells remain to be elucidated. These include RNA methylation, which may play crucial roles in EMT development.

N ${}^{6}$ -methyladenosine (m ${}^{6}$ A) is an abundant, ubiquitous and conserved modification of transcripts in higher eukaryotes that participates in gene regulation through multiple mRNA processes [13, 14]. It is a dynamic and reversible modification process in which mRNA is modified by methyltransferases (including METTL3, METTL14, METTL16, WTAP, RBM15, RBM15B and KIAA1429), reversed by demethylases (including FTO and ALKBH5), and recognized by m ${}^{6}$ A binding proteins (including YTHDC1-2, YTHDF1-3, IGF2BP1-3 and HNRNPA2B1) [14]. Substantial evidence shows that m ${}^{6}$ A methylation participates in the control of RNA transcription [15], processing [16], splicing [17, 18], degradation [19, 20] and translation [21, 22]. Studies have also shown that m ${}^{6}$ A may have different roles in different cell types. For example, high METTL3 expression and high m ${}^{6}$ A levels facilitate EMT in lung cancer cells [23], suggesting that hypermethylation is positively correlated with EMT. However, high ALKBH5 expression and low m ${}^{6}$ A levels promote the development of EMT in glioblastoma, indicating that hypomethylation is positively correlated with EMT in these cells [24]. Therefore, m ${}^{6}$ A modifications may have multiple roles in EMT and complex regulatory networks that involve diverse target genes. m ${}^{6}$ A methylation has been extensively studied in various tumor cell types, but less in retina-associated cells. Increased METTL3 has been reported to reduce apoptosis and pyroptosis triggered by high glucose and inhibit RPE cell proliferation [25]. METTL3 expression is reduced in PVR membranes and RPE cells undergoing EMT, with further studies revealing the Wnt/ $\beta{}$ -catenin signaling pathway plays a role in the inhibition of EMT by METTL3 [26]. In addition, METTL14 participates in regulating the methylation of microtubule-associated protein (MAP)2, leading to an imbalance in neuronal differentiation (NEUROD1) and affecting the activity of RPE cells [27].

Transcriptome-wide m ${}^{6}$ A modifications have rarely been described in RPE cells, and the connection between m ${}^{6}$ A modifications and EMT is not fully understood. Here, RNA-seq and MeRIP-seq were performed to analyze the effects of TGF- $\beta{}$ 2 treatment on mRNA expression and m ${}^{6}$ A modification respectively in RPE cells. This allowed us to analyze the profile of m ${}^{6}$ A methylation and explore the mechanism by which m ${}^{6}$ A methylation regulates EMT in RPE cells.

2. Materials and Methods

2.1 Cell Lines and Culture Conditions

Human ARPE-19 cells were purchased from ATCC (CRL-2302, Manassas, VA, USA) and routinely grown in DMEM medium (Solaibio, 12100, Beijing, China) containing 100 U/mL penicillin, 100 µg/mL streptomycin, and 10% fetal bovine serum (Gibco, 10099141C, Carlsbad, CA, USA). After the cells adhered, the medium was replaced with fresh DMEM medium (1% FBS) containing double antibiotics. No recombinant protein was added to the control group, whereas the TGF- $\beta{}$ 2 group was treated with 20 ng/mL TGF- $\beta{}$ 2 (MedChemExpress, P61812, Monmouth Junction, NJ, USA). Cells in each group were cultured for 72 hours, and three parallel samples were prepared for each group. The ARPE-19 cells used in this study have been identified by short tandem repeat sequence (STR) analysis and mycoplasma testing was performed.

2.2 Isolation of Total RNA

The extraction of total intracellular RNA was performed using TRIzol reagent (Invitrogen, 15596026, Carlsbad, CA, USA) as previously described [28]. The extracted RNA was digested with DNase I (Takara, 2270A, Tokyo, Japan). After the quality, integrity and concentration of the total RNA had been determined, it was sent to Wuhan Seqhealth Technology (Wuhan, China) for sequencing analysis.

2.3 Library Construction and Sequencing

2 µg of total RNA was used for RNA-seq and 50 µg for MeRIP-seq. Polyadenylated RNA was enriched using Oligo d(T) beads (VAHTS, N401-01, Nanjing, China). The mRNAs were broken into 100–200 nt fragments by 20 mM ZnCl ${}_{2}$ . Ten percent of the RNA fragments were set as the input group, while the rest were set as the IP group for immunoprecipitation with m ${}^{6}$ A antibody (Synaptic Systems, 202203, Gottingen, German). The library was constructed using the KC-Digital ${}^{\text{TM}}$ Stranded mRNA Library Prep Kit (Seqhealth, DR08502, Wuhan, China), and the library products were then sequenced on a Novaseq 6000 sequencer (Illumina) [29]. The raw sequencing data were uploaded to the NCBI SRA database (SRA accession: SRP422412, Bioproject: PRJNA934747).

2.4 Analysis of Sequencing Data

The raw data were analyzed for quality with FastQC (version 0.11.5, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and filtered with Trimmomatic (version 0.36, http://www.usadellab.org/cms/index.php?page=trimmomatic) [30]. Clean reads were clustered by UMI tools (version 1.0, default parameters, https://gitHub.com/CGATOxford/UMI-tools) and each sub-cluster was aligned multiple times to eliminate biases and errors [31]. Next, they were aligned with the Homo sapiens genome assembly GRCh38.p13 using STAR software (version 2.5.3a, default parameters, https://github.com/alexdobin/STAR). RSeQC (version 2.6, https://rseqc.sourceforge.net/) was used to assess the distribution, uniformity of coverage and strand specificity of mapped reads [32]. exomePeak (version 3.8, https://bioconductor.org/packages/3.8/bioc/html/exomePeak.html) was used to predict high-intensity and high-resolution m ${}^{6}$ A peaks (peak calling), which were further annotated using bedtools (version 2.25.0, https://bedtools.readthedocs.io/en/latest/index.html) [33]. Motif analysis was performed by Homer (version 4.10, http://homer.ucsd.edu/homer/index.html) to obtain enriched motifs. deepTools (version 2.4.1, http://deeptools.ie-freiburg.mpg.de) was applied to analyze the distribution of m ${}^{6}$ A peaks in each functional region of all transcripts, while bedtools was used to analyze the distribution of m ${}^{6}$ A peaks on the corresponding chromosome [34, 35].

The reads of all genes annotated in the genome were calculated by featureCounts (version 1.5.1, http://subread.sourceforge.net). The RPKM value served as a measure of gene expression and was calculated using the formula shown below. Differentially expressed genes (DEGs) and differentially methylated genes (DMGs) between the two experimental groups were screened by the edgeR package (version 3.12.1, http://bioconductor.org/packages/3.2/bioc/html/edgeR.html) with p-value $<$ 0.05 and $|{}$ fold change $|{}$ $>$ 1.5 [36].

$\text{ RPKM }=\frac{\text{ total exon reads }}{\text{ mapped reads (millions) % }\times\text{ exon length }(\mathrm{KB})}$

2.5 Enrichment Analysis and Graphics Production

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analyses were performed using Omicshare tools to determine the enrichment pathways for DEGs and DMGs. A p-value $<$ 0.05 indicated statistical significance. PCA and volcano graphs were prepared by inputting all read counts for each sample, while the nine-quadrant graph was plotted by inputting the fold change of the DEGs and DMGs. These were performed by Omicshare tools set at the default parameters.

2.6 Determination of Total Intracellular m

{}^{6}

A Levels

Total m ${}^{6}$ A methylation was measured using an m ${}^{6}$ A methylation quantification kit (EpiGentek, P-9005, Farmingdale, NY, USA) and 300 ng of total RNA per reaction. Both negative and positive RNA controls were provided with the kit. The absorbance of each sample at 450 nm was measured with a spectrophotometer (Thermo Scientific, Waltham, MA, USA), and the percentage of m ${}^{6}$ A in the total RNA was calculated [37]. Two replicate wells were analyzed for each sample and the experiments were repeated three times.

2.7 Real-Time Quantitative PCR (qRT-PCR)

qRT-PCR was performed to verify the level of mRNA as previously described [28]. cDNA was synthesized using a Reverse Transcription Kit (Takara, RR047A, Tokyo, Japan) and qRT-PCR was performed using a StepOnePlus PCR System with TB Green reagent (Takara, RR820A, Tokyo, Japan). The relative mRNA levels of target genes were calculated using the 2 ${}^{-\Delta\Delta\text{CT}}$ algorithm, with the glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene used as the reference gene. The primer sequences are listed in Supplementary Table 1.

2.8 Statistical Analysis

Three replicates were set for sequencing, m ${}^{6}$ A quantitative assay and qRT-PCR. SPSS software (version 21.0, IBM Corp., Armonk, NY, USA) was used for statistical analysis, with differences between the control and TGF- $\beta{}$ 2 groups identified by the student’s t-test. A p-value $<$ 0.05 indicated statistical significance.

3. Results

3.1 Analysis of DEGs Induced by TGF-

\beta{}

Fig. 1A shows the gene expression clustering map for each group obtained after RNA-seq analysis. To identify TGF- $\beta{}$ 2-induced changes at mRNA levels, a $|{}$ fold change $|{}$ $>$ 1.5 and p-value $<$ 0.05 were used as the threshold for DEGs. As shown in Fig. 1B,C, 929 DEGs were identified, of which 581 genes showed increased expression and 348 genes showed decreased expression (Supplementary Table 2). Heatmap analysis showed consistency in the gene expression differences observed for the three replicate samples (Fig. 1D).

Fig. 1.

Analysis of DEGs induced by TGF- $\beta{}$ 2 in ARPE-19 cells. (A) PCA clustering map of the gene expression profile in each group. (B) Volcano plot analysis of DEGs. Blue dots represent down-regulated genes, grey dots represent genes with no difference in expression, and red dots represent up-regulated genes. (C) Number of DEGs with increased or decreased expression. (D) Heatmap analysis of DEGs. Red indicates hyper-expressed genes and blue indicates hypo-expressed genes. (E) GO enrichment analysis of DEGs. (F) KEGG enrichment analysis of DEGs. The dot size indicates the number of genes enriched, and the dot color indicates significant enrichment.

GO enrichment analysis revealed the DEGs were primarily related to cell division, the cell cycle, nuclear division, extracellular matrix organization and structural constituent, collagen-containing extracellular matrix, and chromosome segregation (Fig. 1E). KEGG analysis showed the DEGs were involved in the cell cycle, focal adhesion, ECM-receptor interaction, protein digestion and absorption, PI3K-Akt signaling, p53 signaling, and the AGE-RAGE signaling pathway (Fig. 1F).

3.2 Global m

{}^{6}

A Modification Patterns in RPE Cells

We next performed transcriptome-wide MeRIP-seq in the control and TGF- $\beta{}$ 2-treated groups in order to explore the modification of m ${}^{6}$ A methylation in RPE cells. As shown in Fig. 2, 10,025 m ${}^{6}$ A peaks were found in the control group and 10,156 in the TGF- $\beta{}$ 2 group, with 8968 overlapping peaks between the two groups (Supplementary Tables 3,4). Compared with the control group, there were 1188 extra peaks and 1057 missing peaks in the TGF- $\beta{}$ 2 group (Fig. 2A). After annotating the peaks, 7626 genes were found to be m ${}^{6}$ A methylated in the control group, 7685 genes were m ${}^{6}$ A methylated in the TGF- $\beta{}$ 2 group, and 7183 genes were in the intersection between the two groups (Fig. 2B). Moreover, m ${}^{6}$ A methylation frequently occurred within RRACH motifs (Fig. 2C), revealing the m ${}^{6}$ A consensus motif in RPE cells.

Fig. 2.

m ${}^{6}$ A modification patterns in ARPE-19 cells. (A) Venn analysis of m ${}^{6}$ A peaks in the control and TGF- $\beta{}$ 2 groups. (B) Venn analysis of genes modified by m ${}^{6}$ A methylation in the control and TGF- $\beta{}$ 2 groups. (C) Conserved m ${}^{6}$ A consensus motif in the control and TGF- $\beta{}$ 2 groups. (D) Distribution of m ${}^{6}$ A methylation at different sites in transcripts. (E) Distribution of m ${}^{6}$ A peaks in each functional region of genes. 5′UTR, 5′ untranslated region; CDS, coding sequence; 3′UTR, 3′ untranslated region; ncDNA-T, transcripts of noncoding DNA. (F) Distribution of m ${}^{6}$ A methylation in transcripts from different chromosomes. MT represents mitochondrial chromosome.

TGF- $\beta{}$ 2 did not affect the preference for m ${}^{6}$ A-modified regions, which were enriched near stop codons (Fig. 2D). m ${}^{6}$ A methylation in the control and TGF- $\beta{}$ 2 groups was enriched in the CDS, introns and 3′UTR regions of transcripts, with proportions of 23.62%/23.64%, 48.53%/47.83, and 23.66%/24.44%, respectively (Fig. 2E). By counting the number of m ${}^{6}$ A peaks in ARPE-19 cells, it was discovered that m ${}^{6}$ A modifications occurred in transcripts derived from each chromosome. Among them, the number of m ${}^{6}$ A peaks in transcripts from chromosome 1 was 1626 and 1651 in the control group and TGF- $\beta{}$ 2 group, respectively; while chromosome Y was no more than 10 (Fig. 2F). After normalization of chromosome length, we found that m ${}^{6}$ A peaks were the highest on mitochondrial chromosome, followed by chromosome 19, and the smallest on chromosome Y (Supplementary Fig. 1).

3.3 Analysis of DMGs Induced by TGF-

\beta{}

To evaluate the effect of TGF- $\beta{}$ 2 treatment at the m ${}^{6}$ A methylation level of individual genes, we set $|{}$ fold change $|{}$ $>$ 1.5 and p-value $<$ 0.05 to define DMGs. Only 9 genes were found to contain hypermethylated peaks, whereas 7325 genes contained hypomethylated peaks (Fig. 3A, Supplementary Table 5). Table 1 lists the top five genes with hypermethylated peaks (TRO, SLF1, RBM18, MIEF2, CCNF) and the top five genes with hypomethylated peaks (SMIM10L2B, PFKP, TTLL11, RBMS3, BDKRB2). These genes are involved in cell adhesion, regulation of cell cycle, RNA binding and the bradykinin system. DMGs were present on each chromosome, with 754 on chromosome 1 and 5 on chromosome Y (Fig. 3B). According to the normalization of chromosome length, the number of DMGs was the highest on mitochondrial chromosome and the smallest on chromosome Y (Supplementary Fig. 2).

Fig. 3.

Analysis of DMGs induced by TGF- $\beta{}$ 2 in ARPE-19 cells. (A) Number of genes with hypermethylated or hypomethylated peaks. (B) Distribution of DMGs on different chromosomes. (C) GO enrichment analysis of DMGs. (D) KEGG enrichment analysis of DMGs.

Table 1.Top five genes with hypermethylated peaks and top five genes with hypomethylated peaks.

Gene ID	Gene description	Chromosome	Start	End	log ${}_{2}$ FC	p-value	Annotation
TRO	Trophinin, mediates cell adhesion.	X	54930557	54930632	3.1463	0.0211	CDS
SLF1	Smc5/6 localization factor 1, involved in DNA repair.	5	94694906	94695031	1.6002	0.0231	CDS
RBM18	RNA-binding motif protein 18.	9	122238796	122239021	1.1276	0.0322	3′UTR
MIEF2	Mitochondrial elongation factor 2, regulates mitochondrial morphology.	17	18265639	18265689	1.1256	0.0094	3′UTR
CCNF	Cyclin-F, is involved in the regulation of cell cycle and ubiquitination.	16	2457725	2457850	1.0333	0.0055	3′UTR
SMIM10L2B	Small integral membrane protein 10-like protein 2B.	X	135097034	135097084	–3.2745	0	3′UTR
PFKP	Phosphofructokinase, participates in glycolysis.	10	3066332	3066407	–2.8447	0.0014	5′UTR
TTLL11	Tubulin polyglutamylase TTLL11.	9	121989411	121989461	–2.7538	0.0144	CDS
RBMS3	RNA binding motif single stranded interacting protein 3.	3	30004611	30004661	–2.7343	0.0002	3′UTR
BDKRB2	B2 bradykinin receptor, participates in activating the phosphatidylinositol-calcium second messenger system.	14	96242121	96242296	–2.6528	0.0001	3′UTR

GO enrichment analysis revealed the DMGs were associated with the metabolic process and RNA binding (Fig. 3C). KGEE analysis found the DMGs were mainly enriched in cell junction processes (including adherens junction, tight junction and focal adhesion), the RNA processing pathway (including spliceosome and RNA transport), apoptosis, regulation of actin cytoskeleton, and signaling pathways for PI3K-Akt, mTOR, MAPK and AMPK (Fig. 3D).

3.4 Conjoint Analysis of RNA-seq and MeRIP-seq

Conjoint analysis of DEGs and DMGs was performed to identify genes with altered mRNA and m ${}^{6}$ A modification. As shown in Fig. 4A, one gene showed up-regulated mRNA and hypermethylated peaks, 215 genes had up-regulated mRNA and hypomethylated peaks, and 75 genes had down-regulated mRNA and hypomethylated peaks. Venn diagram analysis showed that 290 genes were altered at the m ${}^{6}$ A and RNA levels (Fig. 4B, Supplementary Table 6). Comparing Fig. 4A,B, the CCNF gene was found to be up-regulated at the transcriptional level, but contained both hypermethylated and hypomethylated peaks.

Fig. 4.

Conjoint analysis of DEGs and DMGs. (A) Nine-quadrant plot of DEGs and DMGs. The abscissa represents DMGs and the ordinate represents DEGs. Red dots represent genes up-regulated at the mRNA and m ${}^{6}$ A levels, yellow dots represent genes up-regulated at the mRNA level and down-regulated at the m ${}^{6}$ A level, and green dots represent genes down-regulated at both the mRNA and m ${}^{6}$ A levels. (B) Venn diagram of DEGs and DMGs. The p-value is the hypergeometric distribution, and all genes detected by sequencing was set as the total number. (C) GO enrichment analysis of genes with altered mRNA level and m ${}^{6}$ A modification. (D) KEGG enrichment analysis of genes with altered mRNA level and m ${}^{6}$ A modification.

GO enrichment analysis revealed the conjoint genes were associated with cell division, the cell cycle, chromosome segregation, and the microtubule cytoskeleton (Fig. 4C). KEGG enrichment analysis showed the conjoint genes were related to the cell cycle, apoptosis, focal adhesion, ECM-receptor interaction, protein digestion and absorption, and the signaling pathways for AGE-RAGE, PI3K-Akt, cGMP-PKG, and p53 (Fig. 4D).

3.5 Analysis of Conjoint Genes and EMT-Related Genes

To investigate the connection between m ${}^{6}$ A methylation and EMT, we screened for EMT-related genes that showed changes at both m ${}^{6}$ A and mRNA levels after TGF- $\beta{}$ 2 induction. We retrieved 314 EMT-related genes through a literature review (Supplementary Table 7) [38, 39]. By cross-matching 290 conjoint genes with 314 EMT-related genes, 12 EMT-related genes were found to be altered at both the mRNA and m ${}^{6}$ A levels (Fig. 5A), including CALD1, CDH2, FN1, MMP2, SPARC, KRT7, CLDN3, ELF3, FGF1, LOXL2, SHROOM3, and TGFBI (Table 2). All 12 genes were up-regulated at the transcriptional level and hypomethylated at the m ${}^{6}$ A level, speculating that EMT gene expression may be regulated by m ${}^{6}$ A methylation in RPE cells.

Fig. 5.

Validation of EMT-related gene expression in ARPE-19 cells. (A) Overlap between DMGs-DEGs and EMT-related genes. (B) Total m ${}^{6}$ A level in the control and TGF- $\beta{}$ 2 groups. (C) qRT-PCR validation of the mRNA levels for EMT-related genes. Asterisks represent statistical significance (*p-value $<$ 0.05, t test).

Table 2.EMT-related genes with increased mRNA level and m ${}^{6}$ A hypomethylation.

Gene ID	Gene Description	m ${}^{6}$ A level	log ${}_{2}$ FC (m ${}^{6}$ A)	m ${}^{6}$ A sites	mRNA level	log ${}_{2}$ FC (mRNA)
CALD1	Caldesmon 1, participates in smooth muscle contraction.	Down	–0.8467	5′UTR	Up	0.6046
CDH2	Cadherin 2, is involved in homotypic cell adhesion.	Down	–0.8173	CDS	Up	0.6275
FN1	Fibronectin 1, mediates cell adhesion and cell motility.	Down	–0.7624	CDS	Up	0.9369
MMP2	Matrix metallopeptidase 2, participates in vasculature remodeling, cell invasion and inflammation.	Down	–0.7004	3′UTR	Up	1.7231
SPARC	Secreted protein acidic and cysteine rich, affects cell growth.	Down	–0.7842	5′UTR	Up	0.6247
KRT7	Keratin 7, facilitates DNA synthesis.	Down	–0.7546	CDS	Up	0.8805
CLDN3	Claudin-3, participates in tight junction.	Down	–1.1295	3′UTR	Up	0.6724
ELF3	E74 like ETS transcription factor 3, transcriptional activator.	Down	–0.8627	3′UTR	Up	1.2989
FGF1	Fibroblast growth factor 1, regulates cell division, differentiation, migration and angiogenesis.	Down	–0.9726	3′UTR	Up	1.2966
LOXL2	Lysyl oxidase like 2, interacts with SNAI1 and inhibits CDH1 expression.	Down	–0.7532	3′UTR	Up	0.6011
SHROOM3	Shroom family member 3, promotes apicobasal cell elongation and regulates cell morphology.	Down	–1.4456	CDS	Up	0.5912
TGFBI	TGF- $\beta$ induced protein, participates in cell adhesion.	Down	–0.6216	3′UTR	Up	0.9588

The total intracellular m ${}^{6}$ A level was reduced in RPE cells following treatment with TGF- $\beta{}$ 2, as compared to untreated cells (Fig. 5B). In addition, qRT-PCR of the 12 EMT-related genes confirmed that CALD1, CDH2, FN1, MMP2, SPARC, FGF1, LOXL2, SHROOM3, and TGFBI were significantly up-regulated (Fig. 5C).

4. Discussion

m ${}^{6}$ A methylation is an extensive intracellular RNA modification and part of a complex regulatory network that affects gene expression in various physiological conditions [14]. As mentioned above, the EMT process in RPE cells is crucial for the formation of fibrous membranes. TGF- $\beta{}$ 2 is an important EMT-inducing factor. To investigate the connection between m ${}^{6}$ A methylation and EMT in RPE cells, we performed transcriptome-wide m ${}^{6}$ A-modified mapping of TGF- $\beta{}$ 2-induced RPE cells. Consistent with previous studies, m ${}^{6}$ A modification was prevalent in RPE cells. MeRIP-seq revealed 10,025 m ${}^{6}$ A peaks in the control group and 10,156 peaks in the TGF- $\beta{}$ 2 group, with m ${}^{6}$ A modifications present in transcripts from each chromosome (Fig. 2A,F). m ${}^{6}$ A methylation is conserved at modification sites in different regions of transcripts, preferentially near stop codons and within 3′UTR [40]. Furthermore, 35% of m ${}^{6}$ A modifications are located in the CDS region, where they can positively regulate translation by resolving mRNA secondary structure [41]. The present study found that m ${}^{6}$ A modifications were abundant in the CDS, introns and 3′UTR regions of transcripts, as well as being centered around stop codons (Fig. 2D,E). The METTL3-METTL14 complex preferentially recognizes m ${}^{6}$ A methylation sites in RRACH motifs in mammalian cells [42]. Consistent with this, our study also confirmed that m ${}^{6}$ A methylation identified the RRACH motif in RPE cells (Fig. 2C). These results indicate that m ${}^{6}$ A methylation is a prevalent and conserved RNA modification in RPE cells.

Comparison of the methylation changes between control and TGF- $\beta{}$ 2 groups identified 7334 DMGs, of which only 9 genes (0.12%) had hypermethylated m ${}^{6}$ A peaks (Fig. 3A). In accordance with this observation, the total m ${}^{6}$ A methylation level of intracellular mRNA decreased after treatment with TGF- $\beta{}$ 2 (Fig. 5B). Thus, we speculating that m ${}^{6}$ A methylation is negatively correlated with the EMT process in RPE cells. KEGG analysis found the DMGs were enriched in adherens junction, tight junction, focal adhesion, and regulation of actin cytoskeleton (Fig. 3D). These affect epithelial cell polarity [43], cell migration [44], and cell homeostasis [45]. The DMG-enriched pathways include PI3K/Akt signaling, which mediates cell cycle, cell growth, apoptosis, metabolism, angiogenesis, migration, differentiation and EMT [46]. mTOR signaling also participates in the regulation of EMT by altering the protein kinase C (PKC) phosphorylation level and by directly phosphorylating and activating Akt [47]. These observations suggest a potential relationship between m ${}^{6}$ A methylation and EMT in RPE cells.

Conjoint analysis of DMGs and DEGs revealed that 290 genes were simultaneously altered at the m ${}^{6}$ A and RNA levels. These conjoint genes were enriched in ECM-receptor interactions, AGE-RAGE signaling, PI3K-Akt signaling, cGMP-PKG signaling, and the p53 signaling pathway (Fig. 4B,D). ECM-receptor interactions have profound effects on major cellular programs including growth, differentiation, migration and survival, and have been associated with the EMT process in RPE cells [48]. AGE-RAGE signaling enhances NF-kB signaling, stimulates the generation of reactive oxygen species, and accelerates the EMT process induced by TGF- $\beta{}$ [49]. cGMP/PKG signaling is also involved in EMT. For example, atrial natriuretic peptide inhibits the EMT process through the cGMP/PKG pathway [50]. In addition, TGF- $\beta{}$ 1 has a profibrotic effect by inhibiting soluble guanylate cyclase-cGMP-PKG signaling in systemic sclerosis [51]. p53 signaling is also closely related to the EMT process, with wild-type p53 negatively regulating EMT initiation and metastasis, whereas mutant p53 promotes EMT and metastasis [52]. These findings suggest that m ${}^{6}$ A methylation could regulate EMT through multiple signaling pathways.

Twelve EMT-related genes (CALD1, CDH2, FN1, MMP2, SPARC, KRT7, CLDN3, ELF3, FGF1, LOXL2, SHROOM3, TGFBI) showed both increased mRNA levels and hypomethylation after TGF- $\beta{}$ 2 treatment (Table 2). CALD1, CDH2, FN1, MMP2 and SPARC are known to be highly expressed in EMT [38], and most can be used as prognostic biomarkers for various cancers. Among them, fibronectin (FN) encoded by FN1 is essential in the EMT process and fibrosis. Massive deposition of FN, a component of the extracellular matrix, is known to promote tumor cell growth, migration, invasion and angiogenesis [53]. Furthermore, FN mediates the initial process of fibrous membrane formation, and intravitreal FN concentrations increase with the clinical stage of PVR evolution [54]. The expression and phosphorylation of caldesmon, encoded by CALD1, is up-regulated after TGF- $\beta{}$ treatment, further increasing the number and size of focal adhesions and cell contractility [55]. N-cadherin encoded by CDH2 is ubiquitous in non-epithelial tissues and acts as a biomarker of cells undergoing EMT, as well as mediating cell survival and migration [56]. MMP2 participates in the degradation of collagen and basement membrane and is associated with cell migration [57]. Variants of the matricellular glycoprotein SPARC play a crucial role in the regulation of EMT induced by TGF- $\beta{}$ [58]. The EMT-related genes KRT7, FGF1, LOXL2, SHROOM3 and TGFBI encode for keratin 7, a transcription factor, a growth factor, lysyl oxidase, a shroom family member and TGF- $\beta{}$ induced protein, respectively. They have been shown to facilitate the EMT program in various cell types [59, 60, 61, 62, 63]. The two remaining EMT-related genes, CLDN3 and ELF3, act as negative regulators of EMT. CLDN3 is a major structural component of tight junctions in epithelial cells [64], while ELF3 inhibits EMT by reducing the expression of ZEB1 [65]. Based on the observation that CLDN3 and ELF3 mRNA levels were up-regulated after TGF- $\beta{}$ treatment, we speculate they may undergo further regulation at the translational level.

In previous studies, TGF- $\beta{}$ treatment was reported to significantly up-regulate the m ${}^{6}$ A levels of EMT-associated transcription factors such as Snail, ZEB and JUNB in lung cancer cells and HeLa cells [23, 38]. However, in the present study we found that m ${}^{6}$ A levels of EMT genes were reduced in RPE cells following TGF- $\beta{}$ 2-induced EMT, suggesting m ${}^{6}$ A methylation could have opposite roles during EMT in different cell types. Although we found potential target genes for m ${}^{6}$ A modification, further in-depth and systematic studies are required to fully elucidate the regulatory mechanism of m ${}^{6}$ A. Additional questions that need to be addressed are whether the regulation of m ${}^{6}$ A on EMT-related genes occurs pre- or post-transcription, and whether m ${}^{6}$ A is also involved in regulation at the translational level. Follow up studies should therefore use molecular techniques to verify the role of m ${}^{6}$ A and to identify m ${}^{6}$ A readers that play a key role in regulating EMT target genes. The relationship between m ${}^{6}$ A and EMT also needs further confirmation in experimental animal models and in clinical specimens before treatments for EMT can be developed.

The EMT process of RPE cells plays an important role in the pathogenesis of retinal diseases, such as PVR and age-related macular degeneration (AMD) [66]. However, the pathogenesis of EMT is still under investigation, therefore the non-surgical treatment on retinal diseases is limited at present. In this study, by establishing the EMT model of RPE cells, we found that the transcripts of 12 EMT-related genes were hypomethylated, and the total m ${}^{6}$ A level was down-regulated during EMT. In addition, the m ${}^{6}$ A methyltransferases METTL3 was shown to attenuate EMT [26]. It is speculated that hypermethylation could inhibit EMT in RPE cells. Therefore, it could be desirable to find a way to increase the intracellular m ${}^{6}$ A level by some methods to inhibit EMT and treat retinal diseases. For example, small molecule inhibitors of m ${}^{6}$ A demethylase FTO (FB23 and FB23-2), which have been shown to play an inhibitory role in the development of leukemia by inhibiting FTO activity could be considered [67]. Furthermore, two other inhibitors of FTO, CS1 and CS2, were found to have more potent anti-tumor efficacy and fewer side effects [68]. Whether these FTO inhibitors can affect the EMT process by acting on the transcript of EMT-related genes, and whether they inhibit the development of retinal diseases remains to be revealed.

5. Conclusions

We identified differentially methylated genes and differentially expressed genes using MeRIP-seq and RNA-seq, respectively, following treatment of RPE cells with TGF- $\beta{}$ 2, an inducer of EMT. We also provide preliminary results on the m ${}^{6}$ A modification profile in RPE cells after TGF- $\beta{}$ 2 treatment. m ${}^{6}$ A modification was found to be widespread and the m ${}^{6}$ A motif was conserved. Further association analysis identified 12 EMT-related genes that were altered at both the methylation level and the mRNA level, suggesting they may be targets of m ${}^{6}$ A methylation that affects EMT.

Availability of Data and Materials

The datasets generated and/or analyzed during the current study are available in the NCBI SRA repository (https://www.ncbi.nlm.nih.gov/sra/) with the SRA accession number: SRP422412.

Author Contributions

XiaL designed the research study and provided advice on supervision. XZ and XueL performed the research. XZ, LL and YZ analyzed the data. FW, RY and MY validated the data. XZ wrote the manuscript. XZ and XiaL revised the manuscript. All authors contributed to editorial changes in the manuscript. All authors read and approved the final manuscript.

Ethics Approval and Consent to Participate

Not applicable.

Acknowledgment

Sequencing was performed by Wuhan Seqhealth Technology. The Omicshare platform was used to analyze GO enrichment and KEGG enrichment of differential genes.

Funding

This research was funded by National Natural Science Foundation of China, grant number 81770952 and 81873681; Basic Research Program of Henan Eye Hospital, grant number 20JCQN005; Natural Science Foundation of Henan Province, grant number 222300420194; Henan Provincial Medical Science and Technology Research Joint Co-construction Project, grant number LHGJ20210075.

Conflict of Interest

The authors declare no conflict of interest.

Supplementary Material

Supplementary material.zip

References

[1]

Zhou M, Geathers JS, Grillo SL, Weber SR, Wang W, Zhao Y, et al. Role of Epithelial-Mesenchymal Transition in Retinal Pigment Epithelium Dysfunction. Frontiers in Cell and Developmental Biology. 2020; 8: 501.