Comprehensive analysis of m⁶A methylomes in idiopathic pulmonary arterial hypertension

ABSTRACT

Idiopathic pulmonary arterial hypertension (IPAH) is a serious and fatal disease. Recently, m6A has been reported to play an important role in the lungs of IPAH patients and experimental pulmonary hypertension models. However, the meaning of m6A mRNAs in the peripheral blood of IPAH patients remains largely unexplored. We aimed to construct a transcriptome-wide map of m6A mRNAs in the peripheral blood of IPAH patients. M6A RNA Methylation Quantification Kit was utilized to measure the total m6A levels in the peripheral blood of IPAH patients. A combination of MeRIP-seq, RNA-seq and bioinformatics analysis was utilized to select m6A-modified hub genes of IPAH. MeRIP-qPCR and RT-qPCR were used to measure the m6A levels and mRNA levels of TP53, RPS27A, SMAD3 and FoxO3 in IPAH patients. Western blot was performed to assess the protein levels of m6A related regulators and m6A related genes in experimental PH animal models, hypoxia-treated and PDGF-BB induced PASMCs. We found that the total m6A levels were increased in peripheral blood of IPAH patients and verified that m6A levels of RPS27A and SMAD3 were significantly elevated and m6A levels of TP53 and FoxO3 were significantly reduced. The mRNA or protein levels of RPS27A, SMAD3, TP53 and FoxO3 were changed in human blood samples, experimental PH animal models and PDGF-BB induced PASMCs. Moreover, METTL3 and YTHDF1 were increased in the hypoxia induced pulmonary hypertension rat model, hypoxia-treated and PDGF-BB induced PASMCs. These finding suggested that m6A may play an important role in IPAH.

KEYWORDS:

• N6-methyladenosine
• idiopathic pulmonary arterial hypertension
• METTL3
• YTHDF1
• TP53
• RPS27A
• FoxO3
• SMAD3

Introduction

Pulmonary arterial hypertension (PAH) is a complex vascular disease characterized by elevated pulmonary vascular resistance due to vasoconstriction or pulmonary vascular remodeling, with a low survival and high disability rate [1,2]. According to the clinical and pathophysiological characteristics of patients, the World Health Organization (WHO) classifies pulmonary hypertension (PH) into five groups [3]. Fibrosis, vascular hypertrophy and reticular lesions are the main pathologic features of PAH. The term idiopathic pulmonary arterial hypertension (IPAH) is used to categorize patients with pre-capillary pulmonary hypertension of unknown pathogenic factors. The study found that the average life expectancy of patients with IPAH after targeted therapy was 2.8 years, while the 1 year, 3 years and 5 years survival rates of IPAH patients were only 68%, 48% and 35%, respectively, suggesting a poor prognosis [4]. Existing IPAH targeted therapy drugs mainly reduce pulmonary circulation pressure through vasodilation. The progressive increase of pulmonary circulation resistance in the end stage of IPAH is mainly caused by pulmonary vascular remodeling, but the mechanism is not clear, resulting in the lack of specific treatment for IPAH [2,5]. Therefore, despite the development of targeted drugs for pulmonary arterial hypertension in the past 20 years, the prognosis of IPAH is still poor. Hence, there is a pressing need to identify and discover new therapeutic targets for IPAH.

In the 1970s, N6-methyladenosine (m⁶A) was found to be the most significant modification of polyadenylated mRNAs and long noncoding RNAs in higher eukaryotes [6]. m⁶A modification mainly interacts with three types of proteins. The first type is writers, such as methyltransferase-like 3 (METTL3)/methyltransferase-like 14 (methyltransferase-like 14). METTL14 complex and The Wilms tumor 1-associated protein (WTAP) can catalyze the formation of mRNA m⁶A. The second type is erasers that can demethylate mRNA m⁶A, including fat mass and obesity-associated protein (FTO) and AlkB homolog 5 (ALKBH5) [7]. The third type is readers, m⁶A also indirectly affects RNA processes by recruiting specific readers. Several m⁶A reader proteins have been identified in mammalian cell extracts that can directly or indirectly bind to the m⁶A methylation site in the mRNA, such as YTHDF1/2/3, YTHDC1/2, IGF2BP1–3, to perform its function. The combination of YTHDF1 and eIF3 with m⁶A can improve the efficiency of mRNA translation [8], while YTHDF2 and m⁶A combine to promote mRNA degradation [9]. IGF2BP1–3 binds to m⁶A to maintain mRNA stability [10.] Above all, m⁶A-related proteins have been rediscovered to play a significant role in multiple aspects of RNA metabolism, including mRNA stability and splicing, translation efficiency, nuclear export, alternative polyadenylation, as well as microRNA processing [11].

RNA m⁶A modification has recently been reported to play a regulatory role in human cancer [12]. Although IPAH is not a cancer, the similarities between these two disorders suggest that lessons learned from cancer may provide important insights into PAH pathogenesis and thus open new therapeutic avenues [13]. Interestingly, studies have identified transcriptomics of m⁶A circRNAs in hypoxia-mediated pulmonary hypertension and found that m⁶A levels in circRNAs decreased under hypoxia [14]. Elevated m⁶A levels were observed in lung of IPAH patients and experimental PH models [15]. But the change of m⁶A levels in peripheral blood of IPAH patients remains largely unexplored. In this study, we measured m⁶A levels and obtained the transcription profile of IPAH patients through sequencing results analysis.

Materials and methods

Subjects

Eight patients diagnosed with idiopathic pulmonary arterial hypertension between August 2013 and May 2022 were retrospectively collected. The diagnostic criteria for IPAH were the mean pulmonary artery pressure (mPAP) ≥25 mmHg, pulmonary artery wedge pressure (PAWP) ≤ 15 mmHg and pulmonary vascular resistance > 3 wood units as measured by standard right heart catheterization. At the same time, pulmonary hypertension caused by congenital heart disease, connective tissue disease, pulmonary embolism, primary pulmonary disease and other causes were excluded. In the healthy control group, pulmonary arterial hypertension was excluded by echocardiography.

Collection of human blood samples and establishment of experimental PH models

