Alteration of the m⁶A Methylation Landscape in a Mouse Model of Scleroderma

Abstract

Aim: To explore the N⁶-methyladenosine (m⁶A) methylation of mRNAs and its roles in a mouse model of scleroderma. Materials & methods: To evaluate whether the mouse model of scleroderma could meet the experimental requirements, we examined skin tissue specimens by pathological staining and identified the related indicators by quantitative PCR  (qPCR). m⁶A-tagged mRNAs were identified via m⁶A epitranscriptomic microarray, and m⁶A-RNA-immunoprecipitation qPCR and qPCR were performed to confirm microarray data. Results: There were differences in m⁶A methylation among 843 mRNAs. Further, there were significant differences among Hras, Saa1, Ccl3, Ccl9 and Il1b in terms of methylation and expression. Conclusion: The m⁶A methylation spectrum in a mouse model of scleroderma may explain the occurrence of scleroderma.

Lay abstract

Scleroderma is an autoimmune disease with an unknown etiology. At present there is no effective treatment, and the prognosis is poor. Although m⁶A methylation is the most common RNA modification and affects cell biological function by influencing mRNA expression levels, the role of m⁶A methylation in scleroderma is unknown. Based on a mouse model of scleroderma, we screened and validated mRNAs related to m⁶A regulation, including Hras, Lama3, Stat3, Tnc, Ccl3, Ccl9, Saa1 and Il1b. While the results are encouraging, they are only preliminary. We plan to explore the altered pattern of m⁶A methylation in skin specimens from patients with scleroderma. This study is important for developing biomarkers for diagnosis and treatment of scleroderma in the future.

Keywords::

• bleomycin
• methylation
• mRNA
• N⁶-methyladenosine
• scleroderma

Scleroderma is an autoimmune connective tissue disease characterized by organ fibrosis, immune abnormalities and vascular injury [1]. Depending on the extent of the lesion involved, the condition can present as localized scleroderma or systemic scleroderma [2]. Systemic scleroderma involves other tissues in addition to the skin; the lung, liver and kidney can also be involved, and destruction of the standard structure of the tissue may be observed. Fibrosis is the most typical characteristic of scleroderma; however, immune disorder and vascular injury precede and aggravate fibrosis [3]. Although scleroderma is a rare disease with low prevalence and morbidity, its mortality remains high [4]. Interstitial lung disease and pulmonary hypertension are the most common causes of death [5]. In addition, treatment for the condition has not changed significantly in the past 40 years and is still dominated by traditional symptomatic treatment [6]. As the pathogenesis of scleroderma remains unknown, it cannot be targeted to reverse the progression of fibrosis, highlighting unmet medical needs [7,8].

Epigenetics is a branch of genetics that deals with heritable changes in gene expression which do not involve any changes in the DNA sequence [9]. The genome contains two types of genetic information, namely DNA genetic information and epigenetic information, explaining why the limited genetic information has so many biological functions. For example, individuals can be genetically susceptible to scleroderma, but environmental factors such as viral infection and drug use are associated with the development of the condition. These factors may affect gene expression profiles through epigenetic mechanisms, thus altering various cells [10–12].

RNAs can be chemically modified to more than 100 distinct types that can enrich genetic diversity, such as N¹-methyladenosine (m¹A), N⁶-methyladenosine (m⁶A), N⁸-methyladenosine (m⁸A), N¹-methylguanosine (m¹G) and N²-methlyguanosine (m²G) [13–16]. RNA methylation affects every aspect of RNA regulation, including RNA stability, localization, transport, splicing and translation [17]. RNA methylation is a reversible process, and its dynamic equilibrium is dependent on the regulation of methyltransferases (METTL3, METTL14 and WTAP) and demethylases (FTO and ALKBH5), which are called ‘writers’ and ‘erasers’, respectively [18]. These enzymes regulate binding proteins (‘readers’), such as YTHDF1, YTHDF2, YTHDF3, YTHDC1 and YTHDC2, which can further identify RNA methylation modification information [19]. When the N⁶ position of adenine is methylated, it is referred to as m⁶A. m⁶A methylation is the most common RNA modification [20], accounting for about 50% of methylated ribonucleotides [21]. It is located near stop codons and in 3′ UTRs [22].

Although m⁶A methylation was discovered in the 1960s [23,24], research on m⁶A methylation has lagged behind compared with that of other modified bases, the main reason behind this being the lack of feasible methods to determine m⁶A sites [25]. Similarly, although high-throughput genomic technologies have revealed changes in noncoding RNAs (ncRNAs; miRNA, circRNA and long ncRNA [lncRNA]) and mRNA expression among scleroderma [26–29], few high-throughput epigenomic studies have reported on RNA m⁶A methylation in scleroderma [30]. Excitingly, scientists have found a way to combine m⁶A antibodies with m⁶A-containing RNA and perform high-throughput sequencing to identify the areas of RNA with high levels of m⁶A; this technique is called m⁶A-RNA-immunoprecipitation sequencing (MeRIP-seq) [20,25].

