N6-methyladenosine-methylomic landscape of lung tissues of mice with chronic obstructive pulmonary disease

Chronic obstructive pulmonary disease (COPD), a common respiratory disease, can be divided into stable phase and acute exacerbation phase (AECOPD) and is characterized by inflammation and hyper-immunity. Methylation of N6-methyladenosine (m6A) is an epigenetic modification that regulates the expression and functions of genes by influencing post-transcriptional RNA modifications. Its influence on the immune regulation mechanism has attracted great attention. Herein, we present the m6Amethylomic landscape and observe how the methylation of m6A participates in the pathological process of COPD. The m6A modification of 430 genes increased and that of 3995 genes decreased in the lung tissues of mice with stable COPD. The lung tissues of mice with AECOPD exhibited 740 genes with hypermethylated m6A peak and 1373 genes with low m6A peak. These differentially methylated genes participated in signaling pathways related to immune functions. To further clarify the expression levels of differentially methylated genes, RNA immunoprecipitation sequencing (MeRIP-seq) and RNA-sequencing data were jointly analyzed. In the stable COPD group, 119 hypermethylated mRNAs (82 upregulated and 37 downregulated mRNAs) and 867 hypomethylated mRNAs (419 upregulated and 448 downregulated mRNAs) were differentially expressed. In the AECOPD group, 87 hypermethylated mRNAs (71 upregulated and 16 downregulated mRNAs) and 358 hypomethylated mRNAs (115 upregulated and 243 downregulated mRNAs) showed differential expression. Many mRNAs were related to immune function and inflammation. Together, this study provides important evidence on the role of RNA methylation of m6A in COPD.

1 Introduction

Chronic obstructive pulmonary disease (COPD) is a common respiratory disease that causes persistent respiratory symptoms and airflow restriction, owing to airway and/or alveolar abnormalities (1). The World Health Organization survey in 2019 reported 380 million COPD cases in the world and more than 3 million patients annually die from this pathology, which is very alarming (2, 3). The pathogenesis of COPD involves inflammation orchestrated by various inflammatory cells and chemokines. According to clinical symptoms, COPD can be divided into a stable phase and an acute exacerbation phase (AECOPD) (4).

Emerging studies have revealed the role of epigenetic regulation in inflammation, especially the modification of RNA methylation regulatory factors. N6-methyladenosine (m6A) methylation, the most common and abundant post-transcriptional modification of eukaryotic RNA, is dynamically reversible and regulated by methyltransferase, demethylase, and methylated reading protein in eukaryotic cells (5). The methylation process of m6A is mainly regulated by methyltransferases METTL3 and METTL14 that modify m6A by catalyzing adenylate on the mRNA. The demethylase FTO can demethylate the m6A-modified base. The methylation reading proteins YTHDF and YTHDC activate the downstream regulatory pathway by recognizing the base where m6A occurs (6). It was found that m6A was specifically modified at the RRACH site where R is adenine and guanine, H is adenine, cytosine, and uracil, and A is an m6A sequence (7–9). Meanwhile, m6A modification is also related to inflammation in alcoholic kidney injury, coronary artery disease, spontaneous colitis, and other diseases (10–12). A few studies have analyzed the relationship between m6A methylation and COPD. To study the differences in the m6A modification mode between a healthy and COPD model and further explore the relationship between m6A modification and lung tissue inflammation disorder, herein we determine differentially methylated mRNA transcripts in COPD and normal mouse lungs by methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq). In addition, we analyze the expression of differentially methylated mRNAs in combination with MeRIP-seq and RNA-seq data to investigate if mRNA methylation affects the corresponding gene expression. Overall, this study reveals the potential mechanism of m6A methylation in COPD.

2 Materials and methods

2.1 Animals and tissue collection

Wild-type C57BL/6J male mice (6-8 weeks old) were purchased from Animal Center of Xinjiang Medical University (China). Mice were housed in a specific pathogen-free environment at Xinjiang Medical University Animal Experimentation Facility under a 12h light/12h dark cycle at 23 ± 1°C and 50 ± 10% relative humidity, fed a non-purified diet and drinking water ad libitum. The study was approved by the Institutional Review Board and Animal Experimentation Committee of Xinjiang Medical University (license number: SYXK-2018-0003) and conducted in accordance with the Xinjiang Medical University Guide for the Care and Use of Laboratory Animals.

Animals were randomly divided into 3 groups: control group (CTL, inhalation of cigarette smoke free room air [CS] exposure, n=10), CS exposure group (CS, n=10), CS exposure with lipopolysaccharide (LPS) infusion group (CS + LPS, n=10). Mice in the CS and CS+LPS groups were placed in a 60 x 40 x 30 cm³ fume box for 2h per exposure, twice daily, 6 days per week, and a mouse COPD model was established after 90 days of CS exposure (13). The AECOPD mouse model was established on day 91 by injection of LPS (1 μg/each, containing 50 μL of saline) (14). Equal amounts of saline were injected into the CS-exposed mice and the control group. Lung function, body weight and neutrophil counts in alveolar lavage fluid were assessed in all animals prior to execution on day 92. The cigarettes used in the mouse model establishment were commercial Honghe brand cigarettes (Hongyun Honghe Tobacco Industry, China). Each cigarette yields 10 mg tar, 1.0 mg nicotine and 11 mg CO. The total particulate matter concentration of 361.3 ± 49.2 μg/m³/day.

2.2 Hematoxylin and eosin (H&E) staining

Lung tissues were fixed in 4% formaldehyde for more than 24 h, paraffin-embedded, sliced and hematoxylin and eosin (H&E) stained. The sections were rinsed in running water, and then placed into acid water and ammonia water for color separation. After a few seconds, the sections were rinsed with running water for 1 h and then placed into distilled water. The slides were then immersed in 70% and 90% alcohol and dehydrated for 10 min. The slides were treated with alcohol eosin dye for 2 min, dehydrated, and sealed.