Blood samples were taken from both eight IPAH patients and five healthy controls. The clinical information of healthy controls and IPAH patients is shown in Table A1. Male Sprague – Dawley (SD) rats (weighing 180–200 g) were obtained from the Hunan Slake Jingda Experimental Animal company (Changsha, China). The animal model was divided into two parts. In the first part, the rats were randomly divided into a normoxic group and a hypoxic group (six rats per group). The rats in the hypoxia group were kept in a hypoxia chamber at 10% FiO₂ for 4 weeks 16], while the normoxic group were maintained under normal pressure and oxygen (FiO₂ = 21%). In the second part, rats in the MCT group were intraperitoneally injected with monocrotaline (60 mg/kg) on the first day to induce PAH. The control group was treated with normal saline at the same time on the first day (three rats per group).

Hemodynamic test

Four weeks later, the right external jugular vein dissection was performed to allow the distal end of the catheter to enter the right atrium and the ventricle. Right ventricular systolic pressure (RVSP) was recorded by a pressure sensor connected to the proximal end of the catheter [17].

Rats under deep anesthesia were sacrificed by cervical dislocation. The heart tissue was harvested and used to evaluate the degree of right ventricular hypertrophy [The weight of the right ventricle/the left ventricle plus the ventricular septum, RV/(LV+S)]. The lung tissue was preserved in liquid nitrogen or 4% buffered paraformaldehyde solution for further experiments. All experimental methods related to animals were approved by the Institutional Animal Care and Use Committee of Hunan Children’s Hospital.

Cell isolation, culture and treatment

Pulmonary arterial smooth muscle cells (PASMCs) were isolated from the pulmonary arteries of 6-week-old male SD rats. Briefly, lungs were collected, and all the pulmonary arteries were separated, removing the connective tissues and endothelium carefully 16], and then placed into PBS and Serum-free Dulbecco’s Modified Eagle medium for cleaning. PASMCs were cultured in Dulbecco’s modified Eagle’s medium supplemented with 15% fetal bovine serum (FBS) and incubated in a humidified incubator at 37°C with 21% O₂ and 5% CO₂. When PASMCs grown to 70–80% confluence, cells were starved for 24 h with serum free DMEM/F12 medium. For the normoxia group, cells were cultured in 21% O₂ and 5% CO₂ for 24 h. For the hypoxia group, cells were maintained at 1% O₂, 5% CO₂% and 94% N₂ for 24 h. Cells were treated with PDGF-BB (40 ng/ml) for 24 h and the control group was no PDGF-BB. PASMCs between passages 3 and 5 were used for the experiments.

Measurement of m⁶A modification

RNA was extracted from human peripheral blood samples by TRIzol LS Reagent (Invitrogen life technologies, America) following the manufacturer’s protocol. Next, Total RNA was assessed for m⁶A content using EpiQuik ™ m⁶A RNA Methylation Quantification Kit (Epigentek, America). The specific procedure was as follows: the binding solution and 200 ng RNA were added to the detection plate in turn. In this assay, positive control and negative control were set simultaneously. Then, capture antibody solution, detection antibody solution and enhancer solution were added and incubated for the corresponding time. The m⁶A level was calculated from the absorbance value at 450 nm on the microplate reader.

MeRIP enrichment and library quality control

The RNA concentration in each human peripheral blood sample was measured using Nanodrop ND-1000 (Thermo Fisher Scientific, Waltham, MA, USA). The OD260/OD280 values were used as RNA purity indexes. If the OD260/OD280 values range from 1.8 to 2.1, RNA purity is qualified. RNA integrity and gDNA contamination were then measured using denatured agarose gel electrophoresis. The library quality was checked with an Agilent 2100 BioAnalyzer.

MeRIP library preparation and sequencing

The MeRIP-seq technique was developed on the basis of published experimental methods 18]. Briefly, fragmented RNA was incubated with anti-m⁶A polyclonal antibody (Synaptic Systems, 202003) in IPP buffer for 2 hours at 4°C. The reaction mixture was further immunoprecipitated with Protein A magnetic beads (Thermo Fisher) at 4°C for 2 hours. Then, bound RNA was eluted from the beads with N6-methyladenosine (BERRY & ASSOCIATES, PR3732) in IPP buffer and then extracted with Trizol reagent (Thermo Fisher) according to the manufacturer’s instruction. Purified RNA was used for RNA-seq library generation with the NEBNext® Ultra™ RNA Library Prep Kit (NEB). Both the input sample without immunoprecipitation and the m⁶A IP samples were subjected to 150 bp paired-end sequencing on Illumina HiSeq sequencer.

Preparation and sequencing of RNA libraries

The high throughput sequencing service was provided by CloudSeq Biotech (Shanghai, China). Total RNA (1 μg) was used for removing the rRNAs using Ribo-Zero rRNA Removal Kits (Illumina, San Diego, CA, USA) following the manufacturer’s instructions. RNA libraries were constructed by using rRNA-depleted RNAs with TruSeq Stranded Total RNA Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. Libraries were controlled for quality and quantified using the BioAnalyzer 2100 system (Agilent Technologies, Inc., USA). Ten pM libraries were denatured as single-stranded DNA molecules, captured on Illumina flow cells, amplified in situ as clusters and finally sequenced for 150 cycles on Illumina HiSeq Sequencer according to the manufacturer’s instructions [19].

m⁶A methylated RNA immunoprecipitation sequencing (MeRIP) qPCR

To examine m⁶A modifications of individual genes in each human peripheral blood sample, total RNA was extracted using GenSeq® m⁶A MeRIP Kit (CloudSeq, Shanghai, China). RNA was incubated overnight with A/G immunomagnetic beads and m⁶A antibody premix at 4°C. Next, he binding RNA was eluted from the magnetic beads by adding RLT buffer for the next step of the experiment. The expression of SMAD3, TP53, FoxO3 and RPS27A were analyzed by qPCR. The primer sequences are provided in Table A2.

RT-qPCR

To determine the expression of target genes, total RNA was extracted and isolated from human peripheral blood samples according to the TRIzol kit instructions (Servicebio, Wuhan, China). After that, the isolated RNA was reverse transcribed to cDNA using a reverse transcription kit (Thermo Fisher Scientific, MA, United States). Quantitative real-time RT-PCR was performed using the FastStart Universal SYBR Green Master Mix (Rox) (Roche, Switzerland). GAPDH was used as the reference. The relevant primer information is listed in Table A3.

Histological analysis

