Comprehensive analysis of transcriptome-wide m⁶A methylome in the anterior capsule of the lens of high myopia patients

ABSTRACT

To assess the m⁶A methylome in the anterior capsule of the lens of high myopia patients. MeRIP-seq and RNA-seq were performed to identify differences in the m⁶A methylomes and gene expression between anterior capsule of the lens of simple nuclear cataract patients (N) and nuclear cataract patients with high myopia (G). Expression of m⁶A-related enzymes was confirmed by quantitative real-time-PCR. ALKBH5 was downregulated in G. The observed m⁶A peak was identical to the conserved RRACH gmotif and was markedly correlated with two distinct coordinates. Differentially methylated genes were enriched in some pathways regulating the formation of extracellular matrix. These findings suggest that upregulation of m⁶A methylation may change fundus anatomy by regulating the composition of the extracellular matrix through encoding protein.

KEYWORDS:

• High myopia
• m⁶A methylome
• extracellular matrix
• MeRIP-seq

Introduction

With the change in lifestyle, great modifications have taken place in human eye-use habits, and the incidence of high myopia is increasing. At present, 163 million of people, who account for 2.7% of the world total population, are suffering from high myopia [1–3]. Worryingly, this number is rising sharply [2–4]. High myopia, defined as myopia exceeding 6.00 dioptres or axis length ≥26 mm, is a disorder that affects almost the entire human eye, from the anterior pole to the posterior pole, including high ametropia, cataract, open-angle glaucoma, and retinopathy [5]. Excessive axial elongation of the eye in high myopia can cause mechanical stretching of the outer coats of the eyeball, resulting in various pathologic changes such as staphyloma, chorioretinal atrophic lesions, lacquer cracks, and choroidal neovascularization.

N⁶-methyladenosine (m⁶A), as the most prevalent internal form of modification in polyadenylated mRNAs and long noncoding RNAs in higher eukaryotes, was first identified in the 1970s [6]. Studies in bioscience and medicine have revealed that RNA m⁶A plays important biological roles in the regulation of cellular metabolic processes. In 2014, Batista et al. have found that RNA m⁶A controls cell transition fate in mammalian embryonic stem cells [7]. The next year, Zhao et al. have discovered that RNA m⁶A regulates pluripotency in murine stem cells [8]. In addition, Shen et al. have studied the role of m⁶A in shoot stem cell fate in Arabidopsis in 2016 [9]. Furthermore, m⁶A methylation plays important roles in human disease such as control of HIV-1 replication and interaction with the host immune system during HIV-1 infection of T cells [10], promoting translation of oncogenes in human lung cancer [11] and induction of breast cancer stem cell phenotype [12]. Thus, more medical fields will reveal the role of m⁶A methylation. In mammals, m⁶A methylation modification involves three enzymes: methylase (METTL3, METTL14, and WTAP), demethylase (FTO and ALKBH5), and methylation recognition enzyme (YTHDF1, YTHDF2, and YTHDF3). Enzyme abnormality causes a series of diseases, including tumour, neurological diseases, and circulatory disorders [13–15]. However, the role of m⁶A methylation in high myopia is still unclear.

To investigate the functions of m⁶A, facilitate future studies of mammalian m⁶A, and explore the pathogenic mechanism of the eye caused by high myopia, we detected the m⁶A methylomes of anterior lens capsule in high myopia patients. This enabled us to acquire a set of transcriptome-wide m⁶A profiles in patients with high myopia and investigate tissue and breed generality and selectivity of methylated genes and their functional implications.

Materials and methods

Acquisition of biological samples

Twelve anterior lens capsule samples were enrolled and divided into two groups. One group was characterized by simple nuclear cataract patients (N, 6 samples), while the other group was characterized by nuclear cataract patients with high myopia (G, 6 samples). Each patient underwent conventional cataract phacoemulsification combined with intraocular lens implantation. During continuous circular capsulorhexis, the anterior capsule was collected into a cryopreservation tube and stored in −80ºC. Patients with nuclear cataract and an axis length ≥29 mm were included in this study. Whereas patients with other eye diseases and other systemic diseases such as hypertension, diabetes, coronary heart disease, and autoimmune disease were excluded from the investigation. Patient information is shown in Table 1. Informed consent was obtained from all individual participants included in the study. All procedures performed in studies involving human participants were in accordance with the ethical standards of the Institutional and National Research Committee and with the 1964 Helsinki Declaration and its later amendments or comparable ethical standards.