2.3 RNA extraction and quantitative polymerase chain reaction analysis

Total RNA was extracted from normal lung tissue, stable COPD lung tissue, and acute exacerbation COPD lung tissue (three biological replicates for each condition) using Trizol (Invitrogen, USA) following the manufacturer’s protocol. Total RNA was extracted using Trizol (Invitrogen, USA) following the manufacturer’s protocol. We used TaqMan^® RNA-to-CT™ 1-Step Kit from Thermo Fisher (METTL3, METTL14, FTO and YTHDF1) and PrimeScript RT Reagent Kit (MMP12, IL-1β, CXCL12, NFKBIA and CEBP-β) from Takara. Real-time quantitative PCR was performed on an ABI 7500 fast real-time PCR system (Thermo Fisher Scientific, USA) to determine the target RNA expression levels. The PCR system of TaqMan^® RNA-to-CT™ 1-Step Kit was as follows:48°C for 15 min, 95°C for 10min, followed by 40 cycles of 95°C for 15s, and 60°C for 1min. The PCR system of PrimeScript RT Reagent Kit was as follows: 95°C for 3 min, followed by 40 cycles of 60°C for the 30s, and 72°C for 10s. The relative quantity of the target gene was calculated with the 2^−△△Ct method and normalized GAPDH. The expression levels of METTL3, METTL14, FTO and YTHDF1 were assessed by TaqMan technology. Their primer sequence has not been listed because of intellectual property protection. However, other primer sequences are shown in Table 1.

TABLE 1

Table 1 Primers of reverse transcription PCR analysis for genes.

2.4 MeRIP-seq

MeRIP-seq was performed by Novogene Technology Co., Ltd. (Beijing, China). Briefly, a total of 300 µg RNAs were extracted from the lung tissue. The integrity and concentration of extracted RNAs were detected using an Agilent 2100 bioanalyzer (Agilent, USA) and simpliNano spectrophotometer (GE Healthcare), respectively. Fragmented mRNA (~100 nt) was incubated for 2h at 4°C with anti-m6A polyclonal antibody (Synaptic Systems, USA) in the immunoprecipitation experiment. Then, immunoprecipitated mRNAs or Input was used for library construction with NEBNext^® Ultra™ RNA Library Preparation Kit for Illumina (New England Biolabs, USA). The library preparations were sequenced on an Illumina Novaseq or Hiseq platform with a paired-end read length of 150 bp according to the standard protocols. The sequencing was carried out with 3 independent biological replicates. The data has been uploaded to NCBI’s BioProject: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA853736.

2.5 RNA-seq

RNA-seq was performed by Novogene Technology Co., Ltd. (Beijing, China). A total amount of 1 μg RNA per sample was used as input material for the RNA sample preparations. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in First Strand Synthesis Reaction Buffer. First-strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase. Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ ends of DNA fragments, Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 370~420 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index Primer. At last, PCR products were purified and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Novaseq platform and 150 bp paired-end reads were generated. The data has been uploaded to NCBI’s BioProject: https://www.ncbi.nlm.nih.gov/bioproject/?term=919075.

2.6 Data processing and statistical analysis

Adapters, duplicates and low quality sequences were removed from the raw reads and input samples of the IP using FASTP software to obtain clean data in fastq format. Me-RIP used BWA and RNA-seq used HISAT2 to compare the clean data to the small house mouse genome in Ensembl version 101 to obtain BAM files. Me-RIP uses exomePeak (Version: 2.16.0) software to analyze the distribution of peaks on the chromosomes and the differential genes in each group. Peak identification was performed for each experimental group with a threshold of FDR < 0.05 and Fold Enrichment > 1. Motif analysis we used the HOMER software to identify motifs on mRNA regions where m6A methylation occurred. RNA-seq using Stringtie software in results based on comparison to the genome on, the reads were spliced into transcripts and quantified, then analyzed for significant differences in expression using edgeR software. GO enrichment analysis was implemented by the GO seq R package. KEGG pathway analysis was implemented by KOBAS software (version 3.0). Correlation was analyzed using Pearson’s correlation. Peaks were visualized using IGV software and annotated with AI software. Weight, Neutrophil count and Gene expression levels of the samples were statistically analyzed using the Student t-test of GraphPad Prism 8 (GraphPad Software Inc, USA), and P < 0.05 was considered statistically significant.

3 Results

3.1 Construction of a mouse model of stable COPD and AECOPD and quantification of key m6A regulatory factors

Analyses of weight, neutrophil count, and H&E staining can determine the successful establishment of this experimental model. The results showed that from the third week, the weights of the mice from stable COPD and AECOPD groups were significantly lower than those of the mice from the control group (Figure 1A). Further, the number of neutrophils in the alveolar lavage fluid was counted and analyzed by Giemsa staining. The neutrophil counts in stable COPD and AECOPD groups were higher than the count in the control group (Figure 1B). The RT-qPCR results showed that the expression levels of MMP12 mRNA were significantly higher in the stable COPD and AECOPD groups than in the control group (Figure 1C). H&E staining revealed the thinner and fused alveolar walls and immune cell infiltration into the alveoli of the mice from the stable COPD and AECOPD groups. These results demonstrate the successful construction of stable COPD and AECOPD mouse models (Figure 1D).

FIGURE 1

