Comprehensive bioinformatic analysis identi ed m6A methylation as an important regulator in mediating immune response during idiopathic pulmonary arterial hypertension

Na Liu Department of Cardiology/Cardiac Catheterization Lab, Second Xiangya Hospital, Central South University Yunhong Zeng Department of Cardiology, Hunan Children's Hospital Ting Huang Department of Cardiology, Hunan Children's Hospital Wanyun Zuo Department of Cardiovascular Medicine, The Second Xiangya Hospital of Central South University, Changsha, Hunan, China Yunbin Xiao Department of Cardiology, Hunan Children's Hospital, Changsha, Hunan, China Biao Li Department of Cardiovascular Medicine, The Second Xiangya Hospital of Central South University, Changsha, Hunan, China Zhenghui Xiao Department of Cardiology, Hunan Children's Hospital, Changsha, Hunan, China Yaozhong Liu Department of Cardiovascular Medicine, The Second Xiangya Hospital of Central South University, Changsha, Hunan, China Qiming Liu (  qimingliu@csu.edu.cn ) Second Xiangya Hospital


Background
Pulmonary arterial hypertension (PAH) is a condition characterized by increased pulmonary arterial pressure and vascular remodeling of the small pulmonary arteries. PAH develops either as an idiopathic condition or in association with various underlying diseases such as collagen vascular disease, portal hypertension, or HIV infection (1,2). The incidence of PAH ranges from 2.0 to 7.6 cases per million adults per year, and its prevalence varies from 11 to 26 cases per million adults. The incidence of PAH is fourfold higher in women than in men, but survival is paradoxically worse in men with PAH (3,4). Despite striking progress in the development of diagnosis and therapeutics in the last three decades, PAH remains relatively incurable and the limited e cacy of current treatment options possibly results from an incomplete understanding of the pathophysiology of this cardiopulmonary disorder (5).
Epigenetic modi cations,such as DNA methylation, histone modi cations, and noncoding RNAs, may play important roles in regulating the pathogenesis of PAH (1). As a bridge to pass the genetic information from DNA to proteins RNA, RNA is an important part of the central dogma, and its various chemical modi cations are involved in the regulation of various biological processes (6). N6-methyladenosine (m 6 A) of RNA transcripts is the most prevalent modi cation found in many classes of RNA (7,8). The abundance of m 6 A on RNA is determined by the dynamic interplay between its methyltransferases ("writers", such as METTL3-METTL14 complex(9)), and demethylases ("erasers", including FTO(10) and ALKBH5 (11)). The binding proteins ("readers", such as YTH domain-containing proteins (12,13)) then mediate the effect of m 6 A modi cations on RNA processing or metabolisms, including alternative splicing (14,15), export (16,17), stability (18,19), and translation (20,21).
Recent studies have suggested the involvement of the immune response during PAH. In ammatory cell in ltrations can be observed around the pulmonary artery, such as macrophages dendritic cells, mast cells, T lymphocytes, and B lymphocytes (22). These immune cells play an important role in the vascular remodeling characteristic of PAH and might be important targets for PAH therapy. Nevertheless, how the immune system is activated and regulated remains poorly understood. m 6 A RNA modi cations might be a novel immunoregulatory factor (23).
Weighted gene co-expression network analysis (WGCNA) is a system biology algorithm for describing the correlation patterns among genes across samples (24). WGCNA can identify and cluster highly correlated genes into the same module. Furthermore, this method can also be used for relating modules to external clinical traits. By far, WGCNA has been validated as a valuable method to identify underlying mechanisms, potential biomarkers, or therapeutic targets in different types of diseases.
In our study, we rst constructed a PAH rat model and conducted the m 6 A methylation sequencing analysis. We then analyzed 18 idiopathic PAH (IPAH) patients microarray data from the GEO database to identify m 6 A methylation correlated modules by using the method of WGCNA. Afterward, we applied the CIBERSORT algorithm to investigate the effect of m 6 A methylation on immune cell in ltration during PAH. Findings from our study may contribute to novel regulators and therapeutic targets in PAH from an epigenomic perspective.