MeRIP-seq has helped scholars to reveal the role of abnormal m⁶A modification in diseases and particularly in oncology, including hepatocellular cancer [31], breast cancer [32] and pancreatic cancer [33]. In the same way, it is of great value to explore the role of m⁶A methylation in diseases outside of malignancy, especially fibrosis. Although the way in which m⁶A modification regulates fibrosis is not well understood, multiple studies have suggested that abnormal m⁶A modifications are involved in multiple organ fibrosis. In idiopathic pulmonary fibrosis, m⁶A methylation is upregulated in patients, a bleomycin (BLM)-induced pulmonary fibrosis mouse model and fibroblast-to-myofibroblast transition-derived myofibroblasts, and silencing METTL3 can suppress fibroblast-to-myofibroblast transition in vitro and in vivo [34]. In a mouse model of renal interstitial fibrosis, the m⁶A level is time-dependent and decreased within 1 week [18]. Similar alterations are observed in METTL3, METTL14 and WTAP [18]. In a mouse model of liver fibrosis, differential m⁶A-modified peaks are involved in oxidative stress and cytoplasmic metabolism during fibrogenesis, while differential m⁶A-modified peaks are enriched in immune response and apoptosis during fibrosis regression [9]. In a mouse model of myocardial infarction, Fto overexpression alleviates myocardial fibrosis [35]. However, its effects on skin fibrosis remain unknown. m⁶A methylation is closely associated with autoimmune diseases such as psoriasis [36], systemic lupus erythematosus [37,38] and rheumatoid arthritis [39], but little is known about how it affects immune abnormalities in scleroderma. Based on these two perspectives, we are interested in m⁶A methylation in scleroderma.

Although human scleroderma samples can be obtained from pathological examination, there may be significant differences between samples based on the degree of the disease and the site of skin lesions; therefore a reliable animal model to limit confounders is necessary. We selected the BLM-induced scleroderma model [40] and performed m⁶A epitranscriptomic microarray for analyzing the m⁶A methylation profile in a mouse model of scleroderma. Compared with MeRIP-seq, m⁶A epitranscriptomic microarray can estimate the levels of different types of RNAs, such as mRNA, lncRNA and circRNA. Also, MeRIP-seq requires 120 μg of total RNA, which may make it unsuitable for valuable clinical samples. Our study will provide a basis for further research on how m⁶A methylation regulates the occurrence and development of scleroderma, and provide potential therapeutic directions.

Materials & methods

Establishment of a mouse model of scleroderma

Female Balb/c mice aged 6–8 weeks were acquired from Charles River Laboratories (Shanghai, China). Sixteen mice were randomly divided into two groups: the control and a BLM-induced (BLM) group. In the BLM group, 100 μl BLM (1 mg/ml; Thermo Fisher Scientific, MA, USA) was injected subcutaneously into the back of the mice daily for 28 consecutive days. Similarly, mice in the control group were injected with 100 μl phosphate-buffered saline daily for 28 days. In the BLM group, a circular depilation area with a diameter of 0.1 cm was observed at the injection site, along with scar formation. Compared with the control group mice, the BLM group mice gained less weight (Figure 1F).

Figure 1. A mouse model of scleroderma induced by bleomycin was established.

(A) Skin specimens were stained with hematoxylin and eosin (top) and Masson’s trichrome (bottom). Original magnification ×100. Scale bar = 100 μm. (B) Dermal thickness and (C) collagen content were quantified from the hematoxylin and eosin and Masson’s trichrome staining data, respectively (n = 8 per group). (D)Col1α1 and (E)Col3α1 mRNA levels were examined by quantitative PCR (n = 4 per group). (F) Changes in body weight of mice in the control group and BLM-treated group within 4 weeks. Mean ± standard deviation.

*p < 0.05; **p < 0.01; ****p < 0.001 compared with the control group.

^#No significance compared with the control group.

BLM: Bleomycin.

Histological analysis

All mice were sacrificed by cervical dislocation to allow for histological analysis. The skin samples were obtained from the lower right of the back, where BLM was injected. Skin samples (n = 8, each group) were fixed with 4% paraformaldehyde for 24 h at 4°C and then paraffin-embedded. The skin samples were stained with hematoxylin and eosin and Masson’s trichrome stain to assess the dermal thickness and collagen content, respectively. ‘Dermal thickness’ refers to the distance between the epidermal–dermal junction and the dermal–subcutaneous fat junction. ‘Collagen content’ refers to the ratio of collagen deposition area to the tissue area. Sections were observed using a light microscope (Olympus, Tokyo, Japan) and photographed using a digital camera. The evaluation was conducted by two independent examiners using ImageJ software (NIH, MD, USA), and three areas in each section were randomly selected.

Quantitative real-time PCR

Total RNA extraction was performed using TRIzol™ reagent (Sigma Aldrich, MO, USA), following the manufacturer’s instructions. cDNA was reverse transcribed using SuperScript™ III Reverse Transcriptase kit (Invitrogen, CA, USA), and quantitative PCR (qPCR) was performed with 2X PCR master mix (Arraystar, MD, USA) on the ViiA 7 Real-time PCR System (Thermo Fisher Scientific). The mRNA levels of target genes were normalized to the mRNA levels of ACTB or GAPDH via the 2^-ΔΔCt method. The primers used are presented in Supplementary Table 1. We conducted qPCR experiments twice; the first time, Col1α1 and Col3α1 were normalized to ACTB (n = 4, each group), and the second time, the other indicators were normalized to GAPDH (n = 4, each group).