Figure 1 Detection of m6A methylation regulatory factor expression levels between the control and experimental groups: (A) Body weights of mice measured every week; (B) neutrophil count in the alveolar lavage fluid; (C) Detection of the expression level of MMP12 mRNA; (D) H&E staining of the lung tissue; (E) detection of expression levels of m6A methylation regulatory factor mRNAs (*P < 0.05; **P < 0.01, ***P < 0.001).

To explore whether m6A methylation participates in the occurrence and development of COPD, the expression levels of m6A methylation regulatory factors, including METTL3, METTL14, FTO, and YTHDF1, were compared by RT-qPCR. In comparison with the control group, stable COPD and AECOPD groups showed a significant decrease in METTL3, METTL14, FTO, and YTHDF1 levels (Figure 1E). Hence, one might assume that there are differences in m6A methylation during the occurrence and development of COPD, which will be further discovered through MeRIP-Seq.

3.2 Distribution difference of m6A methylation on chromosomes during the occurrence and development of COPD

We performed MeRIP-Seq to explore the role of m6A methylation in the occurrence and development of COPD (Figure 2). With three replicates and one INPUT per group, a total of 12 libraries were sequenced with m6A in the transcriptome range. All differentially methylated N6-methyladenosine (DMM) in mRNAs and long-noncoding RNAs (lncRNAs) were mapped to the chromosome to observe their distribution (Figure 3A). Analysis of the distribution of m6A methylation on different chromosomes showed that the chromosomes with the most m6A methylation were chromosome 2 (2338 m6A methylation peaks), chromosome 7 (2306 m6A methylation peaks), and chromosome 11 (2198 m6A methylation peaks) (Figures 3B, C). As shown in Figure 3C, in comparison with the control group, the stable COPD group had a significant increase and the AECOPD group had a significant decrease in the number of m6A methylation.

FIGURE 2

Figure 2 Schematic of m6A-seq analysis in C57 mice: C57 mice were stimulated with cigarette smoke or LPS, and their lung tissues were excised and used to extract total RNA. The RNA was fragmented and m6A RNA separated according to Me-RIP analysis. Construction of m6A-seq library and sequencing for analysis.

FIGURE 3

Figure 3 Distribution of differentially methylated m6A sites: (A) Distribution of all DMM sites in mRNAs and lncRNAs on chromosomes; (B) density distribution of m6A peak along chromosomes; (C) quantification of m6A peak on each chromosome. DMM: differentially methylated N6 methyladenosine.

3.3 Overall m6A methylation modification on the transcriptome during the occurrence and development of COPD

In the R package, A total of 14274 m6A peaks were detected for the control group, including 8709 gene transcripts, and 13231 m6A peaks were found for the stable COPD group, including 8181 gene transcripts. Similarly, 13366 m6A peaks corresponding to 8444 gene transcripts were observed for the AECOPD group. In addition, 8162 peaks that corresponded to m6A modification of 6319 genes were observed in three groups. The control, stable COPD, and AECOPD groups had their own 2735, 2944, and 2196 new peaks, respectively, which reflected significant differences in the overall m6A modification trend between them (Figures 4A, B). To determine whether the m6A peak contains an RRACH conservative sequence motif, samples from the three groups were tested and further analyzed by the HOMER software. The results showed that the consistent motif of m6A RRACH was AAACC (Figure 4C).

FIGURE 4

Figure 4 m6A-Seq of the transcriptome width reveals the overall m6A modification mode in the COPD group: (A) Number of common and specific m6A peaks in the control, stable COPD, and AECOPD group; (B) Venn diagram of transcript m6A peaks for the three groups; (C) motifs enriched from m6A peaks identified from the control, stable COPD, and AECOPD groups.

3.4 Distribution characteristics of transcriptome m6A methylation peaks during the occurrence and development of COPD

The meta-genomic model of the m6A peak was analyzed to determine the differential distribution of m6A in the full transcriptome. The results showed that the m6A peak was mainly concentrated at the beginning of the 3′-untranslated region (3′-UTR) (Figure 5A). Then, we explored the number of m6A modification peaks for each gene, and found that almost 60% of the methylated genes exhibited only one m6A peak while most genes contained 1-3 m6A peaks (Figure 5B).

FIGURE 5

Figure 5 Distribution of m6A methyl groups in the transcript: (A) m6A peaks enriched with transcripts, wherein a transcript was separated into 5′-UTR, coding sequence, and 3′-UTR; (B) proportion of gene-processing m6A peaks with different numbers; (C) pie chart shows the percentage of m6A peak in the five segments of transcript. The m6A peak was most enriched in the CDS segment.

To analyze the distribution of m6A peaks in mRNAs, we divided the transcript into five transcriptional fragments as follows: Transcription start site (TSS), 5′-untranslated region (5′-UTR), coding sequence (CDS), stop-codon, and 3′-UTR. The results demonstrated that m6A was most commonly detected in CDS and was also enriched in 3′-UTR (Figure 5C).

3.5 DEGs of m6A methylation during the occurrence and development of COPD

To explore the role of m6A in the occurrence and development of COPD, we compared the abundance of m6A peak between the stable COPD and control group as well as between AECOPD and control group. The stable COPD group had 430 upregulated and 3995 downregulated genes as compared with the control group. The top 10 upregulated and downregulated mRNAs are shown in Table 1. The AECOPD group has 740 upregulated and 1373 downregulated genes as compared with the control group. The top 10 upregulated and downregulated mRNAs are shown in Table 2 (Figures 6A, B; Tables 2, 3). Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were used to investigate the differentially methylated genes (DMGs) and the biological significance of m6A methylation in the occurrence and development of COPD. GO analysis showed that during the occurrence and development of COPD, DMGs were mainly related to the regulation of Wnt signaling pathway, B cell receptor signaling pathway, natural killer cells regulating immunity, G-protein coupled receptor activity, and macrophage activation (Figure 6C). KEGG enrichment results showed the enrichment of different methylation genes in cancer pathway, PI3K signaling pathway, MAPK signaling pathway, and B cell receptor signaling pathway (Figure 6D). These results indicate that m6A methylation modifies many immune-related genes during the occurrence and development of COPD.