After the distal lung tissues were excised, they were fixed in 4% paraformaldehyde overnight, placed in paraffin blocks, and then stained with hematoxylin and eosin. After that, the thickness of the pulmonary artery medium could be measured. Images of distal pulmonary arterioles for each rat (diameter between 30 and 100 μm) were captured using a light microscope. The pulmonary artery media thickness (PAMT) is defined as the distance between inner and outer elastic lamina. Vessel external diameter (ED) was determined. The relative PAMT (%) was calculated as 100 × 2PAMT/ED [20].

Western blot

Lung tissues and PASMCs were lysed with RIPA (Beyotime, China) and PMSF (Beyotime, China). After measuring the protein concentration of these samples, they were loaded on 10% and 12.5% SDS – PAGE gels and then transferred from gels to a PVDF membrane. Next, the membranes were incubated in blocking solution, namely 5% non-fat milk, for 1 hour and incubated with primary antibodies overnight. Membranes were trimmed according to the predicted molecular weight of the antibody. Finally, after membranes being incubated in the solutions with secondary antibody for 1 hour and being photographed via developer, the band intensities were analyzed using Image J software. Relative densitometry refers to β-actin levels. All the blots were cut prior to hybridization with antibodies during blotting.

Immunofluorescence staining

Paraffin sections were submerged in EDTA antigen repair buffer (pH 9.0) in a repair box to retrieve antigen using a microwave oven. The primary antibody was added to the slice and incubated overnight at 4°C in a wet box. The sections were washed three times and dried slightly. Then, the fluorescence-conjugated secondary antibody against the primary antibody was added to the section to fully submerge the tissue. The sections were incubated at room temperature for 50 min, followed by DAPI nuclear counterstaining and washing. After sealing with an anti-fluorescence quenching tablet, the sections were observed under laser scanning confocal microscopy.

Data analysis

Paired-end reads were harvested from Illumina HiSeq 4000 sequencer and were quality controlled by Q30. After 3’ adaptor-trimming and low-quality reads removal, cutadapt software (v1.9.3) was used 21]. First, clean reads of all libraries were aligned to the reference genome (UCSC hg19) by Hisat2 software (v2.0.4) [22. Methylated sites on RNAs (peaks) were identified by MACS software 23]. Diffreps software was used to identify differentially methylated genes (P value ≤ 0.000001 and Foldchange ≥ 4) 24]. Then, under the guidance of the Ensembl GTF gene annotation file, cuffdiff software [19] (part of cufflinks) was used to obtain the gene level Fragments per kilobase of exon per million fragments mapped (FPKM) as the expression profiles of mRNA, fold change and p-value were calculated based on FPKM. Fold change ≥ 2, P<0.05 and FPKM ≥ 0.1 at least one sample as a threshold for difference screening.

The Gene Ontology (GO) project developed a structured, controlled vocabulary for annotated genes, Gene products, and sequences in communities. It is divided into three parts: Molecular Function (MF), Biological Process (BP) and Cell Component (CC). GO functional analysis was performed using differentially expressed mRNA to annotate and infer the function of these mRNA. P value ≤ 0.05 is considered statistically significant.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis is a process in which the molecular data sets in genomics, transcriptomics, proteomics and metabolomics are mapped to the KEGG pathways graph, in order to interpret their biological functions [25–27]. Pathway analysis of differentially expressed mRNA can indicate the pathways in which they participate and their biological functions. P value ≤ 0.05 was used as the threshold of significant enrichment.

Student’s t-test (two-tailed) was used for statistical significance. P value < 0.05 was considered statistically significant. All experiments were independently repeated at least three times.

Result

The transcriptome-wide map of m⁶A mRNAs was constructed in the peripheral blood of IPAH patients

We found that the m⁶A levels of total RNA were increased in peripheral blood of IPAH patients compared with the levels in healthy controls (Figure 1a) and took advantage of MeRIP-seq to map the m⁶A modification in healthy controls and IPAH patients. In total, we identified 2732 m⁶A sites in healthy controls and 2079 m⁶A sites in IPAH patients, 2197 m⁶A sites were detected within both normal and IPAH patients (Figure 1b). As Table 1 displayed, compared with healthy controls, 506 differentially upmethylated m⁶A sites and 2170 differentially downmethylated m⁶A sites were selected in IPAH patients. To recover the distribution of differentially methylated m⁶A sites (DMMSs) in transcript, we systematically classified m⁶A sites into five transcript regions-5’UTR, 3’UTR, stop code, start code and coding sequences (CDs) and found that they were mainly enriched in 3’UTR and CDs (Figure 1c), which was consistent with previous reports. Altered m⁶A sites were discovered on all chromosomes, chr1 showed the greatest number of DMMSs, while chr19 showed the highest density of DMMSs altered (Figures 1d,e).

Figure 1. Construction a transcriptome-wide map of m⁶A mRnas identified in the peripheral blood of healthy controls and IPAH patients. (a) the m⁶A content in total RNA form the peripheral blood of CON group and IPAH group. (b) Venn diagram showing the number of methylated m⁶A sites in CON and IPAH. (c) Pie chart showing the proportion of differentially methylated m⁶A sites in the five transcriptional sequences-5’UTR, 3’UTR, stop codon, start codon, CDs. (d) Bar chart showing the distribution of differentially methylated m⁶A sites on chromosomes 1–22 and sex chromosomes. (e) Relative density of DMMSs in each chromosome normalized by length of the respective chromosome. Orange, upmethylated m⁶A sites in IPAH patients; Green, downmethylated m⁶A sites in IPAH patients.

Table 1. Total numbers of differentially methylated N6-methyladenosine sites and associated genes.

Item	upmethylated site	upmethylated gene	downmethylated site	downmethylated gene
mRNA	506	479	2170	1858

GO analysis and KEGG pathway of DMMSs

A total of 2337 protein coding genes containing DMMSs were detected in IPAH patients, including 479 differentially upmethylated genes and 1858 differentially downmethylated genes (Table 1). To discover the role of these genes containing DMMSs, GO and KEGG pathway analyses were performed. For the BP category, we observed that differentially downmethylated genes were highly enriched in SMAD protein complex assembly, positive regulation of histone deacetylation, mRNA destabilization, cytoskeletal anchoring at plasma membrane (Figure 2a & Table 2), while differentially upmethylated genes were significantly enriched in microtubule anchoring (Figure 2b & Table 3). For the MF category, differentially downmethylated genes were mainly enriched in protein kinase B binding, co-SMAD binding, cytoskeletal adaptor activity (Figure 2a & Table 4), while differentially upmethylated genes were mainly enriched in oxidized DNA binding (Figure 2b & Table 5).