RNA extraction & quality control

Skin tissue was cut into small pieces and homogenized in TRIzol reagent to isolate total RNA. The quantification and quality of RNA were examined using NanoDrop ND-1000 (Thermo Fisher Scientific, Shanghai, China). The extracted total RNA was stored at -80°C for the subsequent experiments. The quantification and quality of RNA are shown in Supplementary Figure 1 & Supplementary Table 2.

m⁶A immunoprecipitation

The extracted RNA was immunoprecipitated with an anti-m⁶A antibody according to the manufacturer’s protocol. For this, 1–3 μg RNA and m⁶A spike-in control mixture was added to 300 μl 1× immunoprecipitation (IP) buffer containing 2 μg anti-m⁶A rabbit polyclonal antibody (Synaptic Systems, Göttingen, Germany). The reaction was incubated with head-over-tail rotation at 4°C for 2 h. Following this, 20 μl of Dynabeads™ M-280 sheep anti-rabbit IgG suspension (Thermo Fisher Scientific) per sample was blocked with freshly prepared 0.5% bovine serum ALB at 4°C for 2 h, washed thrice with 300 μl 1× IP buffer and subsequently resuspended in the prepared RNA–antibody mixture. RNA binding to the m⁶A antibody beads was carried out with head-over-tail rotation at 4°C for 2 h. The beads were then washed thrice with 500 μl 1× IP buffer and twice with 500 μl wash buffer. Finally, the enriched RNA was eluted with 200 μl elution buffer at 50°C for 1 h. The RNA was extracted by the acid phenol–chloroform and ethanol precipitation method.

Labeling & hybridization

An equal amount of calibration spike-in control RNA was added to the immunoprecipitated RNA (IP) and supernatant RNA (Sup), following which they were amplified and labeled with Cy3 (for Sup) and Cy5 (for IP) using Arraystar Super RNA Labeling Kit (Arraystar). The synthesized cRNAs were purified using RNeasy Mini Kit (Sigma Aldrich). The cRNA concentration and specific activity were measured using NanoDrop™ ND-1000 (Thermo Fisher Scientific; Supplementary Table 3).

Next, 2.5 μg of the Cy3- and Cy5-labeled cRNAs were mixed. The cRNA mixture was fragmented by adding 5 μl 10× blocking agent and 1 μl of 25× fragmentation buffer, heating at 60°C for 30 min and combining with 25 μl 2× hybridization buffer. Next, 50 μl hybridization solution was dispensed into the gasket slide and assembled to the m⁶A epitranscriptomic microarray slide. The slides were incubated at 65°C for 17 h in an Agilent hybridization oven (Agilent, CA, USA). Finally, the hybridized arrays were washed, fixed and scanned using an Agilent G2505C scanner (Agilent).

Epitranscriptomic microarray data analysis

The array images were analyzed using Agilent Feature Extraction software (Agilent). The raw intensities of IP and Sup were normalized with an average of log2-scaled spike-in RNA intensities. Following normalization, probe signals with Present or Marginal QC flags in at least four out of eight samples were retained for further m⁶A methylation level, m⁶A quantity and expression level analyses. For determining the m⁶A methylation level, the percentage modification was calculated based on the normalized intensities of IP and Sup. For determining the m⁶A quantity, the m⁶A methylation amount was calculated based on normalized intensities of IP. The expression level was calculated based on the total normalized intensities of IP and Sup. An additional quantile normalization method in the limma package was used to normalize the RNA expression level between arrays before probe flag screening.

Filtering with fold change ≥1.5 and p < 0.05 (unadjusted p-values) was used to identify differentially m⁶A-methylated RNAs and differentially expressed RNAs between the two comparison groups.

MeRIP-qPCR

IP (n = 4, each group) was used to validate the m⁶A epitranscriptomic microarray data. Then MeRIP-qPCR was used to quantify the RNA enrichment via 2^-ΔΔct analysis. The primers used are described in Supplementary Table 1.

Statistical analysis

All data were presented as mean ± standard deviation. A two-sided t-test was used to compare the two unpaired groups using GraphPad Prism 8.0 (GraphPad, CA, USA). A p-value of <0.05 was considered statistically significant.

Results

Progression of Balb/c mouse skin fibrosis induced by BLM

After the subcutaneous injection of BLM for 4 weeks, a mouse model of scleroderma was successfully established. Compared with the control group, hematoxylin and eosin staining (Figure 1A) showed that the distinction between epidermal, dermal and subcutaneous tissues was not clear in the BLM group. Further, dermal thickness had increased in the BLM group (Figure 1B). In the BLM group, Masson’s trichrome staining (Figure 1A) showed that excessive collagen was deposited in the dermis, which further extended to the subcutaneous fat and even substituted fat tissues. Collagen content analysis also revealed the same result (Figure 1C). To evaluate the changes in extracellular matrix (ECM) production, we selected Col1α1 and Col3α1 as indicators, because they are mainly distributed in the skin. Compared with the control group, the mRNA expression of Col1α1 and Col3α1 was significantly increased in the BLM group (Figure 1D & E). In summary, a mouse model of scleroderma was successfully established, making it possible for the subsequent experiments to succeed.