TABLE 2

Table 2 Top 10 differentially expressed mRNAs between stable COPD and control group based on log₂ (FC) value.

FIGURE 6

Figure 6 Gene Ontology and Kyoto Encyclopedia of Gene and Genomes analyses: (A, B) Heat map and volcano map show the m6A-modified differentially expressed genes (DEGs) between the control group and stable COPD and AECOPD groups. (C) GO terminologies of m6A-modified DEGs between the control group and stable COPD and AECOPD groups. (D) Twenty pathways of m6A-modified DEGs between the control group and stable COPD and AECOPD groups.

TABLE 3

Table 3 Top 10 differentially expressed mRNAs between AECOPD and control group based on log₂ (FC) value.

3.6 Differentially expressed genes (DEGs) during the occurrence and development of COPD

To determine the expression profiles of DEGs in COPD lung tissues, DEGs between the control group and stable COPD and AECOPD groups were compared. The results showed 983 upregulated genes and 363 downregulated genes in the stable COPD group as compared with the control group. In comparison with the control group, the AECOPD group exhibited 164 upregulated and 56 downregulated genes (Figure 7A). For DEGs, the lung tissues from the same group showed a similar expression trend (Figure 7B). Further, GO and KEGG enrichment analyses were conducted for these DEGs. The results of GO enrichment analysis showed that these genes were related to functions such as immune response and regulation, including inflammation, activation of the immune response, B cell receptor signaling pathway, natural killer cell regulating immunity, T cell regulating immunity, NF-κB signal positive regulation, and G protein-coupled receptor activity (Figure 7C). KEGG analysis results revealed their involvement with a series of inflammation-related signaling pathways, such as natural killer cells-induced cytotoxicity, MAPK signaling pathway, HIF-1 signaling pathway (Figure 7D).

FIGURE 7

Figure 7 DEGs during the occurrence and development of COPD: (A, B) Volcano map and thermal map show DEGs in lung tissues (|log2FC |>1, P < 0.05). (C) GO enrichment analysis of DEGs. (D) KEGG enrichment analysis of DEGs.

3.7 Joint analysis of m6A methylation and gene expression during the occurrence and development of COPD

MeRIP-seq and RNA-seq data were jointly analyzed and 986 and 445 differentially methylated mRNA transcripts were detected in the stable COPD and AECOPD group, respectively. Significant DEGs were found in these differentially methylated mRNA transcripts. In particular, 82 upregulated and 37 downregulated mRNAs were observed in the stable COPD group from 119 hypermethylated mRNAs and 419 upregulated and 448 downregulated mRNAs were obtained from 867 hypomethylated mRNAs. In the AECOPD group, 71 upregulated mRNAs and 16 downregulated mRNAs were obtained from 87 hypermethylated mRNAs and 115 upregulated and 243 downregulated mRNAs were detected from 358 hypomethylated mRNAs (Figure 8A). To clarify the relationship between methylation level and expression level of mRNAs, Pearson’s correlation analysis was carried out. The results showed a positive correlation between mRNA methylation level and mRNA expression level in both stable COPD and AECOPD groups (Figure 8B). This observation indicates that m6A methylation plays a key function in the regulation of gene expression in COPD. To determine the role of these differentially methylated mRNAs in the development of COPD, GO and KEGG enrichment analyses were conducted. These differentially methylated mRNAs were significantly enriched in GO terms such as protein-binding, inflammation, and immune system processes (Figure 8C). KEGG results showed that these differentially methylated mRNAs were significantly enriched in interleukin (IL)-17 signaling pathway, tumor necrosis factor (TNF) signaling pathway, natural killer cell-induced cytotoxicity, nuclear factor kappa B (NF-κB) signaling pathway, and other aspects, suggesting that m6A methylation may have an immunomodulatory effect in COPD (Figure 8D). Considering the pathogenesis of COPD involving immune disorders and excessive inflammation, it is necessary to explore whether there are mRNAs with m6A differential modification related to pathological process. We found many mRNAs related to immune regulation. Several representative mRNAs were detected, including IL-1β and CXCL12 in the stable COPD group and NFKBIA and CEBP-β in the AECOPD group. IGV was used to visualize the distribution of m6A modification in these mRNAs (Figure 9A). We examined the lung tissues of stable-COPD mice and AECOPD mice by RT-qPCR and showed that compared with the control group, the expression levels of IL-1β and CXCL12 mRNA were significantly increased in the stable-COPD group and CEBP-β mRNA in the AECOPD group, while the NFKBIA mRNA The differences in expression levels were not statistically significant (Figure 9B). The validation results were generally consistent with the predicted results.

FIGURE 8

Figure 8 Joint analysis of MeRIP-seq and RNA-seq data: (A) Four quadrant graph shows the differentially methylated peaks in differentially expressed mRNAs (| log2FC>1, P < 0.05). (B) There was a positive correlation between differentially methylated mRNAs and the mRNA expression level (0.3<Pearson R <0.5, P < 1e-10). (C, D) GO and KEGG enrichment analyses of mRNAs with differential methylation and expression.

FIGURE 9

Figure 9 Changes in IL-1β, CXCL12, NFKBIA, and CEBP-β m6a modifications and mRNA expression levels during the development of COPD: (A) IGV tracks display the regional distribution of peaks in IL-1β, CXCL12, NFKBIA, and CEBP-β genes. (B) Detection of IL-1β, CXCL12, NFKBIA and CEBP-β mRNA expression levels by RT-qPCR (*P < 0.05). ns, P=0.2078.