Figure 2. Function annotations of genes containing DMMSs. (a) Top 10 biological processes and molecular function of differentially downmethylated genes. (b) Top 10 biological processes and molecular function of differentially upmethylated genes. (c) Top 10 pathway analysis of differentially downmethylated genes. (d) Top 10 pathway analysis of differentially upmethylated genes.

Table 2. Biological process of downregulated m⁶A methylation.

GO.ID	Term	Count	Fold. Enrichment	P-value	FDR
0097201	negative regulation of transcription from RNA polymerase II promoter in response to stress	6	5.431692308	0.000341824	0.011978003
0045945	positive regulation of transcription by RNA polymerase III	5	4.937902098	0.001889292	0.04367133
0007183	SMAD protein complex assembly	6	4.655736264	0.00094497	0.025665989
0097284	hepatocyte apoptotic process	5	4.526410256	0.002994627	0.060818804
0043922	negative regulation by host of viral transcription	5	4.526410256	0.002994627	0.060818804
0010457	centriole-centriole cohesion	5	4.526410256	0.002994627	0.060818804
0031065	positive regulation of histone deacetylation	7	4.473158371	0.000469004	0.015333207
0061157	mRNA destabilization	9	4.444111888	7.6092E–05	0.003477097
0033148	positive regulation of intracellular estrogen receptor signaling pathway	4	4.345353846	0.009544623	0.138817811
0031848	protection from non-homologous end joining at telomere	4	4.345353846	0.009544623	0.138817811
0031053	primary miRNA processing	4	4.345353846	0.009544623	0.138817811
0007016	cytoskeletal anchoring at plasma membrane	4	4.345353846	0.009544623	0.138817811
0002664	regulation of T cell tolerance induction	4	4.345353846	0.009544623	0.138817811

‘Count’ stands for the number of DE genes associated with the listed GOID.

Table 3. Biological process of upregulated m⁶A methylation.

GO.ID	Term	Count	Fold. Enrichment	P-value	FDR
0034453	microtubule anchoring	6	10.23362319	1.94662E–05	0.003761417
0045063	T-helper 1 cell differentiation	4	10.03296391	0.000556815	0.067245369
0034134	toll-like receptor 2 signaling pathway	3	9.840022297	0.003074075	0.189130533
0045625	regulation of T-helper 1 cell differentiation	2	8.528019324	0.021798932	0.580418014
0042908	xenobiotic transport	2	8.528019324	0.021798932	0.580418014
0038003	opioid receptor signaling pathway	2	8.528019324	0.021798932	0.580418014
0034058	endosomal vesicle fusion	2	8.528019324	0.021798932	0.580418014
0032836	glomerular basement membrane development	2	8.528019324	0.021798932	0.580418014
0014894	response to denervation involved in regulation of muscle adaptation	2	8.528019324	0.021798932	0.580418014
0014877	response to muscle inactivity involved in regulation of muscle adaptation	2	8.528019324	0.021798932	0.580418014
0003356	regulation of cilium beat frequency	2	8.528019324	0.021798932	0.580418014

‘Count’ stands for the number of DE genes associated with the listed GOID.

Table 4. Molecular function of downregulated m⁶A methylation.

GO.ID	Term	Count	Fold. Enrichment	P-value	FDR
0098505	G-rich strand telomeric DNA binding	5	5.336982968	0.001208508	0.028476835
0043422	protein kinase B binding	6	5.336982968	0.000376501	0.010843234
0070530	K63-linked polyubiquitin modification-dependent protein binding	9	4.57455683	5.62867E–05	0.002026321
0043047	single-stranded telomeric DNA binding	5	4.447485807	0.003237143	0.060753992
0008641	ubiquitin-like modifier activating enzyme activity	5	4.269586375	0.002994627	0.060818804
0001206	transcriptional repressor activity, RNA polymerase II distal enhancer sequence-specific binding	4	3.881442159	0.014790908	0.19169017
0098847	sequence-specific single stranded DNA binding	5	3.812130692	0.006980259	0.11128415
0001222	transcription corepressor binding	5	3.812130692	0.006980259	0.11128415
0008536	Ran GTPase binding	12	3.767282095	3.39091E–05	0.001292537
0070410	co-SMAD binding	4	3.557988646	0.020567791	0.220954474
0042975	peroxisome proliferator activated receptor binding	4	3.557988646	0.020567791	0.220954474
0008093	cytoskeletal adaptor activity	4	3.557988646	0.020567791	0.220954474

Count’ stands for the number of DE genes associated with the listed GOID.

Table 5. Molecular function of upregulated m⁶A methylation.

GO.ID	Term	Count	Fold. Enrichment	P-value	FDR
0070840	dynein complex binding	5	9.061241351	0.000182742	0.032471479
0032356	oxidized DNA binding	2	8.336342043	0.022748279	0.737044237
0017166	vinculin binding	2	8.336342043	0.022748279	0.737044237
0098748	endocytic adaptor activity	2	7.578492766	0.027367914	0.738933667
0035615	clathrin adaptor activity	2	7.578492766	0.027367914	0.738933667
0000774	adenyl-nucleotide exchange factor activity	2	7.578492766	0.027367914	0.738933667
0043422	protein kinase B binding	2	6.946951702	0.0323281	0.805715735
0008093	cytoskeletal adaptor activity	2	6.946951702	0.0323281	0.805715735
0019789	SUMO transferase activity	3	6.946951702	0.008557751	0.553923826
0003841	1-acylglycerol-3-phosphate O-acyltransferase activity	3	6.581322665	0.009984023	0.553923826

‘Count’ stands for the number of DE genes associated with the listed GOID.

Regarding KEGG pathway analysis, we found that differentially downregulated genes were enriched in the NF-κB signaling pathway, Toll-like receptor signaling pathway, MAPK signaling pathway, FoxO signaling pathway, TGFβ signaling pathway (Figure 2c & Table 6). Differentially upregulated genes were enriched in the NF-κB signaling pathway, Wnt signaling pathway, AMPK signaling pathway and NOD-like receptor signaling pathway (Figure 2d & Table 7). Most of these pathways were associated with pulmonary vascular remodeling.

Table 6. KEGG pathway analysis of downregulated m⁶A methylation.

