2.1.1 Establishment and verification 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 (specific 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 flow acceleration time (PAAT). After echocardiographic examination, We performed a right heart catheter to measure pulmonary artery pressure.
Rats were sacrificed 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 fixed in 4% buffered paraformaldehyde solution.
2.1.2 Histological analysis
The lung tissues obtained in each group were placed in 4% buffered paraformaldehyde solution overnight, and then dehydrated and embedded in paraffin. Then, all lung tissues were sliced into 5 µm-thick sections, fixed on a glass slide and baked dry. The staining procedures were performed according to the instructions. Briefly, 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.
2.1.3 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.).
2.1.4 Methylated RNA immunoprecipitation sequencing (MeRIP-Seq) and library construction
MeRIP-Seq was performed following a previously reported procedure with slight modifications. Briefly, m6A RNA immunoprecipitation was performed with the GenSeqTM m6A RNA IP Kit (GenSeq Inc., China) by following the manufacturer’s instructions. Both the input sample without immunoprecipitation and the m6A 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.
2.2 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 profiles from 18 IPAH patients. The median expression values among all multiple probe IDs were selected to represent the corresponding gene symbol, leading to the identification of 18612 unique genes across 18 samples. Human gene annotation file (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.
2.3 Estimation of m6A methylation score and implementation of weight correlation network analysis (WGCNA)
We constructed gene signatures of m6A writer (METTL3, METTL14, METTL16, WTAP, KIAA1429, ZC3H13, RBM15/RBM15B), and m6A eraser (ALKBH5, FTO) as suggested by Yang et.al (25). We then computed the Gene Set Variation Analysis (GSVA) enrichment score of m6A writer and m6A eraser across the 18 samples using the ‘GSVA’ package(26) in R software. The estimated m6A methylation score was calculated by substrating m6A eraser score from m6A 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 first principal component of all gene expression levels in the module, and therefore, it was reasonable to consider that MEs represented all genes within a specific module. According to Pearson's correlation tests, we further identified the association between MEs and m6A methylation score. Within the most relavant module, those positively correlative with m6A methylation score were made subjected to further analysis.
2.4 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 m6A methylation correlated genes identified from WGCNA and upmethylated genes identified from m6A-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 m6A 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 first applies enrichment analysis to each gene list individually and identifies 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. Briefly, 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 significant (lowest P-value) term within each of them to represent it.
2.5 Protein-protein interaction (PPI) network construction of m6A methylation correlated genes
To find out the functional associations among the m6A 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 significant clusters with strong protein-protein linkages with default parameters.
2.6 Estimation of immunocyte infiltration
The CIBERSORT algorithm (https://cibersort.stanford.edu/) was applied to estimate the abundance of infiltrated 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 m6A methylation’ and ‘Low m6A methylation’ by the median value of their m6A methylation score. Comparisons between two groups were tested by the Wilcox rank test. P-value < 0.05 was considered significant.
2.7 Statistical analysis
Briefly, 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 identified by MACS software. Differentially methylated sites with a fold change cutoff of ≥ 2 and false discovery rate cutoff of ≤ 0.0001 were identified with the diffReps differential analysis package. These peaks identified by both softwares overlapping with exons of mRNA were figured out and choosed by home-made scripts.