4 Discussion

According to the World Health Organization, COPD is the third leading cause of death globally and costs are estimated to be in excess of 100 billion dollars per year (15, 16). COPD can be divided into stable COPD and AECOPD depending on the aggravation of clinical symptoms and its pathogenesis mainly involves inflammation (17). m6A, a dynamic RNA modification, is a key regulator of gene expression and plays an important role in the pathogenesis of various inflammatory diseases, and may also be involved in the development of COPD (18–22). Herein, we found that the expression levels of m6A key regulatory factors, namely, METTL3, METTL14, FTO, and YTHDF1 mRNAs, significantly decreased in the COPD group. In summary, m6A is involved in the occurrence and development of COPD. However, studies comprehensively examining m6A modification characteristics in COPD animal models and exploring its effects on the immune response in COPD are warranted. To supplement the existing literature, we performed MeRIP-seq and RNA-seq to reveal the m6A methylation characteristics of mice with COPD at the whole transcriptome level for the first time, and analyzed the methylation genes related to COPD and their potential functions in detail.

4.1 m6A methylation model of COPD

There were obvious m6A modifications in the lung tissues of normal and COPD mice. We found that the common RRACH sequence of m6A modification was AAACC in COPD through Homer; future studies should investigate the underlying reason. In addition, the enrichment of m6A modification in 3′-UTR was much higher than that in 5′-UTR, which was consistent with the previously reported distribution trend of m6A modification in lung tissues of acute allergic asthmatic mice. m6A usually appeared near the termination codon and 3′-UTR (23). High levels of m6A modification at 3′-UTR may be related to mRNA stability, selective polyadenylation, signal transduction, and translocation (24, 25). Further, m6A modification at 3′-UTR plays a regulatory role in protein translation by recruiting specific factors to these m6A modification sites for RNA transport or protein synthesis, which may be one of the reasons for the potential positive correlation between the degree of m6A methylation and level of transcription (26, 27). In this study, approximately 80% of methylation transcripts included one or two m6A peaks and 20% of m6A transcripts included three or more m6A peaks. Meanwhile, the distribution characteristics of the m6A peak were roughly the same in the three groups, and were mainly concentrated in the CDS and 3′-UTR areas. This observation is consistent with the peak distribution trend for peripheral nerve injury, allergic asthma, and rectal cancer reported in mice (28–30). In summary, these results show that m6A modification plays a role in COPD.

4.2 The relationship between differentially methylated mRNAs and airway inflammation in COPD

We observed several unique m6A peaks and methylated mRNAs in lung tissues of COPD mice at each stage. To determine the biological functions of differentially methylated mRNAs in COPD, GO and KEGG enrichment analyses were carried out. GO enrichment revealed the biological functions of most differentially methylated mRNAs in both stable COPD and AECOPD groups that were closely related to immune functions such as regulation of Wnt signaling pathway, B cell receptor signaling pathway, natural killer cell regulating immunity, G-protein coupling receptor activity, and macrophage activation. Previous studies have shown that the Wnt signaling pathway exerts immunomodulatory effects and is closely related to inflammation in COPD patients. For instance, the Wnt/β-catenin pathway promotes airway inflammation and Th17/Treg imbalance in COPD (31). Wnt4 can increase IL-8, IL-6, and monocyte chemoattractant protein-1 levels in bronchial epithelial cells, which may lead to neutrophil infiltration and inflammation in COPD (32, 33). Li et al. found that upregulation of Wnt signaling can weaken Toll-like receptor signaling and mediate inflammation of alveolar epithelial cells (34). KEGG enrichment results showed that the differentially methylated mRNAs were mainly enriched in the B cell receptor signaling pathway, Toll-like receptor signaling pathway, and other immune regulatory pathways. Janus kinase (JAK) signaling pathway and signal transducer and activator of transcription(STAT) signaling pathway participates in the occurrence of multiple inflammatory diseases and contributes to the pathogenesis of COPD by generating cytokines such as IL-1β and IL-6 (35). It is speculated that these differentially methylated mRNAs play a role in the pathogenesis of COPD by regulating inflammation-related pathways. m6A metabolism affects the stability and degradation of mRNAs and consequently regulates their expression (6, 36).

4.3 m6A methylation regulates gene expression

To investigate whether m6A modification affects gene expression in COPD patients, we evaluated differentially methylated loci gene expression. First, the gene expression spectrum was analyzed, and 983 upregulated and 363 downregulated genes were observed in the stable COPD group and 164 upregulated and 56 downregulated genes were detected in the AECOPD group. These DEGs participate in the immune regulation mechanism of COPD by modulating pathways such as inflammation, immune response activation, B cell receptor signaling pathway, and natural killer cells regulating immunity.

We jointly analyzed DEGs and DMGs in COPD lung tissues, and found that the stable COPD group had 82 upregulated hypermethylated mRNAs, 37 downregulated hypermethylated mRNAs, 419 upregulated hypomethylated mRNAs, and 448 downregulated hypomethylated mRNAs. The AECOPD group also showed 71 upregulated hypermethylated mRNAs, 16 downregulated hypermethylated mRNAs, 115 upregulated hypomethylated mRNAs, and 243 downregulated hypomethylated mRNAs. Pearson’s correlation analysis revealed a positive correlation between the methylation level of these genes and their expression levels in the two COPD groups. Downregulated hypomethylated mRNAs accounted for the largest proportion of COPD lung tissues. Hence, m6A methylation plays a dominant role in maintaining mRNA stability in COPD. Enrichment analysis revealed that some of these mRNAs are related to immune functions, including inflammation, immune system processes, IL-17 signaling pathway, and TNF signaling pathway. We selected several representative differentially methylated mRNAs from the above pathways (IL-1β, CXLC12, NFKBIA, CEBP-β), reviewed their roles in immunity and validated them by RT-qPCR.