Alteration of m⁶A-methylated mRNAs in the BLM-induced scleroderma

The results of three-dimensional principal components analysis revealed distinct clustering and segregation between the control and BLM groups. The percentages of variance associated with principal component (PC) 1, PC2 and PC3 were 43.75, 32.94 and 23.77%, respectively (Figure 2A). Moreover, a corresponding MeRIP enrichment %(MeRIP/Input) of positive and negative spike-in controls was present between the control and the BLM groups, indicating a rigorous analysis of the m⁶A epitranscriptomic microarray (Figure 2B). Hierarchical clustering analysis revealed the intragroup similarities and intergroup differences in methylated mRNA patterns (Figure 2C). A volcano plot revealed that compared with the control group, 843 mRNAs in the BLM group had significantly different m⁶A modification, including 733 hypermethylated mRNAs and 110 hypomethylated mRNAs, all with fold change ≥1.5 and p < 0.05 (Figure 2D). The differentially m⁶A-methylated mRNAs were mapped to all chromosomes except Y (Figure 2E), because the Balb/c mice were female. The top five chromosomes with the most hypermethylated and hypomethylated m⁶A mRNAs (Figure 2F) were: chromosomes 2 (63), 4 (54), 7 (73), 9 (79), 11 (48); and 1 (8), 3 (8), 6 (8), 7 (8), 11 (10), 13 (9) and 15 (10), respectively.

Figure 2. The change of N⁶-methyladenosine-methylated mRNAs in the bleomycin-induced mouse model of scleroderma.

(A) PC analysis was conducted to classify the differentially methylated mRNAs of the control group and the BLM group. (B) Amplification of positive control primer set (Pos) and negative control primer set (Neg) by quantitative PCR assay indicated comparable m⁶A-RNA-immunoprecipitation enrichment % (m⁶A-RNA-immunoprecipitation/Input) between the control group and the BLM group. (C) The heat map shows opposite m⁶A-methylated mRNAs between the control group and the BLM group. (D) The volcano plot shows the differently m⁶A-methylated mRNAs between the control and BLM groups with fold change ≥1.5 and p < 0.05. (E) Chromosomal distribution of m⁶A-methylated mRNAs. n = 4, each group. Mean ± standard deviation.

^#No significance compared with the control group.

BLM: Bleomycin; Chr: Chromosome; m⁶A: N⁶-methyladenosine; Neg: Negative; PC: Principal componant; Pos: Positive.

Association analysis of m⁶A methylation & mRNA expression in the BLM-induced scleroderma

The microarray data revealed significant differences in the global mRNA expression between the control and BLM groups. As shown in Supplementary Figure 2, the heat map and volcano plot showed 1991 mRNAs with increased expression and 1993 mRNAs with decreased expression (fold change ≥1.5 and p < 0.05).

The differentially expressed mRNAs were analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses to explore the general biological functional pathways in scleroderma of the mRNAs. GO analysis showed that the differentially upregulated mRNAs were significantly associated with immune system process (biological process), cellular anatomical entity (cellular component) and protein binding (molecular function; Supplementary Figure 3A). The differentially downregulated mRNAs were significantly associated with receptor cellular process (biological process), intracellular (cellular component) and binding (molecular function; Supplementary Figure 3B). As shown in Supplementary Figure 3C, the differentially upregulated mRNAs were most significantly associated with malaria, cell adhesion molecules, Kaposi sarcoma-associated herpes virus infection and cytokine–cytokine receptor interaction. The differentially downregulated mRNAs were most significantly associated with tight junction, Hippo signaling pathway, melanogenesis and basal cell carcinoma (Supplementary Figure 3D).

Compared with markers for skin fibrosis, the expression of markers for vascular damage and inflammation was increased in all tested samples, indicating a successful mouse model (Table 1). The top 20 differentially expressed mRNAs are listed in Table 2.

Table 1. Specific gene expression associated with scleroderma.

Gene	Fold change	p-value	Chromosome	Association
Tgfb3	3.5257147	2.35E-03	12	The most important signal transduction pathway
Tgfb1	1.7898921	4.99E-03	7	The most important signal transduction pathway
Col2a1	3.1446831	1.89E-03	15	Skin fibrosis (collagen deposition)
Col1a1	2.9964116	1.91E-02	11	Skin fibrosis (collagen deposition)
Col3a1	2.5785271	2.93E-02	1	Skin fibrosis (collagen deposition)
Acta2	3.6369475	1.85E-03	19	Skin fibrosis (epithelial–mesenchymal transition)
Vim	2.4594609	2.82E-02	2	Skin fibrosis (epithelial–mesenchymal transition)
Pecam1	1.688203	4.03E-02	11	Vascular damage
Vcam1	5.104913	2.63E-05	3	Vascular damage
Il1b	4.037706	2.04E-03	2	Inflammatory reaction
Il6	9.747919	1.88E-04	5	Inflammatory reaction
Il33	1.8198779	3.26E-02	19	Inflammatory reaction

