Comprehensive Analysis of the Transcriptome-Wide m6A Methylome in Pterygium by MeRIP Sequencing

Highlights

- Both level of m6A and expression of METTL3 were decreased in pterygium.

- 2,145 Upregulated and 804 downregulated m6A peaks were found in pterygium.

- M6A-targeted genes were enriched in Hippo and Notch signaling pathway.

- 72 Hypermethylated and 15 hypomethylated m6A peaks were differently in mRNA.

- Five hub genes were most related to the development of pterygium.

Introduction

Post-transcriptional modification plays an important regulatory role in various physiological processes and in disease progression (He et al., 2019). To date, more than 100 RNA modifications have been confirmed (Boccaletto et al., 2018). Among them, N6-methyladenosine (m6A), the result of methylation at the sixth nitrogen of adenine, is the most abundant mRNA modification. m6A has an RRACH consensus motif (wherein each R represents A or G, and H represents A, C, or U) (Zhang et al., 2020a), and its related proteins include “Writer” (METTL3, METTL14, and WTAP), “Eraser” (FTO and ALKBH5), and “Reader” (YTHDF1, 2, and 3; YTHDC1-2; IGF2BPs) (Hu et al., 2020). During the transcription of DNA into RNA, adenosine is modified by the action of the writer proteins METTL3, METTL14, and WTAP at the sixth N position (Schöller et al., 2018), and these m6A methylated bases can be recognized by the reader proteins recognition by the YTH Family Proteins (including YTHDF1, 2, 3; YTHDC1, 2), which are involved in downstream translation (Liu et al., 2020), mRNA translation (Wang et al., 2015) and degradation (Lee et al., 2020), and accelerated mRNA exit from the nucleus (Roundtree et al., 2017). In addition, bases that have undergone m6A modification are de-methylated in the presence of two enzymes, FTO and ALKBH5, making the m6A reaction reversible (Chen et al., 2019). Recent studies have shown that m6A is closely involved in many biological functions, such as the regulation of the Epithelial-to-mesenchymal transition (EMT) process, the DNA damage response, gene translation, autophagy, adipogenesis, and the determination of stem cell fate (Zhou et al., 2020).

Pterygium is a triangular fibrovascular tissue that grows from the bulbar conjunctiva of the palpebral fissure to the cornea, and it often occurs on the side of the nose (Chen et al., 2014). According to the latest epidemiological survey, the prevalence of pterygium in Asia is 7% (Ang et al., 2012; Chen et al., 2015; Khanna et al., 2020), while in Africa and South America, it is higher (Nemesure et al., 2008; Fernandes et al., 2020). Therefore, it is necessary to explore the specific mechanism of pterygium. At present, no such research has been published. The literature supports the important role of EMT and DNA damage in the pterygium mechanism (Cimpean et al., 2013; Meshkani et al., 2019). m6A is involved in the regulation of the EMT process and DNA damage response and other biological functions, which suggests that it may be involved in the mechanism of pterygium development. Therefore, we investigated whether there is a connection between m6A and the occurrence of pterygium.

To this end, we used high-throughput sequencing to perform a genome-wide analysis of m6A-tagged transcripts. Through m6A-modified RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq), we identified differentially methylated peaks and differentially expressed genes in mRNA. Further analyses determined that five genes might be associated with pterygium.

Materials and Methods

Sample Collection

The study was approved by the Ethics Committee of Yangpu Hospital, Tongji University and was conducted in accordance with the principles of the Declaration of Helsinki. The severity of pterygium is graded by YC and YJ based on slit lamp examination and slit lamp photography.

The severity of pterygium is graded according to the degree of extension and is divided into four grades as follows (Chien et al., 2013):

Grade 1: The pterygium is still involved in the conjunctiva.

Grade 2: Pterygium involvement includes only the corneal limbus.

Grade 3: Pterygium involvement between the corneal limbus and the pupillary limbus.

Grade 4: Pterygium involvement extends centrally to the pupillary margin.

A total of 24 patients with pterygium diagnosed as grade III or above were enrolled and they all provided informed consent. People with ocular diseases other than pterygium and those who had undergone ophthalmic surgery were excluded from the study and all patients underwent pterygium excision combined with autologous conjunctival transplantation. Pterygium tissue was obtained by surgical excision, and a small rectangular piece of normal conjunctival tissue was taken from the autograft at the contralateral corneal limbus of the same eye as a control. The collected sample is transferred to a 1.5 mL tube and stored at −80°C for subsequent analysis. Detailed steps are shown in Figure 1A. Clinical information on gender, age, pterygium size, and whether UV exposure were included in the 24 pterygium patients is summarized in Supplementary Table 1.

FIGURE 1

Figure 1. Workflow of the study. (A) Workflow of the experimental process. (B) Workflow of the data analysis.

Western Blot Analysis

Six pairs of pterygium and conjunctival tissues were lysed for 1 min in cold RIPA buffer [P0013C, Beyotime, the main components are 50 mM Tris (pH 7.4), 150 mM NaCl, 1% NP-40, 0.5% sodium deoxycholate, 0.1% SDS] supplemented with 1 mM PMSF (ST505, Beyotime) using the GeneReady Standard (BSH-2, Life Real). The cells were scraped off and collected into 1.5 mL EP tubes, centrifuged at 12,000 rpm for 15 min and the supernatant was removed for protein quantification using the BCA Protein Concentration Assay Kit (P0012, Beyotime). After denaturation, 30 mg total protein was separated on 10% SDS-PAGE gels and transferred to PVDF membranes. The gels were blocked for 1 h using 5% bovine serum albumin and then incubated overnight at 4°C in the following primary antibodies: anti-METTL3 (Abcam, 1:1000, cat. No ab195352) (Pan et al., 2021), anti-METTL14 (Abcam, 1:1000, cat. No ab220030) (Yi et al., 2020), anti-WTAP (Abcam, 1:1000, cat. No ab195380) (Chen S. et al., 2020), anti-FTO (Abcam, 1:5000, cat. No ab126605) (Baquero-Perez et al., 2019), anti-ALKBH5 (Abcam, 1:1000, cat. No ab195377) (Zhu et al., 2020) [with GAPDH (Proteintech, 1:2000, cat. No 60004-1-Ig) (Lin et al., 2019) as a loading control]. After three washes with TBS-Tween, the proteins were incubated with the corresponding secondary antibodies conjugated to horseradish peroxidase and then detected via enhanced chemiluminescence (ECL) using the iBright CL1000 device (Invitrogen).

RNA Isolation and Quantitative Real-Time PCR

Using TRIzol reagent (Invitrogen), total RNA was isolated from 18 pairs of pterygium and normal conjunctival tissues (discovery phase n = 12; validation phase n = 6) from 18 pterygium patients. Twelve pairs of pterygium and normal conjunctival tissue RNA from the discovery phase were used independently for m6A and m6A modification-related regulators [include “Writer” (METTL3, METTL14, and WTAP), “Eraser” (FTO and ALKBH5), and “Reader” (YTHDF1)] expression level detection. For the remainder of 12 pairs of pterygium and normal conjunctival tissue RNA, every four equal independent RNA samples were pooled into two types of samples for MeRIP and RNA-seq analysis. Six pairs of pterygium and normal conjunctiva tissue RNA from the validation phase were used independently for hub gene validation. After quantifying the samples using a NanoDrop ND-1000 instrument, total RNA was reversed to cDNA using the Qiagen RNeasy kit according to manufacturer’s protocol (TaKaRa, 6110). Real-time qRT-PCR was performed on an ABI Prism 7500 Sequence Detection system using SYBR green reagent (Applied Biosystems); GAPDH was used as the endogenous control gene, the primers are listed in Supplementary Table 2. The relative gene expression levels were calculated using the 2-ΔΔCt method.