As a member of the cytokine superfamily, the C-X-C motif chemokine ligand 12 (CXCL12) plays a chemotactic role in immunocytes. It was first found that CXCL12 exerts an important function in the regulation of proliferation and differentiation of pre-B cells. CXCL12 is closely related to inflammatory diseases and acts as a co-stimulator for the activation and proliferation of CD4⁺T cells in patients with chronic lymphoblastic leukemia. CXCL12 is also necessary for the differentiation of B-lymphoid and plasmacytoid dendritic cells (37, 38). IL-1β is an important inflammatory cytokine involved in the development of inflammation, and participates in the occurrence of various inflammatory diseases such as atherosclerosis, type II diabetes, and various autoimmune diseases (39). IL-1β also plays an important role in COPD by mediating the release of inflammatory substances such as colony-stimulating factors, cell adhesion factors, and hypersensitive response proteins (40). NFKBIA is a member of the NF-κB family that is closely related to head and neck, nasopharyngeal, and other cancers. It is a therapeutic target of small cell lung cancer (41–43). NFKBIA mutations promote the formation of IL-1β, and give rise to severe immune deficiency in liver diseases (44). A previous meta-analysis found that NFKBIA was closely related to autoimmune diseases and susceptibility to inflammatory diseases (45). As a transcription factor, CCAAT/enhancer-binding protein beta (CEBP-β) participates in the expression of inflammatory factors and the occurrence of various inflammatory diseases. It is one of the core genes of inflammation and infection prevention that can participate in the occurrence of COPD by regulating TNF-α and IL-17. Further, CEBP-β also participates in the mesenchymal transformation of bronchial epithelium by regulating IL-17 expression and thereby affecting the development of COPD (46, 47). Recent studies have shown that m6A reading protein IMP2 participates in the occurrence of autoimmune kidney injury by enhancing m6A modification of CEBP-β, promoting its stability, and upregulating the CEBP-β–mediated transcription of IL-17 (47). We then performed RT-qPCR on lung tissues from stable-COPD and AECOPD mice, and the results showed that compared with the control group, the expression levels of IL-1β and CXCL12 mRNA were significantly increased in the stable-COPD group and CEBP-β mRNA in the AECOPD group, while the NFKBIA mRNA The differences in expression levels were not statistically significant. The validation results were generally consistent with the predicted results. Refer to predictions, validation results and literature above, we speculate that m6A regulates the stability of RNAs such as CXCL12, IL-1β, and CEBP-β, increases their expression, and contribute to COPD.

5 Conclusions

A series of differentially methylated mRNAs that may be important regulators in COPD pathogenesis were determined. The first mouse m6A modification map in COPD lung tissues was drawn. This study provides a new direction for understanding the mechanism of m6A modification affecting inflammation and screening potential therapeutic targets for COPD.

6 Limitations

This study had potential limitations, although the target genes modified by m6A were known, the process of methylation readers, erasers or writers regulating the target genes was not delineated. Further studies should be conducted to investigate whether readers, erasers or writers regulate the stability, translation efficiency, or degradation of target genes. We found the above results in mice, but it is unclear whether the results are consistent in human lung tissue. In the future, we plan to collect samples from COPD patients and perform m6a-related assays on them.

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: PRJNA853736 and PRJNA919075 (NCBI Bioproject).

Ethics statement

The animal study was reviewed and approved by Experimental Animal Ethics Committee of Xinjiang Medical University (approval no. SYXK-2018-0003).

Author contributions

JW conceived and designed the experiments. TH, LX, MJ, FZ, QL, ZL, CW, JD, and FL performed the experiments. TH and LX analyzed the data. JW, TH, LX collected the fund. TH wrote the manuscript. MJ and JW revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was jointly supported by the Outstanding Young Scientists Fund of Xinjiang Uygur Autonomous Region Natural Science Foundation project (Project No: 2022D01E28), the Young Scientists Fund of Xinjiang Uygur Autonomous Region Natural Science Foundation project (Project No: 2022D01C176 and 2022D01C179), the Xinjiang Key Laboratory of Respiratory Disease Research Open Subjects (Project No: ZYYHX202104), the National Natural Science Foundation Regional Fund Project (Project No: 82160844 and 81960848), and the Xinjiang Uyghur Autonomous Region Key Laboratory Open Subjects[Project No: 2021D0423].

Acknowledgments

Thanks to Zheng Li, Dan Xu and Jing Jing for their help with experiments and specimen collection.

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1137195/full#supplementary-material

References

1. López-Campos JL, Tan W, Soriano JB. Global burden of COPD. Respirology (2016) 21(1):14–23. doi: 10.1111/resp.12660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Rabe KF, Watz H. Chronic obstructive pulmonary disease. Lancet (2017) 389(10082):1931–40. doi: 10.1016/S0140-6736(17)31222-9

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bontsevich RA, Adonina AV, Vovk YR, Batisheva GA, Cherenkova OV, Ketova GG, et al. Management of chronic obstructive pulmonary disease. Arch Razi Inst (2022) 77(1):439–47.

PubMed Abstract | Google Scholar