Table 1. Patient characteristics of G and N

	Gender		Ages		Axial		Eye
Samples	N	G	N	G	N	G	N	G
1	F	M	73	60	23.69	29.20	OS	OS
2	F	M	59	78	23.67	29.69	OD	OS
3	F	M	73	59	23.17	29.95	OD	OD
4	F	M	58	59	25.42	29.71	OS	OS
5	F	M	68	78	22.92	29.56	OS	OS
6	M	F	61	70	24.68	29.32	OS	OD
P			0.680		0.000
Mean			65.33 ± 6.99	67.33 ± 9.24	23.93 ± 0.95	29.57 ± 0.27

N: simple nuclear cataract; G: nuclear cataract complicated with high myopia;

F: female; M: male; OD: right eye; OS: left eye.

Preparation and sequencing of RNA library

RNA high-throughput sequencing was performed by Cloud-Seq Biotech (Shanghai, China). Briefly, total RNA was used for removing the rRNAs using the NEBNext® rRNA Depletion kit (New England Biolabs, Inc., MA, USA) following the manufacturer’s instructions. RNA libraries were constructed by using NEBNext® Ultra™ II Directional RNA Library Prep kit (New England Biolabs) according to the manufacturer’s instructions. Libraries were controlled for quality and quantified using the Bioanalyzer 2100 system (Agilent Technologies, Inc., USA). Library sequencing was performed on an illumina Hiseq instrument with 150 bp paired-end reads.

Preparation and sequencing of MeRIP library

m⁶A RNA-seq service was provided by Cloud-seq Biotech Inc. Briefly, m⁶A RNA immunoprecipitation was performed using the GenSeq™ m⁶A RNA IP kit (GenSeq Inc., China) by following the manufacturer’s instructions. The quality and quantity of total RNA were assessed by using NanoDrop ND-2000 (Thermo Fisher Scientific, USA). RNA purity was qualified when OD₂₆₀/OD₂₈₀ values ranged from 1.8 to 2.1. RNA integrity and gDNA contamination were measured using denatured agarose gel electrophoresis. Agilent 2100 Bioanalyzer was used to detect library quality. RNA quantification and quality assurance were assessed using Qubit3. Both the input sample without immunoprecipitation and the m⁶A IP samples were used for RNA-seq library generation using NEBNext® Ultra II Directional RNA Library Prep kit (New England Biolabs). The library quality was evaluated using the BioAnalyzer 2100 system (Agilent Technologies, Inc.). Library sequencing was performed on an illumina Hiseq instrument with 150 bp paired-end reads.

Data analysis

Data analysis of mRNA sequencing results

Paired-end reads were obtained from the Illumina HiSeq 4000 sequencer and were quality-controlled using Q30. After 3ʹ adaptor-trimming and low-quality reads removing using the cutadapt software (v1.9.3) [16], the high-quality clean reads were aligned to the reference genome (UCSC hg19) using the Hisat2 software (v2.0.4) [17]. Guided by the Ensembl gtf gene annotation file, the cuffdiff software [18] (part of cufflinks) was used to get the gene-level FPKM as the expression profiles of mRNA, fold change and p-value were calculated based on FPKM, and differentially expressed mRNA were identified. Gene ontology (GO) and pathway enrichment analysis were performed based on the differentially expressed mRNAs. GO and pathway enrichment analysis were performed using the Database for Annotation, Visualization, and Integrated Discovery [19,20]. The ontology covers three parts: cellular component (CC), molecular function (MF), and biological process (BP). The p-value denotes the significance of GO term enrichment of the genes. Pathway enrichment analysis is a functional analysis that maps genes to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The Fisher p-value denotes the significance of the pathway correlated to the conditions.

Data analysis of MeRIP sequencing results