Table 2. The top 20 differentially expressed mRNAs in the mouse model of scleroderma.

Gene	Regulation	Fold change	p-value	Chromosome
Clec4d	Up	130.39764	1.04E-06	6
Mmp13	Up	98.831288	5.35E-06	9
Sprr2d	Up	88.734418	2.72E-06	3
Sprr2e	Up	86.931675	5.94E-05	3
Ccl3	Up	81.853399	1.67E-06	11
Ccl4	Up	76.467386	4.66E-06	11
Wfdc12	Up	71.755784	1.30E-06	2
Cenpn	Up	49.673087	1.14E-05	8
Ankrd1	Up	48.647211	1.09E-04	19
Il1b	Up	46.46236	1.88E-05	2
Lce3f	Up	41.530177	1.51E-05	3
Nmu	Up	37.086461	3.82E-07	5
Myh7	Up	27.048769	1.71E-02	14
Spink12	Down	0.060961	4.51E-04	18
Krt72	Down	0.0627622	3.79E-05	15
Inpp5a	Down	0.0712282	3.53E-02	7
Asphd2	Down	0.0741685	3.49E-02	5
Tmc8	Down	0.0828253	1.14E-03	11
Txndc5	Down	0.0828767	1.12E-04	13
Crisp3	Down	0.0828781	1.22E-04	17

The association analysis of m⁶A methylation and mRNA expression with fold change ≥1.5 revealed that the BLM group had 37 m⁶A-hypermethylated and upregulated mRNAs, 1164 m⁶A-hypermethylated and downregulated mRNAs, 27 m⁶A-hypomethylated and downregulated mRNAs and 358 m⁶A-hypomethylated and upregulated mRNAs (Figure 3A). A study indicated that m⁶A methylation could affect the stability of mRNA [41]. Following the steps described in that study, we found more hyper/down and hypo/up mRNAs than hyper/up and hypo/down mRNAs, suggesting that m⁶A methylation was negatively correlated with mRNA expression in the mouse model of scleroderma.

Figure 3. Association analysis of N⁶-methyladenosine-methylated and expressed mRNAs in the bleomycin-induced mouse model of scleroderma.

(A) The four-quadrant and Venn diagram showing the N⁶-methyladenosine-methylated and expressed mRNAs (fold change ≥1.5). (B) The protein–protein interaction network showed the connection between mRNAs in the combined analysis (fold change ≥2).

BLM: Bleomycin.

The protein–protein interaction network was established to show the connection between mRNAs (fold change ≥2) in the combined analysis (Figure 3B).

Differentially m⁶A-methylated mRNAs were involved in important pathways

GO analysis showed that the differentially m⁶A-hypermethylated mRNAs were significantly associated with binding (biological process), cellular anatomical entity (cellular component) and cellular process (molecular function; Figure 4A). In addition, the differentially m⁶A-hypomethylated mRNAs were significantly associated with receptor regulator activity (biological process), extracellular region (cellular component) and multicellular organismal process (molecular function; Figure 4B). As shown in Figure 4C, the differentially m⁶A-hypermethylated mRNAs were most significantly associated with growth hormone synthesis, secretion and action, insulin secretion and amphetamine addiction. In addition, the differentially m⁶A-hypomethylated mRNAs were most significantly associated with rheumatoid arthritis, Toll-like receptor signaling pathway and amoebiasis (Figure 4D).

Figure 4. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analysis of differentially N⁶-methyladenosine-methylated mRNAs.

(A) The top ten Gene Ontology terms of hypermethylated mRNAs, ranked by enrichment score. (B) The top ten Gene Ontology terms of hypomethylated mRNAs, ranked by enrichment score. (C) The top ten significantly enriched pathways of hypermethylated mRNAs, ranked by enrichment score. (D) The top ten significantly enriched pathways of hypomethylated mRNAs, ranked by enrichment score.

Validation of the m⁶A epitranscriptomic microarray by MeRIP-qPCR & quantitative real-time PCR

To validate the m⁶A epitranscriptomic microarray data, we conducted MeRIP-qPCR for four hypermethylated mRNAs (Hras, Lama3, Stat3 and Tnc) and four hypomethylated mRNAs (Ccl3, Ccl9, Saa1 and Il1b). These have been preliminarily studied in fibrosis, but not deeply (Table 3).

Table 3. Review for the selected genes in scleroderma.

Gene	Chromosome	Regulation	Regulatory roles on fibrosis	Ref.
Hras	7	Down	Accelerating fibroblast proliferation rates	[42]
Lama3	18	Up	Destroying the balance of epithelial–mesenchymal transition	[43]
Tnc	4	Up	Inducing collagen deposition and myofibroblast differentiation; a leading-edge gene of the PI3K-Akt pathway in skin fibrosis	[44–46]
Stat3	11	Up	A transcription factor activated downstream of TGF-β	[47,48]
Ccl3	11	Up	A critical chemokine for fibrosis and inflammation in scleroderma	[49]
Ccl9	11	Up	Mediates macrophage regulation of extracellular matrix	[50]
Saa1	7	Up	Promotes macrophage-to-myofibroblast transition	[51]
Il1b	2	Up	An important proinflammatory factor in scleroderma	[52]