Quantification of m6A in Total RNA

Total RNA from 18 pairs of pterygium and normal conjunctiva tissues from discovery phase was used for quantification of m6A modifications using an m6A RNA Methylation Quantification Kit (P-9005, EpiGentek). First, we prepared a negative control and also a positive control using six different concentrations of m6A ranging from 0.01 to 0.5 ng/μL, as recommended by the manufacturer. Then 2 μL positive control, 2 μL negative control, or 200 ng total RNA were bound in the 96-well plates. After removing the binding solution from each well, the samples were washed three times. The wells were incubated for 1 h with the addition of anti-m6a antibody and washed four times. Then, add the detection antibody to the wells, incubate at room temperature for 30 min and wash four times. Add the enhancer solution, incubate at room temperature for 30 min and then wash five times. Finally, add the developer solution and incubate at room temperature until the developer solution turns blue and then add the stop solution, the samples were analyzed on a SYNERGY 2 (BioTek) microplate reader at 450 nm; the absolute amount of m6A was quantified and the percentage of m6A in RNA was established as%m6A in total RNA.

Immunohistochemistry Staining

Six pairs of pterygium and normal conjunctiva tissues from discovery phase were embedded in formalin. Each tissue was cut to 4-μm thick and mounted on a glass slide. Dewaxing sections were performed as previously described (Kawaguchi et al., 2002). Endogenous peroxidase activity was inhibited and blocked with 5% bovine serum albumin for 30 min at 37°C. The slices were incubated in anti-METTL3 (Abcam, 1:400, ab195352) (Zhang et al., 2020b) overnight at 4°C, washed three times with PBS for 5 min, and then incubated with secondary anti-horseradish peroxide at 37°C for 30 min. After three more washes with PBS, the slices were visualized in diaminobenzidine chromogenic solution. Microscopic images were obtained by light microscopy (Carl Zeiss, Oberkochen, Germany).

MeRIP Sequencing, RNA Sequencing, and Bioinformatics Analyses

m6A-modified RNA immunoprecipitation sequencing and RNA-seq were performed by Guangzhou Epibiotek Co., Ltd. Briefly, the Zymo RNA Clean and Concentrator-5 Kit was used to purify RNA from total RNA samples; then these were fragmented into 200 nucleotide lengths. A portion of the purified RNA was used for RNA-seq as input. The rest of the purified RNA was enriched via immunoprecipitation with an anti-m6A antibody (Abcam, cat. No. ab151230) in IP buffer (300 mM NaCl, 0.2% Igepal CA-630, 20 mM Tris-HCl, PH 7.4) overnight at 4°C. Then, the m6A-Ab mixture was incubated with protein-g magnetic beads (Thermo Fisher, cat. No. 88847) at 4°C for 2 h for immunoprecipitation. After washing and elution, IP RNA is obtained. Then, first-strand cDNA about m6A antibody-enriched RNAs and input RNAs were synthesized by SMART. The library of ultramicro-RNA methylated m6A and input RNAs were obtained by PCR amplification and library purification. Bioptic Qsep100 analyzer is used to conduct quality control over the library. Finally, the samples were sequenced using Illumina NovaSeq^TM 6000 in PE150 sequencing mode.

After obtaining clean data from Illumina NovaSeq^TM 6000, quality control was performed using Cutadapt (v2.5), including trimming aptamers and filtering sequences, and then the remaining reads were aligned to human pooled genome GRCh38 under parameters: “rna_strandness RF” using the Hisat2 aligner (v2.1.0). Mapped reads from the IP and input libraries were provided to the exomePeak R package (v2.13.2) to identify m6A peaks under the parameters: “PEAK_CUTOFF_PVALUE = 0.05, PEAK_CUTOFF_FDR = NA, FRAGMENT_LENGTH = 200.” Use the exomePeak R package to identify differential m6A peaks under the parameter: “PEAK_CUTOFF_PVALUE = 0.05, PEAK_CUTOFF_FDR = NA, FRAGMENT_LENGTH = 200.” Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed using clusterProfile R package (v3.6.0). m6A-RNA-related genomic features were visualized using Guitar R package (v1.16.0). Identified m6A peaks which P value <0.05 were chosen for the de novo motif analysis using homer (v4.10.4) under parameters: “−len 6 −rna.”

The expression of gene transcripts was quantified using featureCount (v1.6.3), and differentially expressed genes were identified by DESeq2 R package (1.18.1). Genes were considered as differentially expressed if gene expression is FC ≥ 2 and with a FDR ≤ 0.05.

All codes are available on GitHub (GitHub, Inc., San Francisco, CA, United States) at https://github.com/zxiahong/titre-trajectory-aje. To further explore the essential role of m6A mRNA modification in the development of pterygium, the identified m6A peaks and Differentially expressed genes (DEGs) were subjected to network analyses¹ as previous described (Schuierer et al., 2010) (Figure 1B).

Public Databases and Analysis

Expression profiling data (GSE83627, GSE51995, and GSE2513) were downloaded from the NCBI Gene Expression Omnibus (GEO)² database. The expression levels of genes were analyzed using R software (version 4.0.2).³ The Student’s t-test was used to compare the relative gene expression between pterygium and normal conjunctiva tissues. GO and KEGG pathway analyses were performed based on significantly methylated genes and/or significantly expressed genes [|log2Fold Change (FC)| > 1 and P-value < 0.05] (Huang da et al., 2009).

Results

Levels of m6A and METTL3 Are Downregulated in Pterygium

First, we quantified m6A levels in pterygium and normal conjunctival tissues. Interestingly, as shown in Figure 2A, the level of m6A was significantly lower in pterygium than in normal conjunctival tissues (−0.01117 ± 0.00289, Fold-Change = 0.56892, P-value = 0.0008). Next, using qRT-PCR, we examined the mRNA expression levels of five crucial enzymes needed for m6A modification: METTL3, METTL14, WTAP, FTO, and ALKBH5. We found that the mRNA levels of METTL3 were extraordinarily lower in pterygium versus conjunctive (−0.7155 ± 0.1434, Fold-Change = 0.43728, P-value < 0.0001, Figure 2B). A similar decrease in METTL3 was also found in Western blotting analyses (−0.2752 ± 0.08527, Fold-Change = 0.68424, P-value = 0.0091, Figures 2C,D) and immunohistochemistry staining (especially in epithelium, Figure 2E). Also, decreasing expression of FTO in mRNA and protein were also found qRT-PCR and western blotting. However, this trend violated the change of m6A modification. Taken together, these results indicate that METTL3 is the key m6A modification regulator in pterygium.

FIGURE 2

Figure 2. m6A and METTL3 are downregulated in pterygium. (A) Quantification of m6A in total RNA in pterygium and conjunctiva tissues. N = 12 per group. (B) mRNA expression level of the major enzymes of m6A in pterygium and conjunctiva tissues: METTL3, METTL14, WTAP, FTO, ALKBH5, YTHDF1. N = 12 per group. (C,D) Protein expression levels of these enzymes determined via Western blotting (C,D) the results of corresponding densitometric analyses. N = 12 per group. (E) Immunohistochemistry staining of METTL3 in pterygium and conjunctiva tissues. Scale bar: 25 μm. Results are expressed as the mean ± SEM (NS indicates not significant, **P < 0.01, ***P < 0.001, compared to the conjunctiva group).