Pathway.ID	Definition	Count	P-value	FDR
hsa04064	NF-κB signaling pathway	20	8.41253E–05	0.007949932
hsa04620	Toll-like receptor signaling pathway	20	0.000536203	0.01729145
hsa04010	MAPK signaling pathway	36	0.001646395	0.03186836
hsa04621	NOD-like receptor signaling pathway	12	0.002627967	0.0394195
hsa04068	Foxo signaling pathway	19	0.01853444	0.1426183
hsa04350	TGFβ signaling pathway	13	0.02468286	0.1609752
hsa04330	Notch signaling pathway	8	0.04808801	0.2326839

‘Count’ stands for the number of differentially methylated genes associated with the listed Pathway.ID.

Table 7. KEGG pathway analysis of upregulated m⁶A methylation.

Pathway.ID	Term	Count	P-value	FDR
hsa04668	TNF signaling pathway	8	0.002409407	0.053333
hsa04064	NF-κB receptor signaling pathway	6	0.01457591	1
hsa04152	AMPK signaling pathway	7	0.01698368	1
hsa04310	Wnt signaling pathway	7	0.03271221	1
hsa04068	NOD-like receptor signaling pathway	4	0.03344811	1

‘Count’ stands for the number of differentially methylated genes associated with the listed Pathway.ID.

Conjoint analysis of DMMSs and DGEs

We identified 503 upregulated differentially expressed genes (DEGs) and 1715 downregulated differentially expressed genes (Figure 3a). Conjoint analysis with protein-coding genes containing DMMSs showed that these genes could be divided into four plates, including 57 upmethylated as well as downregulated genes (methy.up-expression.down), 31 upmethylated and upregulated genes (methy.up-expression.up), 408 downmethylated and downregulated genes (methy.down-expression.down), and 52 downmethylated as well as upregulated genes (methy.down-expression.up) (Figure 3b).

Figure 3. Conjoint analysis of differentially methylated and differentially expressed genes. (a) Scatter diagram showing the number of differentially expressed genes (fold change ≥ 2, P<0.05 and FPKM ≥ 0.1 at least one sample). (b) Venn diagram showing the differentially expressed genes containing DMMSs. (c) Ten hub genes identified by cytoHubba in the development of IPAH. (d) m⁶A enrichment of SMAD3, TP53, RPS27A and FoxO3 mRNA in the peripheral blood of CON group and IPAH group. (e) the mRNA expression levels of SMAD3, TP53, RPS27A and FoxO3 were measured in the peripheral blood of CON group and IPAH group. *P<0.05 (different from CON),**P<0.01 (different from CON), ****P<0.0001 (different from CON).

Protein-protein interaction network analysis

We imported 388 genes containing both DMMSs and DEGs into the STRING database for protein-protein interaction network analysis. Then, 125 nodes and 192 edges were screened out by setting the combined score to 0.9. At the same time, the Cytoscape 3.7.2 Cytohubba plug-in was used for network topology properties analysis. The Degree and Betweenness Centrality of nodes were used to describe the importance of network nodes and 10 hub genes were obtained according to the degree algorithm, including TP53, SMAD3, FoxO3, JUN, RPS27A, RUNX3, AGO3, PIAS4, CBX4 and CSNK1E (Figure 3c). The first three genes were found to be associated with the development of PAH. Next, we measured m⁶A levels and relative mRNA expression of some hub genes, such as TP53, SMAD3, FoxO3, and RPS27A, in the peripheral blood of IPAH patients. Elevated m⁶A levels of RPS27A and SMAD3 and reduced m⁶A levels of TP53 and FoxO3 were observed in IPAH patients. The mRNA levels of TP53 and SMAD3 were significantly decreased, and the mRNA levels of RPS27A was significantly increased in IPAH blood samples, while FoxO3 mRNA levels had no obvious change (Figures 3d,e).

Establishment of experimental PH models

To further explore the role of the m⁶A modification in pulmonary arterial hypertension, we constructed a hypoxia induced pulmonary hypertension (HPH) rat model and a monocrotaline (MCT) induced pulmonary hypertension rat model. The right ventricular systolic pressure (RVSP) in HPH group elevating to 54.98 ± 1.72 mmHg compared with 23.74 ± 3.43 mmHg in N group. The right cardiac hypertrophy index was evaluated by the ratio of right ventricle/left ventricle plus ventricular septum (RV/LV+S). The RV/LV+S ratio in HPH group was significantly higher than that in the N group (0.49 ± 0.04 vs 0.25 ± 0.03, P<0.0001) (Figure 4a). HE staining showed that the medium membrane of pulmonary small vessels in the lung tissue of HPH group was significantly thicker than that of the N group (71.67 ± 10.97% vs 38.67 ± 8.7%, P<0.001), and the vascular lumen was significantly narrower (Figures 4b,c). In addition, the proportions of fully muscularized vessels (35 ± 2.11% vs 20.33 ± 0.58%, P<0.01) and partially muscularized vessels (35.67 ± 1.15% vs 21 ± 1.11%, P<0.001) in HPH group were significantly higher than those in the N group, while the proportion of non-muscularized vessels (30.3 ± 4.51% vs 58.67 ± 1.53%, P<0.001) was significantly lower than that in N group (Figure 4d). Similarly, these changes were also observed in the MCT induced pulmonary hypertension rat model (Figure 5a–d). All these results indicated that HPH rat model and MCT induced pulmonary hypertension rat model were successfully established.

Figure 4. The identification of hypoxia-induced pulmonary arterial hypertension rat model. (a) RVSP and the ratio of RV/(LV + S) were measured in Normoxia group and Hypoxia group. (b) H&E staining was performed in the lung of Normoxia group and Hypoxia group. scale bar = 50 μm. Quantitative measurement of media thickness of pulmonary artery (c) and vessel muscularization (d). Data are present as means±SD (n = 6 rats in per group). *P<0.05 (different from Normoxia group), **P<0.01 (different from Normoxia group), ***P<0.001 (different from Normoxia group), ****P<0.001 (different from Normoxia group).

Figure 5. The identification of MCT-induced pulmonary arterial hypertension rat model. (a) RVSP and the ratio of RV/(LV + S) were measured in Control group and MCT group. (b) H&E staining was performed in the lung of Control group and MCT group. scale bar = 50 μm. Quantitative measurement of media thickness of pulmonary artery (c) and vessel muscularization (d). Data are present as means±SD (n = 3 rats per group). *P<0.05 (different from Control group), **P<0.01 (different from Control group), ***P<0.001 (different from Control group).