4. Li Y, Qian H, Yu K, Huang Y. The long-term maintenance effect of remote pulmonary rehabilitation via social media in COPD: A randomized controlled trial. Int J Chron Obstruct Pulmon Dis (2022) 17:1131–42. doi: 10.2147/COPD.S360125

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Liu ZX, Li LM, Sun HL, Liu SM. Link between m6A modification and cancers. Front Bioeng Biotechnol (2018) 6:89. doi: 10.3389/fbioe.2018.00089

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Csepany T, Lin A, Baldick CJ Jr, Beemon K. Sequence specificity of mRNA N6-adenosine methyltransferase. J Biol Chem (1990) 265(33):20117–22. doi: 10.1016/S0021-9258(17)30477-5

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature (2012) 485(7397):201–6. doi: 10.1038/nature11112

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3' UTRs and near stop codons. Cell (2012) 149(7):1635–46. doi: 10.1016/j.cell.2012.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Yu J, Hu X, Chen H, Yang Q, Li H, Dong Y, et al. DNA Methylation of FTO promotes renal inflammation by enhancing mA of PPAR-α in alcohol-induced kidney injury. (2021) 163:105286. doi: 10.1016/j.phrs.2020.105286

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Guo M, Yan R, Ji Q, Yao H, Sun M, Duan L, et al. IFN regulatory factor-1 induced macrophage pyroptosis by modulating m6A modification of circ_0029589 in patients with acute coronary syndrome. (2020) 86:106800. doi: 10.1016/j.intimp.2020.106800

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Lu T, Zheng Z, Zhang L, Sun H, Bissonnette M, Huang H, et al. A new model of spontaneous colitis in mice induced by deletion of an RNA mA methyltransferase component METTL14 in T cells. (2020) 10(4):747–61. doi: 10.1016/j.jcmgh.2020.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Lu W, Li D, Hu J, Mei H, Shu J, Long Z, et al. Hydrogen gas inhalation protects against cigarette smoke-induced COPD development in mice. J Thorac Dis (2018) 10(6):3232–43. doi: 10.21037/jtd.2018.05.93

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Li D, Sun D, Yuan L, Liu C, Chen L, Xu G, et al. Sodium tanshinone IIA sulfonate protects against acute exacerbation of cigarette smoke-induced chronic obstructive pulmonary disease in mice. Int Immunopharmacol (2020) 81:106261. doi: 10.1016/j.intimp.2020.106261

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Guarascio AJ, Ray SM, Finch CK, Self TH. The clinical and economic burden of chronic obstructive pulmonary disease in the USA. Clinicoecon Outcomes Res (2013) 5:235–45. doi: 10.2147/ceor.S34321

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Sandelowsky H, Weinreich UM, Aarli BB, Sundh J, Høines K, Stratelis G, et al. COPD - do the right thing. BMC Fam Pract (2021) 22(1):244. doi: 10.1186/s12875-021-01583-w

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Sun X, Liu Y, Feng X, Li C, Li S, Zhao Z. The key role of macrophage depolarization in the treatment of COPD with ergosterol both in vitro and in vivo. Int Immunopharmacol (2020) 79:106086. doi: 10.1016/j.intimp.2019.106086

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Shen J, Wang W, Shao X, Wu J, Li S, Che X, et al. Integrated analysis of m6A methylome in cisplatin-induced acute kidney injury and berberine alleviation in mouse. Front Genet (2020) 11:584460. doi: 10.3389/fgene.2020.584460

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Teng F, Tang W, Wuniqiemu T, Qin J, Zhou Y, Huang X, et al. N (6)-methyladenosine methylomic landscape of lung tissues in murine acute allergic asthma. Front Immunol (2021) 12:740571. doi: 10.3389/fimmu.2021.740571

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Jiang H, Cao K, Fan C, Cui X, Ma Y, Liu J. Transcriptome-wide high-throughput m6A sequencing of differential m6A methylation patterns in the human rheumatoid arthritis fibroblast-like synoviocytes cell line MH7A. J Inflammation Res (2021) 14:575–86. doi: 10.2147/JIR.S296006

CrossRef Full Text | Google Scholar

21. Guo X, Lin Y, Lin Y, Zhong Y, Yu H, Huang Y, et al. PM2.5 induces pulmonary microvascular injury in COPD via METTL16-mediated m6A modification. Environ pollut (2022) 303:119115. doi: 10.1016/j.envpol.2022.119115

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Huang X, Lv D, Yang X, Li M, Zhang H. m6A RNA methylation regulators could contribute to the occurrence of chronic obstructive pulmonary disease. J Cell Mol Med (2020) 24(21):12706–15. doi: 10.1111/jcmm.15848

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Teng F, Tang W, Wuniqiemu T, Qin J, Zhou Y, Huang X, et al. N -methyladenosine methylomic landscape of lung tissues in murine acute allergic asthma. (2021) 12:740571. doi: 10.3389/fimmu.2021.740571

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Yue Y, Liu J, Cui X, Cao J, Luo G, Zhang Z, et al. VIRMA mediates preferential m(6)A mRNA methylation in 3'UTR and near stop codon and associates with alternative polyadenylation. Cell Discovery (2018) 4:10. doi: 10.1038/s41421-018-0019-0

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Shen L, Liang Z, Gu X, Chen Y, Teo ZW, Hou X, et al. N(6)-methyladenosine RNA modification regulates shoot stem cell fate in arabidopsis. Dev Cell (2016) 38(2):186–200. doi: 10.1016/j.devcel.2016.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Niu Y, Zhao X, Wu YS, Li MM, Wang XJ, Yang YG. N6-methyl-adenosine (m6A) in RNA: an old modification with a novel epigenetic function. Genomics Proteomics Bioinf (2013) 11(1):8–17. doi: 10.1016/j.gpb.2012.12.002

CrossRef Full Text | Google Scholar

27. Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, et al. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature (2014) 505(7481):117–20. doi: 10.1038/nature12730

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Zhang L, Hao D, Ma P, Ma B, Qin J, Tian G, et al. Epitranscriptomic analysis of m6A methylome after peripheral nerve injury. Front Genet (2021) 12:686000. doi: 10.3389/fgene.2021.686000

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Sun D, Cai X, Shen F, Fan L, Yang H, Zheng S, et al. Transcriptome-wide m6A methylome and m6A-modified gene analysis in asthma. Front Cell Dev Biol (2022) 10:799459. doi: 10.3389/fcell.2022.799459

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhang Z, Wang Q, Zhang M, Zhang W, Zhao L, Yang C, et al. Comprehensive analysis of the transcriptome-wide m6A methylome in colorectal cancer by MeRIP sequencing. Epigenetics (2021) 16(4):425–35. doi: 10.1080/15592294.2020.1805684

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhou M, Jiao L, Liu Y. sFRP2 promotes airway inflammation and Th17/Treg imbalance in COPD via wnt/β-catenin pathway. Respir Physiol Neurobiol (2019) 270:103282. doi: 10.1016/j.resp.2019.103282

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Heijink IH, de Bruin HG, van den Berge M, Bennink LJ, Brandenburg SM, Gosens R, et al. Role of aberrant WNT signalling in the airway epithelial response to cigarette smoke in chronic obstructive pulmonary disease. Thorax (2013) 68(8):709–16. doi: 10.1136/thoraxjnl-2012-201667

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Durham AL, McLaren A, Hayes BP, Caramori G, Clayton CL, Barnes PJ, et al. Regulation of Wnt4 in chronic obstructive pulmonary disease. FASEB J (2013) 27(6):2367–81. doi: 10.1096/fj.12-217083

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Li Y, Shi J, Yang J, Ma Y, Cheng L, Zeng J, et al. A wnt/β-catenin negative feedback loop represses TLR-triggered inflammatory responses in alveolar epithelial cells. Mol Immunol (2014) 59(2):128–35. doi: 10.1016/j.molimm.2014.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Hansbro PM, Haw TJ, Starkey MR, Miyake K. Toll-like receptors in COPD. Eur Respir J (2017) 49(5):1–3. doi: 10.1183/13993003.00739-2017

CrossRef Full Text | Google Scholar

36. Meyer KD, Jaffrey SR. The dynamic epitranscriptome: N6-methyladenosine and gene expression control. Nat Rev Mol Cell Biol (2014) 15(5):313–26. doi: 10.1038/nrm3785

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Borge M, Nannini PR, Morande PE, Jancic C, Bistmans A, Bezares RF, et al. CXCL12 is a costimulator for CD4+ T cell activation and proliferation in chronic lymphocytic leukemia patients. Cancer Immunol Immunother (2013) 62(1):113–24. doi: 10.1007/s00262-012-1320-7

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Minami H, Nagaharu K, Nakamori Y, Ohishi K, Shimojo N, Kageyama Y, et al. CXCL12-CXCR4 axis is required for contact-mediated human b lymphoid and plasmacytoid dendritic cell differentiation but not T lymphoid generation. J Immunol (2017) 199(7):2343–55. doi: 10.4049/jimmunol.1700054

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Yazdi AS, Ghoreschi K. The interleukin-1 family. Adv Exp Med Biol (2016) 941:21–9. doi: 10.1007/978-94-024-0921-5_2

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zou Y, Chen X, Liu J, Zhou DB, Kuang X, Xiao J, et al. Serum IL-1β and IL-17 levels in patients with COPD: associations with clinical parameters. Int J Chron Obstruct Pulmon Dis (2017) 12:1247–54. doi: 10.2147/COPD.S131877

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Gong L, Zhang D, Dong Y, Lei Y, Qian Y, Tan X, et al. Integrated bioinformatics analysis for identificating the therapeutic targets of aspirin in small cell lung cancer. J BioMed Inform (2018) 88:20–8. doi: 10.1016/j.jbi.2018.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Li L, Zhang ZT. Genetic association between NFKBIA and NFKB1 gene polymorphisms and the susceptibility to head and neck cancer: A meta-analysis. Dis Markers (2019) 2019:6523837. doi: 10.1155/2019/6523837

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Zhang H, Deng S, Zhang J, Zhu G, Zhou J, Ye W, et al. Single nucleotide polymorphisms within NFKBIA are associated with nasopharyngeal carcinoma susceptibility in Chinese han population. Cytokine (2021) 138:155356. doi: 10.1016/j.cyto.2020.155356

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Tan EE, Hopkins RA, Lim CK, Jamuar SS, Ong C, Thoon KC, et al. Dominant-negative NFKBIA mutation promotes IL-1β production causing hepatic disease with severe immunodeficiency. J Clin Invest (2020) 130(11):5817–32. doi: 10.1172/JCI98882

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Zhang GL, Zou YF, Feng XL, Shi HJ, Du XF, Shao MH, et al. Association of the NFKBIA gene polymorphisms with susceptibility to autoimmune and inflammatory diseases: a meta-analysis. Inflammation Res (2011) 60(1):11–8. doi: 10.1007/s00011-010-0216-2

CrossRef Full Text | Google Scholar

46. Chu S, Ma L, Wu Y, Zhao X, Xiao B, Pan Q. C-EBPβ mediates in cigarette/IL-17A-induced bronchial epithelial-mesenchymal transition in COPD mice. BMC Pulm Med (2021) 21(1):376. doi: 10.1186/s12890-021-01738-6

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Bechara R, Amatya N, Bailey RD, Li Y, Aggor FEY, Li DD, et al. The m(6)A reader IMP2 directs autoimmune inflammation through an IL-17- and TNFα-dependent C/EBP transcription factor axis. Sci Immunol (2021) 6(61):1–34. doi: 10.1126/sciimmunol.abd1287

CrossRef Full Text | Google Scholar