Overview of the m6A Methylation Map in Pterygium

After pool four independent RNA from pterygium and normal conjunctiva to get two types of samples for the MeRIP sequence, we performed MeRIP-seq analysis on three pairs of pterygium and normal conjunctival samples. As shown in Figure 3A, compared to conjunctiva, pterygium had 1155 significantly upregulated m6A peaks and 375 downregulated m6A peaks (fold changes ≥ 2). Table 1 list the top 20 peaks, and Figure 3B shows the distribution of m6A signals around mRNA and lncRNA were comparable in the two classes of tissue samples. The tissues showed different patterns of peaks with a relative increase at the start codon region (4.4 vs. 3.6% for pterygium and conjunctiva, respectively) and in the 3′ untranslated region (3′UTR, 27.9 vs. 26.6%), and a relative decrease in the coding sequence (CDS, 43.8 vs. 44.1%) and at the stop codon (22 vs. 23.7%) (Figures 3C,D). Also, the classic GGAC motif (Fu et al., 2014; Liu et al., 2014) and the top five m6A motifs in pterygium (Figure 3E) and conjunctiva (Supplementary Figure 1A) were observed. Analyses of the distribution of m6A peaks indicated that most mRNAs and genes had m6A peaks (mRNAs with 458 downregulated peaks and 1,300 upregulated peaks; genes with 487 downregulated peaks and 1,368 upregulated peaks, Figures 3F,G). Variant peaks were found in all chromosomes, especially chr1, chr3, chr11, chr12, and chr16 (Figure 3H).

FIGURE 3

Figure 3. Overview of altered m6A-modified transcripts in pterygium. (A) Volcano plots showing significantly different m6A peaks. (B) Average m6A peaks in conjunctiva and pterygium. (C,D) Pie charts showing the distribution of m6A peaks in conjunctiva (C) and pterygium (D). (E) Top five m6A motifs from the altered m6A peaks in pterygium. (F) Distribution of altered m6A peaks per mRNA. (G) Distribution of altered m6A peaks per gene. (H) Distribution of altered m6A peaks in human chromosomes. Fold change ≥2 and P < 0.05.

TABLE 1

Table 1. Top 20 altered m6A peaks of m6A-modified mRNA in the pterygium vs conjunctiva.

Differentially Methylated mRNAs Participate in Important Signaling Pathways

To study the biological significance of m6A modification in pterygium, GO and KEGG pathway analyses of differentially methylated mRNAs were conducted. For GO analyses, we considered biological processes (BP), cellular components (CC), and molecular functions (MF). The main results are given in Supplementary Figure 1B. Using Metacore system, we tested dys-methylated lncRNA which were significantly enriched in dosage compensation by inactivation of X chromosome, dosage compensation and nuclear speck organization (Figure 4A). In addition, the top scored network analyses of these were found to be relative with nervous system development (76.2%), excitatory chemical synaptic transmission (16.7%), intracellular signal transduction (61.9%), response to abiotic stimulus (57.1%), system development (85.7%) (Figure 4B and Supplementary Table 3). In the pathway map analysis of mRNA, the hypermethylated peaks were strongly associated with negative regulation of WNT/Beta-catenin signaling in the nucleus. Besides, the most relevant biological pathways of hypomethylated peaks were WNT/Beta-catenin signaling pathway (Signalosome) and Canonical Notch signaling pathway in colorectal cancer (Figure 4B). In addition, GO process results showed that the hypermethylated peaks were significantly enriched in signal transduction (ESR1-nuclear pathway) and cell adhesion (Cadherins), while the hypomethylated peaks were enriched in regulation of signal transduction and regulation of cellular metabolic process (Figure 4C). Furthermore, the top three scored network analyses of hypermethylated peaks about the development of pterygium were found to be involved in nervous system development (76.2%), excitatory chemical synaptic transmission (16.7%), intracellular signal transduction (61.9%), response to abiotic stimulus (57.1%), system development (85.7%) (Figure 4D and Supplementary Table 3). The top three scored networks of hypomethylated peaks were enriched in canonical Wnt signaling pathway (52.0%), cell-cell signaling by wnt (58.0%), Wnt signaling pathway (58.0%), cell surface receptor signaling pathway involved in cell-cell signaling (58.0%), regulation of multicellular organismal development (88.0%) (Figure 4E and Supplementary Table 3). Additionally, in KEGG pathway analyses, the hippo signaling pathway, endocytosis, and cancer-related pathways were significantly correlated with genes that showed upregulated m6A peaks in pterygium (Figure 5A). Genes with downregulated m6A peaks were significantly correlated with the Notch signaling pathway and Mucin type O-Glycan biosynthesis (Figure 5B).

FIGURE 4

Figure 4. Systematical analyses of differentially methylated mRNA in pterygium. (A) Gene Ontology (GO) processes analysis of differentially methylated lncRNA. (B) Top three scored networks analysis of differentially methylated lncRNA. (C) Pathway maps analysis of differentially methylated peaks. Left panel represents hypermethylated genes and right panel represents hypomethylated peaks. (D) GO process analysis of differentially methylated peaks. Left panel represents hypermethylated peaks and right panel represents hypomethylated peaks. (E) Top three scored networks analysis of differentially methylated peaks. Left panel represents hypermethylated genes and right panel represents hypomethylated peaks.

FIGURE 5

Figure 5. KEGG enrichment analyses of differentially methylated mRNA in pterygium. Top 10 KEGG pathways of genes with (A) upregulated m6A peaks and (B) downregulated m6A peaks. Squares in left semicircle refer to the upregulated or downregulated mRNAs, the squares in right semicircle refer to the pathways. *P < 0.05, **P < 0.01, ***P < 0.001, genes compared to pathways.

Overview of Transcriptome Profiles and Conjoint Analyses of MeRIP-Seq and RNA-Seq Data

First, we tested the transcriptome profiles of altered genes in three pairs of pterygium and normal conjunctival samples using RNA-seq. Compared to normal conjunctival samples, pterygium samples had 394 significantly upregulated genes and 679 significantly downregulated genes (Figures 6A,B). The top 20 differentially expressed genes are listed in Table 2. After the systematically analyzed by Metacore, we found the different expression of lncRNA were significantly enriched in development (YAP/TAZ mediated co regulation of transcription) in pathway maps and serval types of cells differentiation in GO processes (Figure 6C). in addition, the top three scored network analyses of these were found to be strongly associated with regulation of transcription by RNA polymerase II (62.8%), regulation of biosynthetic process (79.1%), regulation of macromolecule biosynthetic process (76.7%), response to lipid (72.3%), positive regulation of nitrogen compound metabolic process (83.0%), cellular response to lipid (53.2%), response to organic cyclic compound (63.8%), positive regulation of cellular metabolic process (83.0%) (Figure 6D and Supplementary Table 4). Additionally, cell cycle (Start of DNA replication in early S phase) in the pathway maps analysis, muscle contraction in process networks and phenol-containing compound metabolic process in GO processed were the top three enrichment in different expression of mRNA (Figure 6E).

FIGURE 6