Establishment and veri cation of rat model of PAH induced by MCT
The study was approved by the Institutional Animal Care and Use Committee of the Second Xiangya Hospital of Central South University and complied with the standards in the Guide for the Care and Use of Laboratory Animals. Sprague-Dawley rats (speci c pathogen-free (SPF), male, 180-200 g, 6 weeks old, n = 5) were obtained from Changsha Tianqin Biotechnology Company (China). The rats were randomized to the control (n = 2) and PAH groups (n = 3). The rats in the PAH group were intraperitoneally injected with monocrotaline (MCT) (60 mg·kg − 1 , Sigma, C2401), while the rats in the control group were injected intraperitoneally with the same volume of saline. All rats were housed on a 12 h light/dark cycle with free access to food and water. After 4 weeks of feeding, rats were anesthetized with 1% sodium pentobarbital (130 mg·kg − 1 ) for echocardiography and right heart catheterization. Echocardiography was used to record the tricuspid regurgitation velocity, tricuspid annulus contraction rate, tricuspid annular plane systolic excursion (TAPSE) and pulmonary artery blood ow acceleration time (PAAT). After echocardiographic examination, We performed a right heart catheter to measure pulmonary artery pressure.
Rats were sacri ced by cervical dislocation after deep anesthesia. Then, heart tissues were removed and segregated. The ratio of right ventricle to left ventricle plus ventricular septum [RV/ (LV + S)] was used as an index of right ventricular hypertrophy. Lung tissues were excised and immediately frozen at liquid nitrogen or xed in 4% buffered paraformaldehyde solution.

Histological analysis
The lung tissues obtained in each group were placed in 4% buffered paraformaldehyde solution overnight, and then dehydrated and embedded in para n. Then, all lung tissues were sliced into 5 µmthick sections, xed on a glass slide and baked dry. The staining procedures were performed according to the instructions. Brie y, the sections were soaked in xylene, ethanol in gradient concentration and hematoxylin, respectively, and sealed with resin. After dryness, pulmonary vascular morphology was observed and photographed under optical microscope. Last, the ratio of pulmonary small artery wall thickness and muscularization was calculated.

RNA preparation
For each group, four biological replicates were selected, of which every two were combined into one. Then, total RNA of tissue was extracted using TRIzol reagent (Invitrogen Corporation, CA, USA) following the manufacturer's instructions. The Ribo-Zero rRNA Removal Kit (Illumina, Inc., CA, USA) was used to reduce the ribosomal RNA content of total RNAs. Then, the RNA was chemically fragmented into fragments of about 100 nucleotides in length using fragmentation buffer (Illumina, Inc.).

Methylated RNA immunoprecipitation sequencing (MeRIP-Seq) and library construction
MeRIP-Seq was performed following a previously reported procedure with slight modi cations. Brie y, m 6 A RNA immunoprecipitation was performed with the GenSeqTM m 6 A RNA IP Kit (GenSeq Inc., China) by following the manufacturer's instructions. Both the input sample without immunoprecipitation and the m 6 A IP samples were used for RNA-seq library generation with NEBNext® Ultra II Directional RNA Library Prep Kit (New England Biolabs, Inc., USA). The library quality was evaluated with BioAnalyzer 2100 system (Agilent Technologies, Inc., USA). Library sequencing was performed on an illumina Hiseq instrument with 150 bp paired-end reads.

Microarray data collection and processing
Materials of GSE15197 were downloaded from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) database. The dataset contains lung tissue transcriptomic pro les from 18 IPAH patients. The median expression values among all multiple probe IDs were selected to represent the corresponding gene symbol, leading to the identi cation of 18612 unique genes across 18 samples. Human gene annotation le (GRCh38 release 99 gene transfer format, ensembl.org/index.html) was applied to annotate genes and the 15934 protein-coding genes were selected for further analysis.