We selected hypermethylated mRNAs from the perspective of KEGG analyses. Because collagen deposition is the most important pathological process in scleroderma, mRNAs involved in the ECM–receptor interaction pathway were analyzed. We sorted fold change from high to low, and then picked two mRNAs (Table 4). Although COL4A1 ranked the second highest, we did not select this mRNA because it has received little attention in skin fibrosis compared with another two collagen proteins, COL1A1 and COL3A1, which have been intensively studied. Instead, we selected novel LAMA3 and TNC, both of which play relatively clear roles in skin fibrosis. As the pathways with the highest enrichment score, growth hormone synthesis, secretion and action were also selected. These mRNAs included Pik3cd, Cacna1s, Pik3r2, Ghrl, Atf6b and Atf2 (Table 5). Except for Pik3cd [53] and Atf2 [54,55], the other four mRNAs are not involved in skin fibrosis. A report linking Pik3r2 with scleroderma [56] has been published since the completion of our study. In contrast, there have been many studies on Stat3 in scleroderma, and there are two publications on DNA methylation [57,58]; hence we chose Stat3 for verification.

Table 4. Differentially hypermethylated mRNAs involved in the extracellular matrix–receptor interaction pathway.

Transcript_ID	Gene	Fold change	p-value
ENSMUST00000030056	Tnc	2.906998354	0.001655676
ENSMUST00000208095	Col4a1	2.752693479	0.00095491
uc008ecg.1	Lama3	1.935457924	0.00273348
uc008kia.1	Itgav	1.900846576	0.019797529
ENSMUST00000166905	Dag1	1.726454757	0.039187672
ENSMUST00000042729	Npnt	1.584309017	0.005438113
ENSMUST00000191899	Dag1	1.466675764	0.017589962
ENSMUST00000141725	Itgav	1.226947816	0.016770676

Table 5. Differentially hypermethylated mRNAs involved in the growth hormone synthesis, secretion and action pathway.

Transcript_ID	Gene	Fold change	p-value
ENSMUST00000026572	Hras	2.113980372	0.032686597
ENSMUST00000105689	Pik3cd	1.906831767	0.040208284
ENSMUST00000112064	Cacna1s	1.803270257	0.041289197
ENSMUST00000154685	Pik3r2	1.725025645	0.049359786
ENSMUST00000204163	Ghrl	1.719122443	0.020129527
uc012aqj.2	Atf6b	1.710327299	0.046871714
ENSMUST00000112016	Atf2	1.677316698	0.009778686
ENSMUST00000127638	Stat3	1.587988361	0.041871594
ENSMUST00000020215	Socs2	1.536058126	0.009307881
NM_001256001	Cacna1c	1.528964604	0.003969553
ENSMUST00000205120	Creb5	1.52871572	0.028358205
ENSMUST00000112793	Cacna1c	1.52271918	0.014032103

We selected hypomethylated mRNAs from the perspective of GO/KEGG analyses. Scleroderma is an autoimmune disease, so we selected the duplicated mRNAs in cytokine–cytokine receptor interaction pathways and the inflammatory response (biological process), namely Il1f6, Il1b, Ccl3 and Ccl9. Although Il1f6 ranked the lowest, there are no publications related to its association with fibrosis. Thus we chose Saa1 because its fold change was second only to that of Ccl3 (Tables 6 & 7).

Table 6. Differentially hypomethylated mRNAs involved in cytokine–cytokine receptor interaction pathways.

Transcript_ID	Gene	Fold change	p-value
ENSMUST00000028361	Il1f6	0.331783132	0.00016362
uc008mhu.2	Il1b	0.474644948	0.01262423
ENSMUST00000001008	Ccl3	0.488395453	0.017125028
ENSMUST00000108648	Tnfsf13	0.549612429	0.005544923
ENSMUST00000019266	Ccl9	0.579709293	0.041642398

Table 7. Differentially hypomethylated mRNAs involved in inflammatory response (biological process).

Transcript_ID	Gene	Fold change	p-value
ENSMUST00000028361	Il1f6	0.331783132	0.00016362
uc008mhu.2	Il1b	0.474644948	0.01262423
ENSMUST00000001008	Ccl3	0.488395453	0.017125028
ENSMUST00000128088	Saa1	0.499300978	0.030405502
ENSMUST00000069927	S100a8	0.499945891	0.020490792
ENSMUST00000088041	Ptgis	0.578128978	0.039478702
ENSMUST00000019266	Ccl9	0.579709293	0.041642398
ENSMUST00000026881	Ly96	0.593689188	0.033149311