Figure 6. Results of significant RNA of pterygium in RNA-seq data. (A) Volcano plots and (B) heatmap plots showing differentially expressed genes in pterygium vs. conjunctiva. Fold change ≥2 and P < 0.05. (C) GO processes and pathway maps analysis of significant lncRNA. (D) Top three scored networks analysis of significant lncRNA. (E) Pathway maps, process networks and GO processes analysis of significant mRNA. Left panel represents pathway maps, middle panel represents process networks and right panel represents GO processes.

TABLE 2

Table 2. Top 20 differentially expressed genes in pterygium vs conjunctiva after RNA-seq.

Furthermore, conjoint analyses of MeRIP-seq and RNA-seq data were performed (Figure 7A). We identified 72 hypermethylated m6A peaks in mRNAs that were significantly upregulated (20) or downregulated (52), while 15 hypomethylated m6A peaks in mRNAs were significantly upregulated (10) or downregulated (5). Next, 87 genes that showed significant changes in both m6A modification and RNA expression levels were subjected to GO, pathway analyses and networks analyses. The GO analyses of processes associated with the four gene sets showed detail in Figure 7B, which apparently linked numerous functional processed and pathways. Interestingly, we found Cytoskeleton remodeling (Keratin filaments) and Statin action on the PI3K/Akt pathway in COPD were enriched in pathway maps, and Cell adhesion (Integrin-mediated cell-matrix adhesion) and Signal transduction (Androgen receptor signaling cross-talk) were enriched in process networks (Figure 7C). In addition, the top 20 GO and KEGG pathways that were up- or downregulated in both m6A modification and RNA expression genes are shown in Supplementary Figures 1C–E. Moreover, the top three scored network analyses indicated that m6A-modified genes of pterygium were significantly enriched in response to growth factor (25.8%), regulation of gene expression (63.4%), regulation of RNA metabolic process (54.8%), cellular response to growth factor stimulus (23.7%), cellular response to peptide hormone stimulus (19.4%) (Figure 7D and Supplementary Table 5). After analyzed the key transcription factors and target genes in m6A-modified genes, we showed the top three transcription factors (GATA-1, FOXP3, and TAL1) and the m6A-modified genes which most enriched in regulation of cell shape (14.7%), regulation of cell morphogenesis (20.6%), plasma membrane repair (5.9%), regulation of extrinsic apoptotic signaling pathway (11.8%), response to wounding (20.6%), regulation of gene expression (60.0%), negative regulation of gene expression (36.7%), regulation of metabolic process (70.0%), positive regulation of transforming growth factor beta1 production (6.2%) (Figure 7E). These observations indicate that genes with m6A modification may play important role in the development of pterygium.

FIGURE 7

Figure 7. Systematical analyses of MeRIP-seq and RNA-seq data in pterygium. (A) Significantly m6A-modified RNA identified after conjoint analyses of MeRIP-seq and RNA-seq data. (B) Bubble diagram of GO biological process categories enriched for DEGs with m6A hyper- or hypo-methylated. (C) Pathway maps analysis (left) and Process networks analysis (right) of differentially m6A-modifed genes. (D) Top three scored networks analysis of m6A-modified genes. (E) Key transcription factors and target genes network analysis of m6A-modifed genes. DEGs, differentially expressed genes.

Validating the 87 Genes That Showed Significant Changes in Both m6A Modification and RNA Expression Levels Against the GEO Datasets

To further validate the most relevant gene correlations with pterygium, we examined 87 genes that showed significant changes in both m6A modification and RNA expression levels in GEO datasets including GSE83627, GSE51995 (Hou et al., 2014), and GSE2513 (Wong et al., 2006) (Figures 8A–C). Finally, five hub genes showed consistent difference in three GEO datasets and our study. Their expression levels in each tissue were validated by qRT-PCR (Figure 8D). The expressions of DSP, MXRA5, TMEM43, and OLFML2A were increased, similar to our sequence data, while the decrease in expression of ARHGAP35 was not found. The patterns of m6A modification in these five genes are visualized Supplementary Figure 2.

FIGURE 8

Figure 8. Validation gene expression in GEO datasets. (A,B) The expression of genes from GEO datasets (GSE2513, GSE51995, and GSE83627) were normalized to remove batch effects. (C) Wayne chart of the three GEO datasets. (D) mRNA expression level of five hub genes (DSP, MXRA5, TMEM43, OLFML2A, and ARHGAP35) determined via qPCR. N = 6 per group. Results are expressed as the mean ± SEM (NS indicates not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001 compared to the conjunctiva group).

Discussion

Epigenetic modification is associated with many diseases (Chen et al., 2017; Rezasoltani et al., 2017). RNA modifications, especially m6A, were reported as a novel territory in epigenetic modification to be much associated with various diseases nowadays, such as tumorigenesis, cardiovascular diseases, inflammation, and Type II diabetes (He et al., 2013; Gustavsson et al., 2014; Rezasoltani et al., 2017; Dorn et al., 2019; Zhang L. et al., 2020). However, the role of m6A modification in pterygium development has not yet been elucidated.

In this study, we provide the first data analysis in pterygium. We found significantly downregulated m6A levels in pterygium compared to conjunctival tissues, especially by METTL3, which may be the pivotal regulator in pterygium. Interestingly, METTL3 expression has been found to play a role in the regulation of cell growth-related pathways as well as cell death-related pathways in cancer (Barbieri et al., 2017; Chen et al., 2018). And, previous studies showed the cell growth-related pathways and cell death-related pathways were important for the development of pterygium.

Further analyses identified significantly up- and downregulated m6A peaks, the potential functions of altered m6A-modified transcripts, and synchronously differentially expressed genes that showed both m6A modification and RNA expression. Overall, the analysis shows an altered landscape of m6A modification is pterygium.

In our MeRIP-seq analysis, we found 458 downregulated m6A peaks and 1,300 upregulated m6A peaks on mRNA. We conducted GO and KEGG analysis to study the biological functions of these methylated mRNAs. We found dys-methylated lncRNA were significantly enriched in dosage compensation by inactivation of X chromosome. Dosage compensation of the X chromosome is the inactivation of one of the two X chromosomes by heterochromatization so that only one chromosome is active, resulting in a balanced dose of gene products on the X chromosome in males and females, although the number of X chromosomes is not the same (Jordan et al., 2019). Dosage compensation of the X chromosome is associated with epigenetics as well as gene expression patterns (Brockdorff and Turner, 2015), and long non-coding RNAs are key players in the regulation of the X chromosome (Sahakyan et al., 2018). Therefore, aberrantly methylated lncRNAs may cause alterations in the epigenetics of pterygium as well as gene expression patterns by regulating the inactivation of the X chromosome.