Upregulation of METTL3 and YTHDF1 in HPH rat model and hypoxia-treated and PDGF-BB induced pulmonary arterial smooth muscle cells

By using western blot, the protein expression of seven major m⁶A related regulators, METTL3, FTO, ALKBH5, IGF2BP2, YTHDF1, and YTHDF2 was verified between the HPH group and the N group. Our results showed that the expression of METTL3, IGF2BP2 and YTHDF1 were upregulated in the lungs of HPH rats, whereas the expression of ALKBH5 was downregulated, and YTHDF2 was not significantly upregulated (Figures 6a,b). Moreover, we analyzed the fluorescence intensity in the pulmonary vascular media by immunofluorescence, YTHDF1 was upregulated in the small pulmonary vessel while ALKBH5 was not significantly changed (Figure 7a–c). Similarly, the protein levels of METTL3, ALKBH5, IGF2BP2, YTHDF1 and YTHDF2 were measured to evaluate the effect of hypoxia treatment on the key m⁶A modulators. Consistent with the previous results in vivo, METTL3 and YTHDF1 were significantly upregulated, while FTO, ALKBH5, YTHDF2 and IGF2BP2 were not significantly changed in rat PASMCs exposed to hypoxia compared with rat PASMCs cultured in normoxia (Figures 8a,b). Furthermore, PASMCs exposed to PDGF-BB that contribute to the development of PAH also showed higher protein level of METTL3 and YTHDF1 and resulted in increased expression of proliferation protein PCNA (Figures 8c,d).

Figure 6. The expression of major m⁶A-related proteins in Hypoxia group and Normoxia group. (a) Western blot image of METTL3, FTO, ALKBH5, IGF2BP1, IGF2BP2, YTHDF1, YTHDF2. β-actin was used as an internal reference. (b) Statistical bar graph of METTL3, FTO, ALKBH5, IGF2BP1, IGF2BP2, YTHDF1, YTHDF2. Data are shown as means±SD (n = 6 rats per group). ^n.s.P>0.05 (different from Normoxia group), *P<0.05 (different from Normoxia group), **P<0.01 (different from Normoxia group), ***P<0.001 (different from Normoxia group).

Figure 7. Immunofluorescence staining of YTHDF1 (a), IGF2BP1 (b), ALKBH5 (c) in pulmonary vessels. α-SMA (red), objective protein (green) and nucleus were dyed with DAPI (blue), scale bar = 40 μm.

Figure 8. The expression of major m⁶A-related proteins in PASMCs exposed to hypoxia and PDGF-BB. (a) Western blot image of METTL3, FTO, ALKBH5, IGF2BP1, IGF2BP2, YTHDF1, YTHDF2, PCNA under hypoxia. β-actin was used as an internal reference. (b) Statistical bar graph of METTL3, FTO, ALKBH5, IGF2BP1, IGF2BP2, YTHDF1, YTHDF2, PCNA. (c) Western blot image of METTL3, YTHDF1, PCNA under PDGF-BB stimulation. β-actin was used as an internal reference. (d) Statistical bar graph of METTL3, YTHDF1, PCNA. Data are shown as means±SD (n = 3 each). ^n.s.P>0.05, *P<0.05, **P<0.01.

TP53 and RPS27A were increased in experimental PH models and PDGF-BB induced PAMSCs

Moreover, we measured the protein expression of TP53, RPS27A, FoxO3 and RPS27A in the lung tissue of HPH rats, MCT-induced pulmonary arterial hypertension rats and PDGF-BB induced PAMSCs proliferation model. A significant increase in TP53, RPS27A and FoxO3 protein levels and a significant decrease in SMAD3 protein levels were observed in experimental PH models. Elevated protein levels of RPS27A and reduced protein levels of SMAD3 in the lungs of HPH rats were observed while no obvious change in the lungs of MCT induced pulmonary arterial hypertension rats (Figure 9a–d). Upregulation of TP53, RPS27A and FoxO3 were also observed in PDGF-BB induced PAMSCs (Figures 9e,f).

Figure 9. The protein expression of m⁶A-regulated genes in lungs of experimental PH rat models and PDGF-BB induced PASMCs. Representative immunoblots of TP53, SMAD3, FoxO3 and RPS27A protein expression in (a) lung tissues of hypoxia-induced pulmonary arterial hypertension rats, (b) lung tissues of MCT-induced pulmonary arterial hypertension rats and (c) PDGF-BB induced PASMCs (n = 3 per group) normalized to β-actin. Relative densitometric analysis of TP53, SMAD3, FoxO3 and RPS27A protein expression in (d) lung tissues of hypoxia-induced pulmonary arterial hypertension rats, (e) lung tissues of MCT-induced pulmonary arterial hypertension rats and (f) PDGF-BB induced PASMCs (n = 3 per group) normalized to β-actin. Data are shown as means±SD. ^n.s.P>0.05, *P<0.05.

Discussion

A number of studies have shown that epigenetics occupied the main position in IPAH. As one of the most common post-transcriptional mRNA modifications in eukaryote, m⁶A plays an important role in circadian rhythm [28], stem cell reprogramming [29], cell proliferation [30] and other life processes. In recent years, increasing evidence has demonstrated that the disorder of m⁶A mRNAs can also be used as one of the markers of cardiovascular diseases such as ischemic heart disease [31], cardiac hypertrophy [32], heart failure [33], and pulmonary arterial hypertension [34–36]. A few studies have reported that m⁶A mRNAs are associated with IPAH [15]. However, it has not been reported whether m⁶A levels are also altered in the blood of IPAH patients. In this study, we first found that m⁶A levels of total RNA were increased in blood samples of IPAH patients compared with the m⁶A levels in healthy controls and identified 2676 differentially methylated m⁶A sites in IPAH patients. Meanwhile, we performed GO and KEGG pathway analysis to reveal the functions of these differentially methylated transcripts. Second, we filtered out 10 core hub gene of PAH combined with RNA-seq, MeRIP-seq and bioinformatics analysis. Moreover, we verified the m⁶A level and mRNA expression of TP53, SMAD3, FoxO3 and RPS27A by MeRIP-qPCR and RT-qPCR. Simultaneously, the protein levels of these hub genes were also detected in experimental PH models and PDGF-BB induced PSAMCs. Finally, we detected several m⁶A related regulators in HPH rat model, hypoxia-treated PASMCs and PDGF-BB induced PASMCs.