As shown in Figure 5A, the results of MeRIP-qPCR revealed a significant methylation increase in Hras, Lama3 and Tnc (but not Stat3: p = 0.052) and a significant methylation decrease in Ccl3, Ccl9, Saa1 and Il1b in the BLM group compared with the control group. As shown in Figure 5B, these results showed an identical trend of m⁶A-modified mRNA changes between microarray and MeRIP-qPCR (except for Lama3 and Stat3). Fold change in Lama3 in the BLM group compared with the control group was 1.32, which was less than 1.5.

Figure 5. N⁶-methyladenosine enrichment and mRNA expression in the candidate mRNAs.

(A) The relative m⁶A level of the selected mRNAs was detected using MeRIP-qPCR. (B) The comparison between microarray data and MeRIP-qPCR of the selected mRNAs. (C) The six mRNAs were analyzed using qPCR. (D) The comparison between microarray data and qPCR of the selected mRNAs. (E) Identification of RNA binding proteins to the validated m⁶A-methylated mRNAs differentially. n = 4, each group. Mean ± standard deviation.

^#No significance.

*p < 0.05; **p < 0.01; ***p < 0.001 compared with the control group.

BLM: Bleomycin; m⁶A: N⁶-methyladenosine; MeRIP: m⁶A-RNA-immunoprecipitation; qPCR: Quantitative PCR.

Next, we analyzed the mRNA expression of Hras, Tnc, Ccl3, Ccl9, Saa1 and Ilb1 to determine whether mRNA expression was regulated by m⁶A methylation. However, qPCR showed no statistically significant differences in the expression levels of Tnc and Ccl9 between the BLM and control groups (Figure 5C). Compared with the control group, the m⁶A modification and mRNA expression of Hras, Ccl3, Saa1 and Ilb1 revealed the opposite trend in the BLM group. The results also showed identical trends of mRNA expression between microarray and qPCR results (Figure 5D).

Finally, we used RMBase 2.0 to explore the potential RNA-binding proteins bound to the validated m⁶A-methylated mRNAs [59]. As shown in Figure 5E, Hras, Tnc, Ccl3 and Ccl9 had corresponding RNA-binding proteins. Among them, Mbnl1, Mbnl3 and Taf15 transcripts were also included in our microarray data, and their mRNA expression levels showed statistically significant differences between the BLM and control groups.

Discussion

Scleroderma is an autoimmune disease with a poor prognosis whose pathogenesis remains unclear. Compared with 0.03% incidence in the general population, the incidence of scleroderma in families with a history of scleroderma is 1.5% [60], revealing that heredity is a significant risk factor. Several studies have shown that scleroderma is closely related to epigenetics, and the main research directions involve the investigation of post-translational modifications, DNA methylation, mRNAs and lncRNAs [61]. However, the RNA m⁶A modification landscape in scleroderma remains largely unknown. This study elucidated the m⁶A landscape in the BLM-induced mouse model of scleroderma. First, we demonstrated that BLM induced skin collagen deposition and inflammation in mice. Next, we performed m⁶A epitranscriptomic microarray studies to identify differentially methylated and expressed mRNAs. GO and KEGG pathway analyses of differentially m⁶A-methylated genes identified their potential roles in physiological and pathological mechanisms. Subsequently, we performed an association analysis of m⁶A methylation and mRNA expression. Finally, we verified the candidate mRNAs via MeRIP-qPCR and qPCR. The selected methyltransferases and demethylases were verified using qPCR. The above findings better our understanding of the role of RNA m⁶A methylation in scleroderma.

KEGG analysis revealed that the differentially hypermethylated m⁶A mRNAs were enriched in pathways involved in diabetes, such as insulin secretion, insulin resistance and Type II diabetes mellitus. A nationwide cohort study in Taiwan found that the incidence of diabetes was 10.02% in 2671 scleroderma patients, compared to 24.81% in 7769 controls [62]. A similar study in Minnesota, USA reported that the incidence of diabetes was 2.56% in 78 scleroderma patients, compared to 12.18% in 125 non-scleroderma patients [63]. These two studies revealed epidemiologically a lower prevalence of diabetes in scleroderma patients. Some studies have identified that an increase of IL-10 [64,65] and IL-13 [66,67] led to immune system abnormalities in scleroderma, which was a protective factor for diabetes. TNF-related apoptosis-inducing ligand (TRAIL) is decreased in diabetic patients, and the inhibition of TRAIL intensifies the apoptosis of pancreatic β cells; TRAIL is elevated in scleroderma and might be beneficial for scleroderma [68]. Metformin, a routine treatment for diabetes, could also reduce skin fibrosis and inhibit inflammatory cytokines such as IL-1β, IL-6, IL-17 and TNF-α [69]. A new technique of metabolomic profiling has found that the excessive deposition of ECM depends on metabolic disturbances [70]. Hesperidin 165 (NCT03377140) [71] and Omacor 166 (NCT03018041) [72] have been identified as candidate antifibrotic drugs given their effects in regulating glucose metabolism. As differential mRNAs were enriched in the pathways of diabetes, we speculate that m⁶A methylation is related to the above-mentioned relationship between diabetes and scleroderma.