In KEGG pathway analyses, we found that upregulated m6A modification genes are related to the Hippo signaling pathway, which was first discovered in Drosophila and is a highly conserved kinase cascade composed of tumor suppressor Hippo (MST and Lats in mammals) and the downstream effector molecule Yki (YAP and TAZ in mammals). The latter is a transcriptional co-activator of target genes involved in cell proliferation and survival (Pan, 2010). Many studies have shown that the Hippo signal transduction pathway can regulate the transcriptional co-activation molecules YAP and TAZ through negative regulation, thereby regulating organ size, cancer occurrence, tissue regeneration, and stem cell function (Mo et al., 2014; Yu et al., 2015). YAP and TAZ are transcriptional co-activators that shuttle between the cytoplasm and the nucleus. They recognize homologous cis-regulatory elements by interacting with other transcription factors such as TEA domain family members, which are involved in cell proliferation, tissue regeneration, and stem cell maintenance (Piccolo et al., 2014; Chen M. et al., 2020). The continuous activation of this system can promote abnormal cell proliferation (Koo and Guan, 2018). YAP can regulate the cell cycle by controlling the cyclin to enter the S phase or inducing other proto-carcinoma transcription factors such as c-Myc (Zanconato et al., 2015). The occurrence of pterygium is related to cell proliferation, cell cycle changes, and a lack of limbal stem cells. The Hippo signal transduction pathway is most relevant to pterygium, suggesting that m6A modification may promote the occurrence and development of pterygium by regulating the Hippo signal transduction pathway to affect cell proliferation, cell cycle changes, and loss of limbal stem cell function. We also found that downregulated m6A modification genes were related to the Notch signaling pathway, which is involved in regulating various cellular processes of normal development, such as cell proliferation, differentiation, the epithelial–mesenchymal transition, angiogenesis, and many biological processes that promote the occurrence and development of cancer (Aster et al., 2017). The Notch pathway is mainly composed of four parts: Notch receptor (Notch1–4), Notch ligand (Delta-like [Dll] 1, Dll3, Dll4, Jagged-1, and Jagged-2), CSL-DNA binding protein, and downstream target genes (Meurette and Mehlen, 2018). Zhang et al. (2017) showed that Notch4 is significantly highly expressed in human prostate cancer (PCa), and can promote the progression and metastasis of PCa through cell growth, apoptosis, migration, invasion, and the EMT. Another study (De Francesco et al., 2018) showed that Notch-related signaling pathways can promote the occurrence and development of breast cancer by regulating the EMT. However, the precise functions of this pathway in pterygium remain to be discovered and the m6A-modified genes related to this pathway in pterygium must be further elucidated.

Combining MeRIP-seq and RNA-seq data, we identified 87 differentially methylated m6A peaks and synchronously differentially expressed in mRNA. The expression of these genes may be regulated by m6A modification. Using GEO databases, we also identified five hub genes (DSP, MXRA5, ARHGAP35, TMEM43, and OLFML2A) that were most correlated with the development of pterygium. DSP encodes the desmoplakin protein, and is a critical component of desmosome structures in cardiac muscle and epidermal cells (Arnemann et al., 1991; Bornslaeger et al., 1996). Other studies have reported different roles of DSP in cancer (Papagerakis et al., 2009; Yang et al., 2012), but it is generally clear that it is a key factor in the EMT. BMP6 significantly upregulates the expression of desmoplakin and desmoglein in human corneal limbal epithelial cells, consistent with its pro-differentiation role (Tiwari et al., 2020). Matrix-remodeling associated 5 protein, which is encoded by MXRA5, is a novel biomarker in colorectal cancer (Wang et al., 2013) and non-small-cell lung carcinoma (Xiong et al., 2012). The high expression of MXRA5 is strongly associated with the remodeling of the extracellular matrix (ECM). MXRA5 expression may also be involved in TGFβ activation, which plays an important role in several developing connective tissues (Robertson et al., 2015; Bertoni et al., 2016). Interestingly, the high expression of MXRA5 in our study suggests that it may remodel the ECM by activating TGFβ in pterygium. Transmembrane protein 43 (TMEM43) plays an important role in maintaining the structure of the nuclear envelope (Wiemann et al., 2001), and mutations in TMEM43 are strongly associated with arrhythmogenic right ventricular cardiomyopathy (Baskin et al., 2013). The high expression of TMEM43 in pterygium, as in cancer, may act as a critical component in the EGFR signaling network to control cell survival, migration, and invasion (Jiang et al., 2017). The ARHGAP35 gene encodes glucocorticoid receptor DNA binding factor-1 (also known as p190RHOGAP), which belongs to the large G protein/GTPase activation protein family. It can regulate rho-dependent signal transduction processes, such as cell adhesion, migration, invasion, and cytokinesis (Héraud et al., 2019). Its downregulation promotes the transformation and growth of epithelial tumor cells (Sun et al., 2014). In pterygium, it increases cell invasion and migration, thereby inducing or promoting the development of pterygium. OLFML2A is a novel regulator of transforming growth factor-β-induced smooth muscle differentiation and a novel Carcinogenic factor in liver hepatocellular carcinoma (Shi et al., 2014; Bai et al., 2020). Its role in pterygium has not been illustrated and further research is required.

The m6A modification plays an important role in the regulation of gene expression, and an abnormal level of m6A modification may be related to human diseases or cancer. It has been found that enzymes related to m6A modification may affect the UV-induced DNA damage response (Lin et al., 2016), tumor formation (Cui et al., 2017), and cell development and differentiation (Li et al., 2017). The decrease in m6A modification level caused by the downregulation of METTL3 may be an important reason for the development of pterygium, but further research is needed to clarify the exact mechanism. Our research has identified many genes that are hypermethylated or hypomethylated and differentially expressed simultaneously. These genes may play an important role in the development of pterygium. Our findings suggest that adjusting the level of m6A modification may become a new treatment strategy for pterygium in the future.

Conclusion

High-throughput sequencing produced a broad map of the m6A transcriptome in pterygium. We conducted joint analyses of MeRIP-seq and RNA-seq data and identified many differentially expresses m6A methylation peaks and differentially expressed genes, indicating that there is a potential connection between m6A methylation and gene expression. In addition, the genes we selected are related to the occurrence of pterygium.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI GEO; GSE167933.

Ethics Statement

This study was performed in accordance with the principles of the Declaration of Helsinki, and the study protocol was reviewed and approved by the Ethics Committee of Yangpu Hospital. Informed consent was provided by all subjects enrolled in this study. The ethical approval number is LL-2018-WSJ-010. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

YC and CY designed the work and prepared the manuscript. YJ collected the samples. YJ, XZ, and XYZ performed the experiments, conducted data analyses, and wrote the manuscript. KZ and JZ analyzed and interpreted the data. All authors discussed the results, and read and approved the final version of the manuscript for publication.

Funding

This study was supported by the Scientific Research Foundation of Shanghai Municipal Commission of Health and Family Planning (Grant numbers 201840336 and 20194Y0238) and The Project of Key Disciplines of Medicine in Yangpu District of Shanghai (Grant number YP19ZA01).

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.2021.670528/full#supplementary-material

Supplementary Figure 1 | Serval enrichment analyses of pterygium. (A) Top five m6A motifs from the altered m6A peaks in conjunctiva. (B) Top 15 BP, MF, and CC GO terms of genes with (up panel) upregulated m6A peaks and (down panel) downregulated peaks. (C,D) Top 20 (C) GO enrichment and (D) KEGG enrichment terms of genes with significant changes in both m6A modification and mRNA levels. GO, Gene Ontology; BP, biological processes; CC, cellular components; MF, molecular functions.

Supplementary Figure 2 | Results of visualization analyses of five hub genes. (A–E) The m6A level and expression of five hub genes in each tissue type were visualized using Integrative Genomics Viewer.

Supplementary Table 1 | The detail of clinical information about 24 pterygium patients.

Supplementary Table 2 | List of primers used for qRT-PCR.

Supplementary Table 3 | Most relevant networks of lncRNA and hyper- or hypo-methylated peaks.

Supplementary Table 4 | Most relevant networks of lncRNA from RNA-seq.

Supplementary Table 5 | Most relevant networks of m6A-modified genes.