Briefly, paired-end reads were obtained from the Illumina HiSeq 4000 sequencer and were quality-controlled using Q30. Then, 3ʹ adaptor-trimming and low quality reads removing using the cutadapt software (v1.9.3) [16] were performed. Clean reads of all libraries were aligned to the reference genome (HG19) using the Hisat2 software (v2.0.4) [17]. Reads were matched to the results on the genome (bam file), and the IGV software [21] was used for visualization to observe the abundance of reads at specific locations in the genome of each sample. Methylated sites on RNAs (peaks) were identified using the MACS software [22]. Differentially methylated sites were identified using diffReps [23]. These peaks identified by both software overlapping with exons of mRNA were identified and selected using home-made scripts. GO and pathway enrichment analysis were performed for the differentially methylated protein-coding genes. MEME software was used for motif analysis [24].

RNA extraction and quantitative real-time PCR

Levels of the mRNA m⁶A-related genes METTL3, METTL14, FTO, ALKBH5, YTHDF1, and YTHDF2 were analysed in G and N groups. TRIzol reagent (Invitrogen, USA) was used to isolate total RNA, which was then used to synthesize complementary DNA by using the SuperScript™ III Reverse Transcriptase (Invitrogen). Real-time PCR was performed by using qPCR SYBR Green Master Mix (CloudSeq, Shanghai) and QuantStudio 5 Real-Time PCR System (Thermo Fisher Scientific, USA). All procedures were performed according to the manufacturer’s protocols. The sequences of primers used are presented in Table 2.

Table 2. Primers used for quantitative real-time-PCR

Gene	Primer	Gene sequence
ALKBH5	Forward	CACGGATCCTGGAGATGG
	Reverse	TCTTCACCTTTCGGGCAG
METTL3	Forward	CAAAGGAGCCAGCCAAGA
	Reverse	CTGAAGTGCAGCTTGCGA
METTL14	Forward	CGAGCCTTGCTTGCAGTT
	Reverse	CAAGGCCCTGTCAGCACT
FTO	Forward	GAGTGGACCCAGCCCTCT
	Reverse	CTTGAAATGTGGCCAGGG
YTHDF1	Forward	TCAATGGGAGTGGGCATT
	Reverse	GTGTCCCGGGAGTTTGTG
YTHDF2	Forward	ACCCAGGGGTCTCCTCAG
	Reverse	CCATTACCATCCACCCCA
GAPDH	Forward	GGCCTCCAAGGAGTAAGACC
	Reverse	AGGGGAGATTCAGTGTGGTG

Statistical analyses

The significance of differences between G and N groups was tested by using unpaired two-tailed Student’s t-test. All statistical analyses were conducted using GraphPad Prism v7.00 software. The differences were considered significant if p-value was <0.05.

Results

Overview of m⁶A methylation map in anterior capsular of G lens

We calculated the methylation sites in the two groups and found that the number of methylation sites in the G group was higher than that in the N group (Figure 1a).

Figure 1. (a) The number of methylation sites and genes of G and N. (b)The methylation difference between N and G. (c) Distribution of differentially methylated m6A sites with significance in chromosomes of high myopia

More than 80% and 70% of the m⁶A peaks were consistently detected in three biological replicates of N and G groups, respectively. These recurrent peaks should be regarded as highly enriched m⁶A peaks for further analysis. The results showed that 13,617 m⁶A recurrent peaks among 8,196 expressed genes in N and 8,569 m⁶A recurrent peaks among 5896 expressed genes in G were detected (Figure 1a, b). We performed genome-wide profiling of m⁶A-modified mRNA in G and N. Compared to N, 2,083 significantly hypermethylated m⁶A peaks and 907 significantly hypomethylated m⁶A peaks were identified in G (fold change ≥2 and p ≤ 0.0001). The top 20 altered m⁶A peaks are listed in Table 3. The hypermethylated m⁶A peaks were identified in all chromosomes (especially chr1), while no hypomethylated m⁶A peak was found in chrY (Figure 1c).

Table 3. The top 20 differently methylated m⁶A peaks