Estimation of m 6 A methylation score and implementation of weight correlation network analysis (WGCNA)
We constructed gene signatures of m 6 A writer (METTL3, METTL14, METTL16, WTAP, KIAA1429, ZC3H13, RBM15/RBM15B), and m 6 A eraser (ALKBH5, FTO) as suggested by Yang et.al (25). We then computed the Gene Set Variation Analysis (GSVA) enrichment score of m 6 A writer and m 6 A eraser across the 18 samples using the 'GSVA' package (26) in R software. The estimated m 6 A methylation score was calculated by substrating m 6 A eraser score from m 6 A writer score. WGCNA was accomplished with the R package 'WGCNA' (24). According to the standard variation value of gene, we ranked them from largest to smallest and only selected the top 2000 as input for WGCNA. A power of β value was introduced so that it could transform the similarity matrix into an adjacency matrix. In this study, β = 5 was used as a soft threshold parameter to ensure a scale-free network. The Dynamic Tree Cut method was applied to generate modules with the following major parameters to avoid the generation of too many modules: deepSplit of 2, minModuleSize of 30, and the height cut-off was set as 0.25 (modules were merged if their similarity was > 0.75). Module eigengenes (MEs) referred to the rst principal component of all gene expression levels in the module, and therefore, it was reasonable to consider that MEs represented all genes within a speci c module. According to Pearson's correlation tests, we further identi ed the association between MEs and m 6 A methylation score. Within the most relavant module, those positively correlative with m 6 A methylation score were made subjected to further analysis.