Footnotes

1. ^ https://portal.genego.com
2. ^ https://www.ncbi.nlm.nih.gov/gds/
3. ^ www.r-project.org

References

Ang, M., Li, X., Wong, W., Zheng, Y., Chua, D., Rahman, A., et al. (2012). Prevalence of and racial differences in pterygium: a multiethnic population study in Asians. Ophthalmology 119, 1509–1515. doi: 10.1016/j.ophtha.2012.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Arnemann, J., Spurr, N. K., Wheeler, G. N., Parker, A. E., and Buxton, R. S. (1991). Chromosomal assignment of the human genes coding for the major proteins of the desmosome junction, desmoglein DGI (DSG), desmocollins DGII/III (DSC), desmoplakins DPI/II (DSP), and plakoglobin DPIII (JUP). Genomics 10, 640–645. doi: 10.1016/0888-7543(91)90446-l

CrossRef Full Text | Google Scholar

Aster, J. C., Pear, W. S., and Blacklow, S. C. (2017). The Varied Roles of Notch in Cancer. Annu. Rev. Pathol. 12, 245–275. doi: 10.1146/annurev-pathol-052016-100127

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, K. H., He, S. Y., Shu, L. L., Wang, W. D., Lin, S. Y., Zhang, Q. Y., et al. (2020). Identification of cancer stem cell characteristics in liver hepatocellular carcinoma by WGCNA analysis of transcriptome stemness index. Cancer Med. 9, 4290–4298. doi: 10.1002/cam4.3047

PubMed Abstract | CrossRef Full Text | Google Scholar

Baquero-Perez, B., Antanaviciute, A., Yonchev, I. D., Carr, I. M., Wilson, S. A., and Whitehouse, A. (2019). The Tudor SND1 protein is an m(6)A RNA reader essential for replication of Kaposi’s sarcoma-associated herpesvirus. Elife 8:261. doi: 10.7554/eLife.47261

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbieri, I., Tzelepis, K., Pandolfini, L., Shi, J., Millan-Zambrano, G., Robson, S. C., et al. (2017). Promoter-bound METTL3 maintains myeloid leukaemia by m(6)A-dependent translation control. Nature 552, 126–131. doi: 10.1038/nature24678

PubMed Abstract | CrossRef Full Text | Google Scholar

Baskin, B., Skinner, J. R., Sanatani, S., Terespolsky, D., Krahn, A. D., Ray, P. N., et al. (2013). TMEM43 mutations associated with arrhythmogenic right ventricular cardiomyopathy in non-Newfoundland populations. Hum. Genet. 132, 1245–1252. doi: 10.1007/s00439-013-1323-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertoni, N., Pereira, L. M., Severino, F. E., Moura, R., Yoshida, W. B., and Reis, P. P. (2016). Integrative meta-analysis identifies microRNA-regulated networks in infantile hemangioma. BMC Med. Genet. 17:4. doi: 10.1186/s12881-015-0262-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Boccaletto, P., Machnicka, M., Purta, E., Piatkowski, P., Baginski, B., Wirecki, T., et al. (2018). MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic Acids Res. 46, D303–D307. doi: 10.1093/nar/gkx1030

PubMed Abstract | CrossRef Full Text | Google Scholar

Bornslaeger, E. A., Corcoran, C. M., Stappenbeck, T. S., and Green, K. J. (1996). Breaking the connection: displacement of the desmosomal plaque protein desmoplakin from cell-cell interfaces disrupts anchorage of intermediate filament bundles and alters intercellular junction assembly. J. Cell Biol. 134, 985–1001. doi: 10.1083/jcb.134.4.985

PubMed Abstract | CrossRef Full Text | Google Scholar

Brockdorff, N., and Turner, B. M. (2015). Dosage compensation in mammals. Cold Spring Harb. Perspect Biol. 7:a019406. doi: 10.1101/cshperspect.a019406

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Maqsood, S., Kaye, S., Tey, A., and Ahmad, S. (2014). Pterygium: are we any closer to the cause? Br. J. Ophthalmol. 98, 423–424. doi: 10.1136/bjophthalmol-2013-304232

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., Wei, L., Law, C. T., Tsang, F. H., Shen, J., Cheng, C. L., et al. (2018). RNA N6-methyladenosine methyltransferase-like 3 promotes liver cancer progression through YTHDF2-dependent posttranscriptional silencing of SOCS2. Hepatology 67, 2254–2270. doi: 10.1002/hep.29683

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Li, Y., Zhi, S., Ding, Z., Wang, W., Peng, Y., et al. (2020). WTAP promotes osteosarcoma tumorigenesis by repressing HMBOX1 expression in an m(6)A-dependent manner. Cell Death Dis. 11:659. doi: 10.1038/s41419-020-02847-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., Huang, B., Zhu, L., Chen, K., Liu, M., and Zhong, C. (2020). Structural and Functional Overview of TEAD4 in Cancer Biology. Oncol. Targets Ther. 13, 9865–9874. doi: 10.2147/ott.s266649

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, T., Ding, L., Shan, G., Ke, L., Ma, J., and Zhong, Y. (2015). Prevalence and racial differences in pterygium: a cross-sectional study in Han and Uygur adults in Xinjiang. China. Invest. Ophthalmol. Vis. Sci. 56, 1109–1117. doi: 10.1167/iovs.14-15994

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X. Y., Zhang, J., and Zhu, J. S. (2019). The role of m(6)A RNA methylation in human cancer. Mol. Cancer 18:103. doi: 10.1186/s12943-019-1033-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Hong, T., Wang, S., Mo, J., Tian, T., and Zhou, X. (2017). Epigenetic modification of nucleic acids: from basic studies to medical applications. Chem. Soc. Rev. 46, 2844–2872. doi: 10.1039/c6cs00599c

PubMed Abstract | CrossRef Full Text | Google Scholar

Chien, K. H., Chen, S. J., Liu, J. H., Woung, L. C., Chen, J. T., Liang, C. M., et al. (2013). Correlation of microRNA-145 levels and clinical severity of pterygia. Ocul. Surf. 11, 133–138. doi: 10.1016/j.jtos.2012.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Cimpean, A., Sava, M., and Raica, M. (2013). DNA damage in human pterygium: one-shot multiple targets. Mol. Vis. 19, 348–356.

Google Scholar

Cui, Q., Shi, H., Ye, P., Li, L., Qu, Q., Sun, G., et al. (2017). m(6)A RNA Methylation Regulates the Self-Renewal and Tumorigenesis of Glioblastoma Stem Cells. Cell Rep. 18, 2622–2634. doi: 10.1016/j.celrep.2017.02.059

PubMed Abstract | CrossRef Full Text | Google Scholar

De Francesco, E. M., Maggiolini, M., and Musti, A. M. (2018). Crosstalk between Notch, HIF-1α and GPER in Breast Cancer EMT. Int. J. Mol. Sci. 19:11. doi: 10.3390/ijms19072011

PubMed Abstract | CrossRef Full Text | Google Scholar