Gene	Gene ID	Fold change	Regulation	Chromsome	Start	End	Peak_length	P-value
C11orf96	4279	44.76811594	Up	chr11	43,964,201	43,965,420	1219	7.14069E-14
CSF1	1365	11.87987988	Up	chr1	110,464,468	110,464,616	148	1.68025E-13
TMEM176B	24,762	35.12587413	Up	chr7	150,493,453	150,493,780	327	1.84248E-13
PTP4A3	26,335	9.891826923	Up	chr8	142,441,561	142,441,620	59	2.11684E-13
CHL1	16,617	36.8	Up	chr3	447,177	447,580	403	2.24665E-13
COL6A3	14,856	9.425531915	Up	chr2	238,258,789	238,258,852	63	2.57265E-13
COL6A3	14,857	10.49201278	Up	chr2	238,259,772	238,259,835	63	3.0255E-13
CHI3L1	1997	6.584013051	Up	chr1	203,152,768	203,152,919	151	3.40939E-13
PXDN	12,802	16.53529412	Up	chr2	1,635,658	1,636,780	1122	3.90611E-13
IGFBP4	10,837	11.05427975	Up	chr17	38,599,741	38,600,336	595	4.01326E-13
RIMS1	22,499	306.1	Down	chr6	72,892,441	72,892,852	411	1.40E-13
PTCHD4	22,367	22.209302	Down	chr6	48,078,784	48,078,800	16	1.757E-13
CABP7	16,282	9.1723077	Down	chr22	30,116,343	30,116,522	179	2.625E-13
OCLN	20,678	334.8	Down	chr5	68,849,396	68,849,820	424	3.016E-13
PTCHD4	22,355	13.119048	Down	chr6	47,847,061	47,847,420	359	3.384E-13
ANKRD24	12,139	6.303719	Down	chr19	4,224,424	4,224,811	387	4.95E-13
ALS2CL	17,103	5.9847036	Down	chr3	46,711,821	46,712,140	319	6.597E-13
HEY2	22,927	4.247495	Down	chr6	126,073,152	126,073,236	84	7.109E-13
FGF10	20,490	10.644231	Down	chr5	44,389,281	44,389,740	459	7.967E-13
HEY2	22,926	7.2871901	Down	chr6	126,070,731	126,071,005	274	8.221E-13

To understand the preferential location of m⁶A in mRNA, we investigated the metagene profiles of m⁶A peaks in the mRNA transcriptome. m⁶A peaks were markedly correlated with two distinct coordinates: immediately following near the end of the 5ʹ untranslated regions (5ʹUTRs) and start of the coding sequence (CDS) and near the end of the CDS and the beginning of the 3ʹ untranslated region (3ʹUTRs) (Figure 2a).

Figure 2. (a) Preferential location of m6A in mRNA. Each transcript is divided into three parts including 5ʹuntranslated region, coding DNA sequence and 3ʹuntranslated region. (b) Pie charts showing m⁶A peaks distribution in different gene context. (c) Data visualization of MYL12A mRNA m⁶A modification in G

The stop codon of peaks was more pronounced than the start codon of peaks (31.63% vs 27.10%). To assess the enrichment methodically, we assigned each m⁶A peak to one of the five non-overlapping transcript segments: 5ʹUTRs, start codon, CDS, stop codon, and 3ʹUTR. As a result, 87% of m⁶A peaks were within genic regions, and more than 60% of genic peaks were localized near the stop codon and CDS, while 12% were found in the 5ʹUTRs and 3ʹUTRs (Figure 2b). The topological patterns distributing within genes were highly similar in both samples, suggesting that recognition of motif for m⁶A methylation was conserved among anterior capsule of human lens. MYL12A, a considerably hypermethylated peak, is shown in Figure 2c.

Notably, more than 3,500 of m⁶A-methylated coding genes contained only one m⁶A peak, while a relatively small number of genes contain two or more peaks (Figure 3a), which was consistent with the trend of the proportions previously reported in the mouse brain and pig liver [25].

Figure 3. (a) Number of m⁶A peaks per gene. more than 3500 m⁶A-methylated coding genes contained only one m⁶A peak. (b) The top five motifs enriched across m⁶A peaks identified from G

An unbiased search for motifs enriched in regions surrounding m⁶A peaks was performed to determine whether the identified m⁶A peak was identical to the conserved RRACH motif (where R represented purine, A was m⁶A, and H was a non-guanine base) [26,27]. Clustering of significantly enriched sequences perfectly confirmed the previously established m⁶A RRACH consensus sequence in both groups (Figure 3b). The identification of a strong consensus reinforces the authenticity of the discovered m⁶A peaks and supports the existence of a predominant methylation machinery.