KEGG pathway analysis suggested that genes with upmethylated m⁶A sites were related to Wnt signaling pathway, which is closely associated with pulmonary vascular homeostasis. In recent years, targeted drugs that restore the Wnt signaling pathway to the physiological level have been considered to be a promising direction in the treatment of pulmonary arterial hypertension 36]. Moreover, genes with upmethylated m⁶A sites were also enriched in NF-κB signaling pathway and TNF-α signaling pathway. A series of studies have shown that these two pathways are associated with inflammation, and activation of these pathways promotes PASMCs proliferation [37–39. IL1β binds to IL1R1, thereby activating NF-κB to induce IL1, IL6 and TNFα synthesis, thus resulting in pulmonary vascular remodeling 40]. Our results suggested that m⁶A modification might suppress pulmonary vascular remodeling by regulating inflammation-related pathways.

Genes with downmethylated m⁶A sites were enriched in the MAPK signaling pathway. The MAPK signaling pathway is a critical mediator of hypertrophic stimuli and pro-inflammatory response in PAH [41,42,]. Importantly, PAH-associated BMPRII mutation eventually activated the MAPK pathway, and dysregulated MAPK signaling has been observed in various experimental models of PAH and in PASMCs from patients with BMPR2 mutations 43]. Nasim MT et al. found that BMPRII dysfunction enhanced the TGFβ-TAK1-MAPK signaling axis, resulting in an imbalance of smooth muscle cell proliferation and apoptosis 44]. These results suggested that m⁶A modification might be an intermediate hub between BMPRII gene mutation and MAPK signaling pathway. In addition, genes with downmethylated m⁶A sites were enriched in Foxo signaling pathway, which mainly plays an important role in pro-proliferative phenotype of PASMCs and oxidative stress 45].

By combining the analysis of MeRIP-seq and RNA-seq data, we discovered 388 genes with DMMSs and synchronously differential expression in IPAH. Based on the combination of biological informatics, we further filtered out 10 core target genes (TP53, SMAD3, FoxO3, JUN, RPS27A, RUNX3, AGO3, PIAS4, CBX4 and CSNK1E). It has been reported that some of these genes are associated with PAH. For example, TP53 is related to cellular processes such as cell proliferation, apoptosis, DNA repair, and autophagy which had gradually attracted attention in the field of pulmonary hypertension in recent years 46]. Mizuno et al. found that TP53 gene deficiency aggravated hypoxia induced pulmonary arterial hypertension and vascular remodeling in mice 47]. Wang Y et al found that MCT upregulated the expression of TP53 in PAH rat lung tissue [48]. However, it is difficult to utilize TP53 as a drug target because of its extensive downstream signaling pathway. Our results indicated that its mRNA was downmethylated and we also found that its mRNA levels were decreased while protein levels were increased in rodent PH models and PDGF-BB induced PASMCs, which undoubtedly expanded its role in the field of pulmonary hypertension research and provided a theoretical basis for discovering new drug targets. Loss of SMAD3 promoted human pulmonary artery smooth muscle (huPASMCs) and endothelial cells (huECs) proliferation and reconciled the induction of a contractile and synthetic smooth muscle cell phenotype via disinhibition of cardiac protein-related transcription factors MRTF to promote pulmonary arterial vascular remodeling 49]. In our study, we revealed that the m⁶A levels of SMAD3 was increased while relative mRNA expression and protein expression were decreased in the blood of IPAH patients and lung tissue of HPH rat. Forkhead box O3 (FoxO3) is a pro-apoptotic factor. Zhao et al. discovered that the expression of FoxO3 in the plasma of patients with PAH was decreased compared with healthy controls 50]. In addition, the activation of T-type Ca²⁺ channels in IPAH-derived pulmonary artery smooth muscle cells could stimulate the Akt signaling pathway, leading to a decrease in the expression of downstream signaling factor FoxO3, thus promoting the proliferation of pulmonary artery smooth muscle 51]. In our study, we found that m⁶A levels in FoxO3 were reduced and that relative protein levels were significantly decreased in rodent PH models while the opposite result was shown in PDGF-BB induced PASMCs proliferation model. RPS27A is not only involved in protein translation but also performs extra-ribosomal functions. Currently, there are no reports elucidating the role of RPS27A in PAH. Nevertheless, some studies have reported that RPS27A is associated with the proliferation and apoptosis of tumor cells 52]. We discovered that a significance increase m⁶A levels and mRNA levels of RPS27A were observed in IPAH patient samples, and higher RPS27A protein levels were observed in experimental PH models and PDGF-BB induced PASMCs. It is worth mentioning that the m⁶A abundance of RPS27A mRNA was also increased in the lung tissue of MCT induced pulmonary hypertension arterial rats. According to the protein-protein interaction diagram, RPS27A may play a certain role in PAH. These genes regulated by m⁶A modification may play an important role in the development of IPAH. However, the specific mechanism is still unclear.

To further clarify the mechanism by which methylation modification proteins regulate PAH, we constructed HPH rat models and measured the major m⁶A-related proteins in the lungs of rats undergone hypoxia or normoxia. Compared with the N group, the protein levels of METTL3, YTHDF1 and IGF2BP2 in the lungs of HPH rats were increased, while the protein levels of ALKBH5 were decreased. Most of these results are consistent with our research group’s previous study in the rat model of MCT-induced pulmonary hypertension 35]. To explore whether hypoxia regulates PASMCs proliferation by affecting methylation modifying proteins, rat pulmonary artery smooth muscle cells were isolated. Our results showed that the expression levels of PCNA, METTL3 and YTHDF1 were increased in PASMCs exposed to hypoxia compared with rat PASMCs cultured in normoxia. Moreover, when PASMCs were stimulated by PDGF-BB, the protein levels of METTL3 and YTHDF1 were also increased.