Dorn, L. E., Lasman, L., Chen, J., Xu, X., Hund, T. J., Medvedovic, M., et al. (2019). The N(6)-Methyladenosine mRNA Methylase METTL3 Controls Cardiac Homeostasis and Hypertrophy. Circulation 139, 533–545. doi: 10.1161/CIRCULATIONAHA.118.036146

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernandes, A., Salomão, S., Ferraz, N., Mitsuhiro, M., Furtado, J., Muñoz, S., et al. (2020). Pterygium in adults from the Brazilian Amazon Region: prevalence, visual status and refractive errors. Br. J. Ophthalmol. 104, 757–763. doi: 10.1136/bjophthalmol-2019-314131

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Y., Dominissini, D., Rechavi, G., and He, C. (2014). Gene expression regulation mediated through reversible m6A RNA methylation. Nat. Rev. Genet. 15, 293–306. doi: 10.1038/nrg3724

PubMed Abstract | CrossRef Full Text | Google Scholar

Gustavsson, J., Mehlig, K., Leander, K., Lissner, L., Bjorck, L., Rosengren, A., et al. (2014). FTO genotype, physical activity, and coronary heart disease risk in Swedish men and women. Circ. Cardiovasc. Genet 7, 171–177. doi: 10.1161/CIRCGENETICS.111.000007

PubMed Abstract | CrossRef Full Text | Google Scholar

He, L., Li, H., Wu, A., Peng, Y., Shu, G., and Yin, G. (2019). Functions of N6-methyladenosine and its role in cancer. Mol. Cancer 18:176. doi: 10.1186/s12943-019-1109-9

PubMed Abstract | CrossRef Full Text | Google Scholar

He, S., Li, X., Chan, N., and Hinton, D. R. (2013). Review: Epigenetic mechanisms in ocular disease. Mol. Vis. 19, 665–674.

Google Scholar

Héraud, C., Pinault, M., Lagrée, V., and Moreau, V. (2019). p190RhoGAPs, the ARHGAP35- and ARHGAP5-Encoded Proteins, in Health and Disease. Cells 8:351. doi: 10.3390/cells8040351

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, A., Lan, W., Law, K. P., Khoo, S. C., Tin, M. Q., Lim, Y. P., et al. (2014). Evaluation of global differential gene and protein expression in primary Pterygium: S100A8 and S100A9 as possible drivers of a signaling network. PLoS One 9:e97402. doi: 10.1371/journal.pone.0097402

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, Y., Wang, S., Liu, J., Huang, Y., Gong, C., Liu, J., et al. (2020). New sights in cancer: Component and function of N6-methyladenosine modification. Biomed. Pharmacother. 122:109694. doi: 10.1016/j.biopha.2019.109694

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, da, W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, C., Zhu, Y., Zhou, Z., Gumin, J., Bengtsson, L., Wu, W., et al. (2017). TMEM43/LUMA is a key signaling component mediating EGFR-induced NF-kappaB activation and tumor progression. Oncogene 36, 2813–2823. doi: 10.1038/onc.2016.430

PubMed Abstract | CrossRef Full Text | Google Scholar

Jordan, W. III, Rieder, L. E., and Larschan, E. (2019). Diverse Genome Topologies Characterize Dosage Compensation across Species. Trends Genet. 35, 308–315. doi: 10.1016/j.tig.2019.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawaguchi, Y., Cooper, B., Gannon, M., Ray, M., MacDonald, R. J., and Wright, C. V. (2002). The role of the transcriptional regulator Ptf1a in converting intestinal to pancreatic progenitors. Nat. Genet. 32, 128–134. doi: 10.1038/ng959

PubMed Abstract | CrossRef Full Text | Google Scholar

Khanna, R., Marmamula, S., Cicinelli, M., Mettla, A., Giridhar, P., Banerjee, S., et al. (2020). Fifteen-year incidence rate and risk factors of pterygium in the Southern Indian state of Andhra Pradesh. Br. J. Ophthalmol. 2020:359. doi: 10.1136/bjophthalmol-2020-316359

PubMed Abstract | CrossRef Full Text | Google Scholar

Koo, J. H., and Guan, K. L. (2018). Interplay between YAP/TAZ and Metabolism. Cell Metab. 28, 196–206. doi: 10.1016/j.cmet.2018.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y., Choe, J., Park, O. H., and Kim, Y. K. (2020). Molecular Mechanisms Driving mRNA Degradation by m(6)A Modification. Trends Genet. 36, 177–188. doi: 10.1016/j.tig.2019.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H. B., Tong, J., Zhu, S., Batista, P. J., Duffy, E. E., Zhao, J., et al. (2017). m(6)A mRNA methylation controls T cell homeostasis by targeting the IL-7/STAT5/SOCS pathways. Nature 548, 338–342. doi: 10.1038/nature23450

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, L., Sun, J., Wu, D., Lin, D., Sun, D., Li, Q., et al. (2019). MicroRNA-186 is associated with hypoxia-inducible factor-1α expression in chronic obstructive pulmonary disease. Mol. Genet. Genomic Med. 7:e531. doi: 10.1002/mgg3.531

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, S., Choe, J., Du, P., Triboulet, R., and Gregory, R. I. (2016). The m(6)A Methyltransferase METTL3 Promotes Translation in Human Cancer Cells. Mol. Cell. 62, 335–345. doi: 10.1016/j.molcel.2016.03.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Dou, X., Chen, C., Chen, C., Liu, C., Xu, M. M., et al. (2020). N (6)-methyladenosine of chromosome-associated regulatory RNA regulates chromatin state and transcription. Science 367, 580–586. doi: 10.1126/science.aay6018

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Yue, Y., Han, D., Wang, X., Fu, Y., Zhang, L., et al. (2014). A METTL3-METTL14 complex mediates mammalian nuclear RNA N6-adenosine methylation. Nat. Chem. Biol. 10, 93–95. doi: 10.1038/nchembio.1432

PubMed Abstract | CrossRef Full Text | Google Scholar

Meshkani, S., Kooshan, N., Moghadam, A., Falanji, F., Adli, A., Baghbani-Arani, F., et al. (2019). Signaling roadmap to epithelial-mesenchymal transition in pterygium, TWIST1 centralized. J. Cell. Physiol. 234, 18146–18155. doi: 10.1002/jcp.28447

PubMed Abstract | CrossRef Full Text | Google Scholar

Meurette, O., and Mehlen, P. (2018). Notch Signaling in the Tumor Microenvironment. Cancer Cell 34, 536–548. doi: 10.1016/j.ccell.2018.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Mo, J. S., Park, H. W., and Guan, K. L. (2014). The Hippo signaling pathway in stem cell biology and cancer. EMBO Rep. 15, 642–656. doi: 10.15252/embr.201438638

PubMed Abstract | CrossRef Full Text | Google Scholar

Nemesure, B., Wu, S., Hennis, A., and Leske, M. (2008). Nine-year incidence and risk factors for pterygium in the barbados eye studies. Ophthalmology 115, 2153–2158. doi: 10.1016/j.ophtha.2008.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, D. (2010). The hippo signaling pathway in development and cancer. Dev. Cell 19, 491–505. doi: 10.1016/j.devcel.2010.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, X., Hong, X., Li, S., Meng, P., and Xiao, F. (2021). METTL3 promotes adriamycin resistance in MCF-7 breast cancer cells by accelerating pri-microRNA-221-3p maturation in a m6A-dependent manner. Exp. Mol. Med. 53, 91–102. doi: 10.1038/s12276-020-00510-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Papagerakis, S., Shabana, A. H., Pollock, B. H., Papagerakis, P., Depondt, J., and Berdal, A. (2009). Altered desmoplakin expression at transcriptional and protein levels provides prognostic information in human oropharyngeal cancer. Hum. Pathol. 40, 1320–1329. doi: 10.1016/j.humpath.2009.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Piccolo, S., Dupont, S., and Cordenonsi, M. (2014). The biology of YAP/TAZ: hippo signaling and beyond. Physiol. Rev. 94, 1287–1312. doi: 10.1152/physrev.00005.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Rezasoltani, S., Asadzadeh-Aghdaei, H., Nazemalhosseini-Mojarad, E., Dabiri, H., Ghanbari, R., and Zali, M. R. (2017). Gut microbiota, epigenetic modification and colorectal cancer. Iran. J. Microbiol. 9, 55–63.