m⁶A-containing genes are involved in important biological pathways

To further determine general biological functional pathways that the significance of m⁶A methylation in high myopia development, the genes containing significantly altered m⁶A peaks (differentially methylated genes, G) were analysed by performing GO and KEGG pathway analysis. We systematically screened these different and common peaks and the related genes and identified the GO terms with the help of the GO consortium database (Figure 4).

Figure 4. (a–c) Major gene ontology terms were significantly enriched for the hypermethylated genes. (d–f) Major gene ontology terms were significantly enriched for the hypomethylated genes

GO analysis showed that the hypermethylated m⁶A peak-related genes encoding m⁶A-containing mRNAs were mainly involved in a variety of biological processes, including anatomical structure morphogenesis (ontology: biological process; Figure 4a), proteinaceous extracellular matrix (ontology: cellular component; Figure 4b), extracellular matrix structural constituent (ontology: molecular function; Figure 4c). The hypomethylated genes were significantly associated with regulation of ion transport (ontology: biological process; Figure 4d), plasma membrane part (ontology: cellular component; Figure 4e), and ion channel activity (ontology: molecular function; figure 4f). GO biological process classification indicated that most of the m⁶A methylation was enriched in the cellular process (Figure 4a). These genes were involved in the formation of extracellular matrix (Figure 4c), suggesting that extracellular matrix was associated with the formation of pathological damage to high myopia.

Outstandingly, KEGG pathway analysis showed the result of unique G and unique N (Figure 5c), in addition to information on hypermethylated genes (Figure 5a) and hypomethylated genes (Figure 5b).

Figure 5. (a) The top ten significantly enriched pathways for the hypermethylated genes. (b) The top ten significantly enriched pathways for the hypomethylated genes. (c) Unique pathway of G and N

Conjoint analysis of RNA-seq and MeRIP-seq

Transcriptome profiles of altered genes in the anterior capsule of patients with high myopia were determined using RNA-seq and shown in Figure 6a. Compared to N, 20,385 genes were differentially expressed in G (fold change ≥1.5 and p < 0.05), including 7,834 upregulated genes and 1,2524 downregulated genes. The top 20 altered genes are listed in Table 3. The top 10 GO and KEGG pathways are shown in Supplementary Figures 1 and 2. Based on conjoint analysis of differentially expressed genes and differentially methylated genes, the higher the transcription expression, the higher the methylation ratio (Figure 6b, c).

Figure 6. (a) Scatter plots showing the differentially expressed genes (fold changes ≥2 and p < 0.05). (b) m⁶A and gene expression of transcript. (c) Data visualization of methylation and gene expression

ALKBH5 and FTO were downregulated in G

Based on quantitative real-time-PCR, the expressions of six major enzymes involved in m⁶A methylation, including METTL3, METTL14, FTO, ALKBH5, YTHDF1, and YTHDF2, were verified between G and N. As shown in Figure 7, in G, the levels of demethylases ALKBH5 and FTO, demethylases acting to remove m⁶A modification, significantly decreased compared to those in N (ALKBH5:G/N = 0.245, P < 0.05; FTO:G/N = 0.234, P < 0.05), whereas METTL3 and METTL14, the key methyltransferases responsible for m⁶A modifications, significantly decreased (G/N = 0.346, P < 0.05) and increased (G/N = 2.565, P < 0.05), respectively. The levels of two methylation recognition enzymes (YTHDF1 and YTHDF2) significantly decreased in G compared to N (YTHDF1: G/N = 0.284, P < 0.05; YTHDF2: G/N = 0.226, P < 0.05).

Figure 7. The expressions of six major enzymes involved with m6A methylation including METTL3, METTL14, FTO, ALKBH5, YTHDF1, and YTHDF2 between G and N

Expression of candidate genes correlated with pathological damage in high myopia

To further investigate the mechanisms of pathological damage in high myopia, we performed statistical analysis of differentially expressed m⁶A methylation-modified genes, including C11orf96 (uncharacterized protein), CSF1 (macrophage colony-stimulating factor 1), PTP4A3 (protein tyrosine phosphatase type IVA 3), COL6A3 (collagen alpha-3(VI) chain), CHI3L1 (chitinase-3-like protein 1), PXDN (peroxidasin homolog). However, COL6A3, CHI3L1, and PXDN were most likely to be our target genes.