Pathway enrichment analysis
Metascape (https://Metascape.org/) is a web-based portal designed to provide a comprehensive gene list annotation and analysis resource for biologists (27). To gain insights into biological roles of m 6 A methylation correlated genes identi ed from WGCNA and upmethylated genes identi ed from m 6 A-seq, we conducted pathway enrichment analysis in Metascape tools using Gene Ontology biological process (GO BP), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome (27). By inputting the lists of m 6 A methylation correlated genes and upmethylated genes simultaneously, Metascape can identify commonly-enriched and selectively-enriched pathways from two levels, which enable a comprehensive assessment of the molecular features of the biological process. For Multi-list Enrichment Analysis, Metascape rst applies enrichment analysis to each gene list individually and identi es terms that are statistically enriched. All gene lists are then combined into one new list, and enrichment analysis is conducted on this combined list. Distinguishing it from many existing portals, Metascape automatically clusters enriched terms into non-redundant groups. Brie y, pairwise similarities between any two enriched terms are computed based on a Kappa-test score. The similarity matrix is then hierarchically clustered and a 0.3 similarity threshold is applied to trim the resultant tree into separate clusters. We selected the top 20 clusters and chose the most signi cant (lowest P-value) term within each of them to represent it.

Protein-protein interaction (PPI) network construction of m 6 A methylation correlated genes
To nd out the functional associations among the m 6 A methylation correlated genes, we used the online Search Tool for the Retrieval of Interacting Genes (STRING database; http://string-db.org/) to construct a PPI network based on uniquely comprehensive coverage and predictive function of genome-wide data (28). A stringent threshold of a combined score of > 0.7 was used to construct the PPI network and Cytoscape software was used to visualize and analyze the biological networks. Plugin Molecular Complex Detection (MCODE) was applied to identify signi cant clusters with strong protein-protein linkages with default parameters.

Estimation of immunocyte in ltration
The CIBERSORT algorithm (https://cibersort.stanford.edu/) was applied to estimate the abundance of in ltrated immune cell subtypes in the lung tissues of 18 samples, based on a deconvolution algorithm in the R software (29). The perm parameter was set as 1000. The 18 samples were grouped to 'High m 6 A methylation' and 'Low m 6 A methylation' by the median value of their m 6 A methylation score.
Comparisons between two groups were tested by the Wilcox rank test. P-value < 0.05 was considered signi cant.

Statistical analysis
Brie y, Paired-end reads were harvested from Illumina HiSeq 4000 sequencer, and were quality controlled by Q30. After 3' adaptor-trimming and low quality reads removing by cutadapt software (v1.9.3). First, clean reads of all libraries were aligned to the reference genome (UCSC RN5) by Hisat2 software (v2.0.4). Methylated sites on RNAs (peaks) were identi ed by MACS software. Differentially methylated sites with a fold change cutoff of ≥ 2 and false discovery rate cutoff of ≤ 0.0001 were identi ed with the diffReps differential analysis package. These peaks identi ed by both softwares overlapping with exons of mRNA were gured out and choosed by home-made scripts.

Construction of rats PAH model
The pulmonary artery velocity diagram appeared as a dagger in PAH group, and PAAT consumedly shortened compared to control group (20.33 ± 2.62 vs 29.50 ± 0.50 ms, p < 0.0001). Then we observed ventricular septum signi cantly shift to the left, simultaneously, right atrium and right ventricle enlarged compared to control group (4.97 ± 0.30 vs 3.42 ± 0.09 mm, p < 0.0001). The range of TAPSE was greatly reduced in PAH (1.33 ± 0.05 vs 1.58 ± 0.04 mm, p < 0.0001) (Fig. 1.A). Four weeks after the MCT injection, the right ventricular systolic pressure (RVSP) of the PAH group was elevated to 49.87 ± 1.17 mmHg compared with 24.84 ± 0.34 mmHg in the control (p < 0.001), and at the same time, we found that RV/ (LV + S) increased in PAH group (0.51 ± 0.03 vs 0.25 ± 0.04 g, p < 0.0001) (Fig. 1.B). The results of HE staining showed that the pulmonary artery smooth muscle was signi cantly thickened in PAH group compared with control group (71.13 ± 4.02% vs. 17.95 ± 1.35%, p < 0.001) (Fig. 1.C). These results indicate that the rat model of PAH has been successfully established.

Methylations pro le of rat lung tissue
MeRIP-Seq analysis identi ed 922 non-overlapping m 6 A sites in control group and 9059 non-overlapping m 6 A sites in PAH group, while 18655 m 6 A sites were overlapped in two groups (Fig. 2.A). Next, differentially methylated m 6 A sites (DMMSs) were identi ed by diffReps with the following default screening criteria: p-value ≤ 0.0001 and fold change ≥ 2 between the groups. We selected 3298 DMMSs in the two groups. A total of 777 m 6 A sites exhibited up-methylation, and 2521 exhibited downmethylation. (Fig. 2.B). These results demonstrated an increased m 6 A modi cation process during PAH and the 2521 up-methylated m 6 A sites representing 1261 unique genes were selected for further analysis.

Construction of weighted gene co-expression network and identi cation of key modules
The top 2000 genes in 18 samples of PAH lung tissue were used to construct the co-expression network. β values of 5 were con rmed to obtain the approximate scale-free topology with a scale-free topology t index > 0.85 and the lowest power (Fig. 3.A). Next, the method of dynamic tree cutting was employed to produce co-expression modules, which led to the identi cation of 11 modules (Fig. 3.B). We then calculated and plotted the relation of each module with the m 6 A methylation score. Among these modules, the yellow module depicting the highest correlation (module-trait weighted correlation = 0.46, p = 0.055) with m 6 A methylation score (Figrue 4). Within the yellow module, a total of 87 genes were positively correlative with m 6 A methylation score and were made subjected to further analysis.

Multi-list pathway enrichment analysis and visualization
We conducted a multi-list pathway enrichment analysis in Metascape by inputting the list of 87 m 6 A methylation correlated genes identi ed from WGCNA and 1261 m 6 A upmethylated genes identi ed from MeRIP-SEq. As Fig. 5 shown, both m 6 A correlated genes and upmethylated genes are involved in immune response-related pathways such as 'cytokine-mediated signaling pathway', 'leukocyte migration', and 'activation of immune response'. Besides, the upmethylated genes also participate in some biological processes already known to be associate with PAH such as 'ECM organization' and 'supermolecular ber organization'. In sum, the results demonstrated that m 6 A methylation may play an important role in PAH pathogenesis, especially in mediating immune response during PAH.

Identi cation of hub genes in the yellow module through PPI analysis
To identi y hub genes, we submitted the gene list of 87 m 6 A correlated genes in the yellow module to STRING for protein-protein interaction analysis, and the cut-off con dence interval was set to 0.7. After MCODE algorithm, the most signi cant cluster containing 10 hub genes was identi ed and enrichment analysis shows that genes in this cluster were mainly involved in 'Chemokine receptors bind chemokines' pathway (Fig. 6). By comparing the results of MeRIP-Seq, CCR5 and CXCL9 were found to be signi cantly upmethylated with a fold change of 2.55 and 6.93, respectively.

Correlation analysis between m 6 A methylation and immunocyte proportion
We have concluded that m 6 A methylation can mediate immune response during PAH process. We then aimed to investigate how it regulates immune cell in ltration. After the Wilcox rank test, monocyte and M1 macrophage were detected differentially in ltrated between the high methylation group and low methylation group (p = 0.03, and p = 0.001, respectively) (Fig. 7). Figure 8 displayed the scale histogram of immune cell fraction among 18 samples. These results indicated that during PAH, the m 6 A methylation can drive monocyte to form M1 macrophage, which plays vital roles in the immune response.
[Insert Fig. 7 and Fig. 8] 4. Discussion m 6 A is the most prevalent modi cation in the mRNA of many eukaryotic species. There is growing evidence that m 6 A dysregulation has a profound impact on the pathogenesis of various human diseases (30)(31)(32)(33). Su et al. rstly identi ed the transcriptome-wide map of m 6 A circRNAs in hypoxic PAH rat model and con rmed that m 6 A methylation can affect the circRNA-miRNA-mRNA network (34).
However, the mechanism of mRNA m 6 A methylation modi cation to regulate IPAH has not been reported. In this study, we took advantage of MeRIP-seq to map the transcriptomic landscape of PAH rat lung tissue and thus quantitatively compared transcriptome-wide changes between PAH and control groups.
The results found that the m 6 A level was signi cantly up-regulated in MCT induced PAH lung tissue, which highlighted the role of m 6 A methylation in the pathogenesis of PAH.
WGCNA is a method for constructing gene co-expression networks based on gene expression data. We applied this method to identi ed key modules associated with the estimated m 6 A methylation score.
Through multi-list enrichment analysis, the 87 m 6 A methylation correlated genes identi ed from WGCNA and 1261 m 6 A up-methylated genes identi ed from MeRIP-seq analysis are both signi cantly involved in immune response-related pathways such as 'cytokine-mediated signaling pathway', 'leukocyte migration', and 'activation of immune response'. These results highly suggested the immunoregulatory role of m 6 A methylation during PAH.
Through PPI analysis of the 87 m 6 A methylation correlated genes, a cluster containing 10 hub genes were identi ed, including CXCL9, CXCL10, CXCR3, PMCH, CXCL11, CXCR6, CCR7, CCL5, CCR5, and HCAR3. Among them, CCR5 and CXCL9 were signi cantly up-methylated. CCR5 is activated on stimulation by the CCR5 ligands CCL3 (macrophage in ammatory protein-1α), CCL4 (macrophage in ammatory protein-1β), and CCL5 (RANTES) and is strongly expressed on the principal cell types implicated in PAH progression, including endothelial cells, SMCs, T cells, and macrophages (35)(36)(37)(38). Amsellem et al. studied the effect of CCR5 receptor antagonists on PASMC and in ammatory response in PAH mouse model and found that the activation of the CCL5-CCR5 axis directly leads to PASMC proliferation and macrophage recruitment (39). CXCL9, a chemokine, is a T-cell chemoattractant that is induced by IFN-γ. CXCL9 is closely related to two other CXC chemokines called CXCL10 and CXCL11, whose genes are located near the gene for CXCL9 on human chromosome 4 (40,41). CXCL9, CXCL10, and CXCL11 are commonly produced by local cells in in ammatory lesions and can attract Th1 cells (42).
A high level of CXCL9 in peripheral liquids can be considered as a marker of host immune response, especially of that involving Th1 cells (43,44).
By applying the CYBERSORT algorithm, we found that the m 6 A methylation can drive monocyte to form M1 macrophage. Circulating monocytes migrate into the majority of tissues in the body, where they differentiate into functionally distinct mature macrophages (45,46). Macrophages undergo classical M1 activation induced by IFN-γ, which mediates Th1 cell-type activation of macrophages (47). The activated macrophage then causes vasoconstriction, increases vascular permeability, and induce PASMCs proliferation (48). Our research sheld light on the regulatory role of m 6 A methylation in macrophage activation during PAH. As we discussed above, The hub 10 m 6 A methylation correlated genes,especially CCR5 and CXCL9 may mediate m 6 A methylation caused M1 macrophage activation.

Conclusions
In conclusion, our research revealed that m 6 A methylation modi cation may play important roles in mediating immune response during PAH. It also caused activation of M1macrophage, which may be mediated by CCR5 and CXCL9. These results will help us to better understand the mechanisms of PAH, and provide candidate therapeutic targets.

Availability of data and materials
All data in our study are available upon request.

Competing interests
The authors declare that they have no competing interests.