Google Scholar

Robertson, I. B., Horiguchi, M., Zilberberg, L., Dabovic, B., Hadjiolova, K., and Rifkin, D. B. (2015). Latent TGF-beta-binding proteins. Matrix Biol. 47, 44–53. doi: 10.1016/j.matbio.2015.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Roundtree, I. A., Luo, G. Z., Zhang, Z., Wang, X., Zhou, T., Cui, Y., et al. (2017). YTHDC1 mediates nuclear export of N(6)-methyladenosine methylated mRNAs. Elife 6:311. doi: 10.7554/eLife.31311

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahakyan, A., Yang, Y., and Plath, K. (2018). The Role of Xist in X-Chromosome Dosage Compensation. Trends Cell. Biol. 28, 999–1013. doi: 10.1016/j.tcb.2018.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Schöller, E., Weichmann, F., Treiber, T., Ringle, S., Treiber, N., Flatley, A., et al. (2018). Interactions, localization, and phosphorylation of the m(6)A generating METTL3-METTL14-WTAP complex. RNA 24, 499–512. doi: 10.1261/rna.064063.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Schuierer, S., Tranchevent, L. C., Dengler, U., and Moreau, Y. (2010). Large-scale benchmark of Endeavour using MetaCore maps. Bioinformatics 26, 1922–1923. doi: 10.1093/bioinformatics/btq307

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, N., Guo, X., and Chen, S. Y. (2014). Olfactomedin 2, a novel regulator for transforming growth factor-beta-induced smooth muscle differentiation of human embryonic stem cell-derived mesenchymal cells. Mol. Biol. Cell. 25, 4106–4114. doi: 10.1091/mbc.E14-08-1255

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Q., Cibas, E. S., Huang, H., Hodgson, L., and Overholtzer, M. (2014). Induction of entosis by epithelial cadherin expression. Cell Res. 24, 1288–1298. doi: 10.1038/cr.2014.137

PubMed Abstract | CrossRef Full Text | Google Scholar

Tiwari, A., Swamynathan, S., Campbell, G., Jhanji, V., and Swamynathan, S. K. (2020). BMP6 Regulates Corneal Epithelial Cell Stratification by Coordinating Their Proliferation and Differentiation and Is Upregulated in Pterygium. Invest. Ophthalmol. Vis. Sci. 61:46. doi: 10.1167/iovs.61.10.46

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, G. H., Yao, L., Xu, H. W., Tang, W. T., Fu, J. H., Hu, X. F., et al. (2013). Identification of MXRA5 as a novel biomarker in colorectal cancer. Oncol. Lett. 5, 544–548. doi: 10.3892/ol.2012.1038

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Zhao, B. S., Roundtree, I. A., Lu, Z., Han, D., Ma, H., et al. (2015). N(6)-methyladenosine Modulates Messenger RNA Translation Efficiency. Cell 161, 1388–1399. doi: 10.1016/j.cell.2015.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiemann, S., Weil, B., Wellenreuther, R., Gassenhuber, J., Glassl, S., Ansorge, W., et al. (2001). Toward a catalog of human genes and proteins: sequencing and analysis of 500 novel complete protein coding human cDNAs. Genome Res. 11, 422–435. doi: 10.1101/gr.gr1547r

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, Y. W., Chew, J., Yang, H., Tan, D. T., and Beuerman, R. (2006). Expression of insulin-like growth factor binding protein-3 in pterygium tissue. Br. J. Ophthalmol. 90, 769–772. doi: 10.1136/bjo.2005.087486

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, D., Li, G., Li, K., Xu, Q., Pan, Z., Ding, F., et al. (2012). Exome sequencing identifies MXRA5 as a novel cancer gene frequently mutated in non-small cell lung carcinoma from Chinese patients. Carcinogenesis 33, 1797–1805. doi: 10.1093/carcin/bgs210

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Chen, Y., Cui, T., Knosel, T., Zhang, Q., Albring, K. F., et al. (2012). Desmoplakin acts as a tumor suppressor by inhibition of the Wnt/beta-catenin signaling pathway in human lung cancer. Carcinogenesis 33, 1863–1870. doi: 10.1093/carcin/bgs226

PubMed Abstract | CrossRef Full Text | Google Scholar

Yi, D., Wang, R., Shi, X., Xu, L., Yilihamu, Y., and Sang, J. (2020). METTL14 promotes the migration and invasion of breast cancer cells by modulating N6-methyladenosine and hsa-miR-146a-5p expression. Oncol. Rep. 43, 1375–1386. doi: 10.3892/or.2020.7515

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, F. X., Zhao, B., and Guan, K. L. (2015). Hippo Pathway in Organ Size Control, Tissue Homeostasis, and Cancer. Cell 163, 811–828. doi: 10.1016/j.cell.2015.10.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Zanconato, F., Forcato, M., Battilana, G., Azzolin, L., Quaranta, E., Bodega, B., et al. (2015). Genome-wide association between YAP/TAZ/TEAD and AP-1 at enhancers drives oncogenic growth. Nat. Cell. Biol. 17, 1218–1227. doi: 10.1038/ncb3216

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Shi, X., Huang, T., Zhao, X., Chen, W., Gu, N., et al. (2020a). Dynamic landscape and evolution of m6A methylation in human. Nucleic Acids Res. 48, 6251–6264. doi: 10.1093/nar/gkaa347

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Zhang, P., Long, C., Ma, X., Huang, H., Kuang, X., et al. (2020b). m(6)A methyltransferase METTL3 promotes retinoblastoma progression via PI3K/AKT/mTOR pathway. J. Cell. Mol. Med. 24, 12368–12378. doi: 10.1111/jcmm.15736

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Qiao, Y., Huang, J., Wan, D., Zhou, L., Lin, S., et al. (2020). Expression Pattern and Prognostic Value of Key Regulators for m6A RNA Modification in Hepatocellular Carcinoma. Front. Med. 7:556. doi: 10.3389/fmed.2020.00556

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Kuang, Y., Wang, Y., Xu, Q., and Ren, Q. (2017). Notch-4 silencing inhibits prostate cancer growth and EMT via the NF-κB pathway. Apoptosis 22, 877–884. doi: 10.1007/s10495-017-1368-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Kong, Y., Fan, W., Tao, T., Xiao, Q., Li, N., et al. (2020). Principles of RNA methylation and their implications for biology and medicine. Biomed. Pharmacother. 131:110731. doi: 10.1016/j.biopha.2020.110731

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, B., Gong, Y., Shen, L., Li, J., Han, J., Song, B., et al. (2020). Total Panax notoginseng saponin inhibits vascular smooth muscle cell proliferation and migration and intimal hyperplasia by regulating WTAP/p16 signals via m(6)A modulation. Biomed. Pharmacother. 124:109935. doi: 10.1016/j.biopha.2020.109935

PubMed Abstract | CrossRef Full Text | Google Scholar