Discussion

High myopia is a blinding disease with a high incidence, which seriously affects the quality of life of patients. In this study, we obtained the first high myopia transcriptome-wide m⁶A modification profile using MeRIP-Seq and discovered that mRNA m⁶A sites were mainly enriched around stop codons, CDS, and 3ʹUTRs, consistent with the distribution characteristics of mammal transcriptomes, although the modifications also occurred in 5ʹUTRs [7,25,28]. The m⁶A distribution in our study was the typical m⁶A topological pattern in most of the mature mRNA. The extensively higher m⁶A signals at the stop codon or 3ʹUTRs contributed to RNA stability, signalling for transport and translocation, or acted as regulatory elements for protein translation through the recruitment of specific factors onto the m⁶A sites for RNA transport or protein synthesis [29,30]. High-throughput m⁶A RNA sequencing databases showed that the distribution of m⁶A modifications on mRNA was sequence-specific and tended to occur in the conserved motif RRACH [25,28,31–35]. Accordingly, in our study, we successfully identified the consensus motif sequence in patients with high myopia transcriptome.

m⁶A mRNA modifications can be dynamically regulated by a multicomponent methyltransferase complex with opposite modifying activities, including m⁶A writers, m⁶A erasers, and m⁶A readers. Increasing evidence suggests that m⁶A modification plays a significant role in the proliferation, migration, and invasion of tumour cells. METTL3 was reported to be an oncogene in human lung cancer cells by promoting the translation of certain mRNAs such as epidermal growth factor receptor and Hippo pathway effector TAZ [11]. A study verified that METTL14 interacted with the microprocessor protein DGCR8 and positively modulated the primary microRNA 126 process in an m⁶A-dependent manner [36]. In ophthalmology, Jia et al. have found that YTHDF1 suppresses ocular melanoma by modulating HINT2 mRNA translation [37]. Luo et al. have confirmed that silencing METTL3 significantly suppressed uveal melanoma cell proliferation and colony formation through cell cycle G1 arrest, as well as migration and invasion [38]. Recently, it was reported that METTL3 silencing promoted the proliferation and alleviated the apoptosis of high glucose-induced human lens epithelial cells [39]. These studies have proved that RNA m⁶A methylation also plays an important role in the development of ophthalmic diseases. In our study, METTL14 was upregulated and METTL3, FTO, and ALKBH5 were downregulated in the anterior capsule of high myopia patient lens as determined using qPCR. This inconsistency of changes disrupts the dynamic equilibrium of writers and erasers of m⁶A modification. Our data suggest that m⁶A methylation is strongly associated with the pathogenic mechanism of high myopia.

RNA m⁶A methylation plays a key role in the regulation of post-transcriptional gene expression [25,31,40,41]. The importance of m⁶A in post-transcriptional regulation of gene expression is further reinforced by the discovery and characterization of mammalian reader proteins that recognize m⁶A modifications of mRNA and subsequently affect the stability of the target transcripts. YTHDF1 and YTHDF2, members of YTH domain-containing proteins, closely resembled each other and were predominantly cytoplasmic, but their mechanisms could be different. YTHDF2 mediated the decay of target mRNAs, whereas YTHDF1 interaction with eIF3 and other translation initiation factors suggested that it might affect translation rather than affect the half-lives of mRNAs [42,43]. In our study, compared to N, we confirmed that the YTHDF1 and YTHDF2 were downregulated in G. This result indicated that two readers might participate in the progress of high myopia through regulation of post-transcriptional gene expression. The gene FGF10 (Figure 6c), as a differentially methylated region, may have a certain impact on its downstream proteins, including nElavl and RBM10, which are RNA-binding proteins that directly interact with target RNAs and regulate several aspects of RNA metabolism. Studies have reported that nELAV proteins are unique to neurons and associated with Alzheimer’s disease [44]. Therefore, we speculate that hypomethylated FGF10 may enhance the ability to bind nELAV protein, playing an important role in the mechanism of high myopia brain injury. RBM10, another FGF10-binding protein, promotes many transformation and hypoxia-associated processes and events, including angiogenesis [45–47]. The occurrence of choroidal neovascularization in high myopia may be related to hypomethylation, weakening the ability of FGF10 to bind to RBM10.