Previous reports indicated that m⁶A could influenced the stability of circRNAs, thus activating Wnt and Foxo signaling pathway to involve in the development of HPH [14]. Consistent with previous reports, our study found that the host genes of mRNA containing DMMS were enriched in the NF-κB, Wnt, Foxo and MAPK signaling pathway. Moreover, these DMMSs were mainly located in 3’UTR and CDs of mRNA which associated with gene expression. Until now, the relationship between m⁶A and mRNA expression has not been clear. Song H et al. discovered that METTL3 methylated TFEB at two m⁶A residues in the 3’UTR and decreased the expression of TFEB mRNA in H/R-treated cardiomyocytes and I/R-treated mouse heart [31]. In contrast, a number of studies suggested that the upregulation of methylation would increase mRNA expression. Niu Y et al. found that FTO increased the mRNA expression level of ALODA under hypoxia in liver cancer 53]. In our study, we revealed a significant increase in m⁶A levels in blood sample of IPAH patients, which was consistent with lung of IPAH patients and experimental PH models [15]. At present, the pathogenesis of IPAH is not clear and lack of specific blood biomarkers. Our findings suggest that the change of m⁶A levels is expected to be the potential novel molecular indicators for PAH diagnosis. An important thing is that we found discordant relationships between mRNA expression, m⁶A mRNA levels and protein expression, suggesting that at least one or several m⁶A related regulators are involved in the development of pulmonary arterial hypertension. This has been verified in vivo and in vitro. Elevated protein expression of YTHDF1 promote translation efficiency and IGF2BP2 regulated mRNA stability. The functional diversity of these m⁶A-related proteins further complicates the pathogenesis of pulmonary hypertension. In addition, perhaps other important m⁶A modified genes may not be detected as DGEs in RNA-seq but may alter the protein levels and even be important in PAH.

Our study was the first to demonstrate m⁶A modification levels and construct the transcriptome-wide map of m⁶A mRNAs in the peripheral blood of IPAH patients. Moreover, we identified that the m⁶A writer ‘METTL3’ and m⁶A reader ‘YTHDF1’ were significantly increased in HPH rat model, hypoxia-treated PASMCs and PDGF-BB induced PASMCs. YTHDF1 might be involved in the development of pulmonary arterial hypertension by binding to these DGEs containing DMMSs, such as SMAD3, TP53, FoxO3, and RPS27A, regulating translation to promoting PASMCs proliferation. In conclusion, our study provides a new direction for elucidating the new mechanism of PAH, m⁶A levels, m⁶A regulated genes and m⁶A related regulators may provide a reference basis for the diagnosis of PAH and disease course monitoring.

Limitation

However, there are some limitations in this study. Firstly, we only detected m⁶A related regulators in the animal and cell models but did not further detect changes in m⁶A related regulators in human sample and the sample size of our study was fairly limited. Secondly, we did not further confirm the specific mechanism of the involvement of these methylated target genes in pulmonary arterial hypertension.

Authors’ contributions

TH, YHZ, YBX and YCD contributed to design the research; TH, YHZ and HQF carried out the experiments; TH, YHZ, YY, WFL, JQL, FY and WJC analyzed the test data and made the figures; TH and YHZ drafted the manuscript; YBX helped polish the manuscript. All authors read and approved the final manuscript.

Availability of data and materials

The datasets generated and analyzed during the current study are available in the Sequence Read Archive (SRA) repository with accession number SRR19118209, SRR19118210, SRR19118211, SRR19118212, SRR19118213, https://www.ncbi.nlm.nih.gov/sra.

Ethics approval and consent to participate

The protocol of this research was submitted to and approved by the Institutional Animal Care and Ethics Committee of Hunan Children’s Hospital, China (permit number: HCHLL-2020-44) and written informed consent was obtained from all subjects and/or their legal guardian(s). All experimental methods related to animals were approved by the Institutional Animal Care and Use Committee of Hunan Children’s Hospital. We confirmed that all methods were performed in accordance with the relevant guidelines and regulations and that the study was reported in accordance with ARRIVE guidelines.

Acknowledgments

We thanked all subjects who participated in this study. We thank Cloud-Seq Biotech Ltd. Co. (Shanghai, China) for the MeRIP-Seq service and the subsequent bioinformatics analysis.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study was funded by the National Natural Science Foundation of China (81500041), Hunan Provincial Health Commission Project (20200483) and Hunan Clinical Research Center for Children’s Cardiovascular Diseases (2021SK4019).

Appendix

Table

	CON	IPAH
Age, years	27.00 ± 3.54	31.25 ± 12.70
Gender, male/female	2/3	4/4
BMI, kg/m2	22.92 ± 4.94	23.15 ± 3.29
WHO-Fc, I-II/III-IV	/	3/5
NT-proBNP, pg/ml	/	1090.47 ± 602.77
6MWT, m	/	460.88 ± 141.37
SPAP, mmHg	/	93.25 ± 23.04
mPAP, mmHg	/	59.13 ± 11.86
PAWP, mmHg	/	7.50 ± 4.12
PVR, Wood units	/	20.31 ± 7.42

BMI, body mass index; WHO-Fc, Word Health Organization functional class; NT-proBNP, n-terminal natriuretic peptide type-B; 6MWT, six-minute walking distance; SPAP, systolic pulmonary pressure; mPAP, mean pulmonary artery pressure; PAWP, pulmonary arterial wedge pressure; PVR, pulmonary vascular resistance.

Table A2. Primers used for MeRIP-qPCR.

Gene Name		Sequence (5’−3’)
SMAD3	F R	CAGTCCAACCAGAAACACCC ACTGTCACAAGCGTCAGAGA
TP53	F R	CAGTCAGATCCTAGCGTCGA CAGAACGTTGTTTTCAGGAAGT
RPS27A	F R	TTTCAACAAACCAGAAGACAAGT GCACTACTAAAAGGAAATGGTGT
FoxO3	F R	TTGCGATGGTTTATGGGACG GCTAGGGGAAGGACGGTTAA

Table A3. Primers used for Qrt-PCR.

Name	Primer sequence (5′-3′)	GenBank accession number	Product size (bp)
TP53	CCCTCCTCAGCATCTTATCCG GCACAAACACGCACCTCAAA	NM_001126115.2	261
RPS27A	TGAACCCTCGGATACGATAGAAA GCACCACCACGAAGTCTCAAC	NM_001135592.2	180
SMAD3	CTACCAGTTGACCCGAATGTGC TCTGTCTCCTGTACTCCGCTCC	NM_001145102.1	74
FoxO3	CGGACAAACGGCTCACTCT GCTCTTGCCAGTTCCCTCATT	NM_001455.4	173
GAPDH	GGAAGCTTGTCATCAATGGAAATC TGATGACCCTTTTGGCTCCC	NM_002046	168