Unlike the hypermethylated mRNAs, the differentially hypomethylated m⁶A mRNAs were enriched in the pathways of immune regulation, such as rheumatoid arthritis, Toll-like receptor signaling pathway, cytokine–cytokine receptor interaction and graft-versus-host disease, revealing that m⁶A methylation plays a role in immunoregulatory alteration. For example, Saa1 plays the following key role in inflammation and fibrosis: Saa1 knockdown can suppress the process of macrophage-to-myofibroblast transition, which can inhibit renal interstitial fibrosis [51]. Elevated SAA1 is associated with a fibrotic pulmonary phenotype [73]. Different plasma levels of SAA1 were related to different grades of inflammation and stages of fibrosis in the livers of patients with autoimmune hepatitis [74]. In our study, Saa1 was enriched in inflammatory response (biological process) and was highly expressed in skin fibrosis, similar to the above findings, and was differentially hypermethylated, indicating that methylation may be the cause of its role in fibrosis.

As shown in Supplementary Figure 4, a KEGG network showed some critical pathways in scleroderma, which are mainly involved in immune response and cell proliferation. Further studies are warranted to explain the exact pathogenesis.

m⁶A methylation is a reversible process, dynamically regulated by key m⁶A enzymes. However, in the present study the expression levels of methyltransferases (writers), demethylases (erasers) and binding proteins (readers) were unaltered between the control group and BLM group (Figure 6A). There were no statistically significant differences in the expression levels of Mettl3 and Alkbh5 between the BLM and control groups (Figure 6B). Our sample size may be too small, warranting further studies with larger sample sizes to explore the reasons underlying m⁶A methylation alterations in the mouse model of scleroderma. In addition, no other studies have investigated the involvement of these key enzymes in scleroderma, making comparison of our results difficult. Therefore functional experiments involving the knockdown/overexpression of key enzymes are needed to confirm the regulatory role of m⁶A methylation in scleroderma.

Figure 6. Analysis of N⁶-methyladenosine methylation regulators in the bleomycin-induced mouse model of scleroderma.

(A) The relative expression of writers, erasers and readers in microarray data. (B)Mettl3 and Alkbh5 were analyzed using quantitative PCR. n = 4, each group. Mean ± standard deviation.

BLM: Bleomycin.

Conclusion

Our study reports on differentially m⁶A-methylated mRNAs in a BLM-induced mouse model of scleroderma, demonstrating a potential strong connection between RNA m⁶A methylation and the regulation of fibrosis and inflammation. This study hence provides novel insights into the pathogenesis of scleroderma.

Future perspective

Our study revealed alterations in Hras, Lama3, Stat3, Tnc, Ccl3, Ccl9, Saa1 and Il1b, which require further investigation. These mRNAs may play a critical role in epigenetic mechanisms, contributing to the identification of diagnostic markers and therapeutic targets for scleroderma.

Summary points

• A total of differentially 733 hypermethylated N⁶-methyladenosine (m⁶A) mRNAs and 110 hypomethylated m⁶A mRNAs were identified in the bleomycin-induced mouse model of scleroderma.

• The hypermethylated m⁶A mRNAs were enriched in pathways of diabetes, such as insulin secretion, insulin resistance and Type II diabetes mellitus and growth hormone synthesis, secretion and action.

• The hypomethylated m⁶A mRNAs were enriched in pathways of immunoregulation, such as rheumatoid arthritis, Toll-like receptor signaling pathway, cytokine–cytokine receptor interaction and graft-versus-host disease.

• m⁶A-RNA-immunoprecipitation quantitative PCR (qPCR) verified that Hras and Tnc were hypermethylated m⁶A mRNAs, while Ccl3, Ccl9, Saa1 and Il1b were hypomethylated m⁶A mRNAs. qPCR verified that Hras was upregulated, while Ccl3, Saa1 and Ilb1 were downregulated.

• Using RMBase 2.0, MBNL1, MBN3 and TAF15 were predicted to be RNA-binding proteins of Hras and Tnc.

• There were no statistically significant differences in the expression levels of methyltransferases, demethylases or binding proteins between the control and bleomycin-treated groups.

Ethical conduct of research

The authors state that they have obtained appropriate institutional review board approval or have followed the principles outlined in the Declaration of Helsinki for all animal experimental investigations, and all experiments were approved and performed following the guidelines of the Tongji University Ethic committee of China (no. TJ-HB-LAC-2020-83).

Supplementary data

To view the supplementary data that accompany this paper please visit the journal website at: www.tandfonline.com/doi/suppl/10.2217/epi-2021-0369

Financial & competing interests disclosure

This work is supported by the grants from the National Natural Science Foundation of China, which is a subsidiary of the Ministry of Science and Technology, PRC. The General Program number is 81874240. The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed.

Writing assistance was provided by Enago and was funded by the National Natural Science Foundation of China (program no. 81874240).

Additional information

Funding

This work is supported by the grants from the National Natural Science Foundation of China, which is a subsidiary of the Ministry of Science and Technology, PRC. The General Program number is 81874240. The authors have no other relevant affiliations or financial involvement with any organization or entity with a financial interest in or financial conflict with the subject matter or materials discussed in the manuscript apart from those disclosed. Writing assistance was provided by Enago and was funded by the National Natural Science Foundation of China (program no. 81874240).