One of the main functions of m⁶A is to mediate mRNA degradation in mammalian cells [29,48–50], suggesting a possible negative relationship between the m⁶A methylation extent and the transcript level. However, this observation somewhat differed from our present result in the anterior capsule of human lens, which showed that most of the highly expressed transcripts were relatively more modified by m⁶A. These differences may be due to different methodologies, different biological species, or different tissue samples. The results further indicated that different tissues might possess different characteristics in m⁶A methylation sites, suggesting a regulatory role for m⁶A in gene expression.

Differential m⁶A methylation has proved responsible for tissue or organ differentiation and development. Genes encoding m⁶A-containing RNAs in adult mouse brain tissue are linked to neurodevelopmental and neurological disorders [25,31]. Moreover, m⁶A RNA methylation in Drosophila and Zebrafish early embryogenesis shows a conserved mechanism of neuronal mRNA regulation contributing to brain function [43]. Our study uncovers regulation roles of m⁶A modification in the anterior lens capsule of high myopia patients. Differentially methylated genes mostly participate in the regulation of anatomical structure morphogenesis pathways. GO analysis showed that the upregulated m⁶A peak-related genes encoding m⁶A-containing mRNAs were mainly involved in the formation of extracellular matrix, consistent with the results shown by Wang [51]. In differentially expressed m⁶A methylation-modified genes, CHI3L1 (chitinase-3-like protein 1), encoding protein YKL-40, is thought to play a role in tissue remodelling and in the capacity of cells to respond to and cope with changes in their environment [52]. In addition, YKL-40 is an inflammatory marker that plays an important role in pathological processes such as cell proliferation, migration, differentiation, and tissue remodelling of inflammatory reactions [53–55]. Therefore, we suspect that the upregulation of methylase (METTL14) and the downregulation of demethylase (FTO and ALKBH5) catalyse the hypermethylation of gene CHI3L1, which may affect the expression level of its encoded protein YKL-40 and promotes the pathological state of high myopia by regulating the composition of the extracellular matrix. Posterior scleral staphyloma is one of the most basic changes in fundus anatomy in a series of interrelated degenerative changes in high myopia [56]. Thus, high expression of m⁶A methylation of mRNA may be involved in the formation of posterior scleral staphyloma by affecting the components of the extracellular matrix.

The choroid circulation receives approximately 95% of blood from the ophthalmic artery and provides oxygen and nutrition to the outer retinal layers. Therefore, the choroidal circulation may play an important role in the retinal dysfunction and vision loss. GO analysis showed that the upregulated m⁶A peak-related genes encoding m⁶A-containing mRNAs were also associated with circulatory system development, blood vessel development, and vasculature development (Figure 4). This suggests that the increase in methylation level of mRNA may cause damage to the fundus of patients with high myopia by affecting the choroidal circulation.

Conclusion

For the first time, we provide the integral human lens transcriptome m⁶A map in high myopia patients. Our map reveals features of m⁶A distribution in the high myopia of human transcriptome and identifies generality as well as selectivity of methylated genes and their functional implications. Additionally, the abnormal expression of methylation-related enzymes destroys the dynamic homoeostasis of methylation. This comprehensive methylome profiling provides a solid basis for the determination of potential functional roles for RNA m⁶A modification in pathological damage to the eye caused by high myopia.

Acknowledgments

The author thanks NewCore BiodataStudio in Shanghai for sequencing date analysis; Jing Sun, Yan Zhang, Yahong Li, and Qing Wang for important clinical and scientific suggestions. The study was supported by the talent project of Tianjin medical university eye hospital. No honorarium, grant, or other forms of payment were given to anyone to produce the manuscript. Each author listed on the manuscript has seen and approved the submission of this version of the manuscript and takes full responsibility for the manuscript.

Disclosure statement

The authors declare that they have no conflict of interest.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by the talent project of Tianjin medical university eye hospital under Grant YDYYRCXM-B2018-01LC; Tianjin clinical key discipline (specialist) construction project under Grant TJLCZDXKM004 and Tianjin clinical key discipline (specialist) construction project under Grant TJLCZDXKQ009.