Identification of immune‐associated signatures and potential therapeutic targets for pulmonary arterial hypertension

Abstract Pulmonary arterial hypertension (PAH) comprises a heterogeneous group of diseases with diverse aetiologies. It is characterized by increased pulmonary arterial pressure and right ventricular (RV) failure without specific drugs for treatment. Emerging evidence suggests that inflammation and autoimmune disorders are common features across all PAH phenotypes. This provides a novel idea to explore the characteristics of immunological disorders in PAH and identify immune‐related genes or biomarkers for specific anti‐remodelling regimens. In this study, we integrated three gene expression profiles and performed Gene Ontology (GO) and KEGG pathway analysis. CIBERSORT was utilized to estimate the abundance of tissue‐infiltrating immune cells in PAH. The PPI network and machine learning were constructed to identify immune‐related hub genes and then evaluate the relationship between hub genes and differential immune cells using ImmucellAI. Additionally, we implemented molecular docking to screen potential small‐molecule compounds based on the obtained genes. Our findings demonstrated the density and distribution of infiltrating CD4 T cells in PAH and identified four immune‐related genes (ROCK2, ATHL1, HSP90AA1 and ACTR2) as potential targets. We also listed 20 promising molecules, including TDI01953, pemetrexed acid and radotinib, for PAH treatment. These results provide a promising avenue for further research into immunological disorders in PAH and potential novel therapeutic targets.

agonists, endothelin receptor antagonists, phosphodiesterase inhibitors and guanylate cyclase agonists, focus on dilatation of the pulmonary arterial vasculature and reducing pulmonary vascular resistance, which can only ameliorate symptoms but cannot completely reverse pulmonary vascular remodelling and improve overall survival. 2,3Although current treatments for PAH can improve survival, they are not considered disease-modifying as the mean 5-year survival rate of PAH patients still remains at only 57%-59%. 4 The reason for this is that pathological mechanisms of inflammation and immune dysregulation have not been adequately addressed. 5As a result, there is an urgent need for more insights into the pathogenesis of PAH, immune dysregulation and the development of specific anti-remodelling regimens to explore new therapies.
7][8] In particular, it is now widely accepted that immunological disorders contribute to both disease susceptibility and the progression of vascular remodelling in PAH. 6There is a positive correlation between the degree of pulmonary perivascular inflammation, vascular intimal/medial/ adventitial thickness and mean pulmonary arterial pressure.Various inflammatory and immune markers are correlated with the hemodynamics or prognosis of PAH patients. 9,10These findings highlight that inflammation could be involved in pulmonary vascular remodelling and PAH development. 11Increasing evidence of the pathological role of immune cells in innate and adaptive immunity leads to many promising preclinical studies.IL-6-specific antagonist treatment reversed sugen/hypoxia-induced PAH and MCT-PAH in rat models. 12cell lymphopenia in association with vascular endothelial growth factor receptor 2 (VEGFR2) blockade resulted in periarteriolar inflammation with macrophages and B cells even prior to vascular remodelling and elevated pulmonary pressures.And immune reconstitution prevented early inflammation and attenuated PAH development. 13Activation of endogenous retrovirus sequences and their gene products could contribute to the chronic inflammatory and altered immune state related to the vascular remodelling in PAH.This activation may occur through activating monocytes, which leads to elevated production of cytokines such as TNFα, IL6 and IL1β, as well as SAMHD1.Additionally, this activation may also activate B cells, which could be responsible for the production of SAMHD1 antibodies in tertiary lymphoid tissue in the lungs. 5These preclinical studies, in turn, are contributing to innovative clinical trials.The TRANS-FORM-UK trial (NCT 02676947) has been completed, in which PAH patients are being treated with anti-IL-6, and 21 patients received the drug with good safety. 14A randomized multicentre placebocontrolled clinical trial testing B-cell depletion with rituximab to treat systemic sclerosis-associated PAH 15 has also been conducted.
To lead the way for specific targeted therapeutic approaches intervening in the affected homeostasis of innate and adaptive immunity, a deeper insight into the immune-related molecular mechanisms and interactions of the depicted contributors to the pathophysiology of PAH is essential.Gene sequencing, high-throughput technologies and bioinformatics methods are widely applied to screen new biomarkers or potential drug targets.For example, two new candidate genes-FBLN2 and PDGFD-with known functions in vasculogenesis and remodelling were identified by variant analysis by a large international consortium, and trio analysis predicted that 15% of paediatric IPAH may be explained by de novo variants. 16re, we investigated the signatures of tissue-infiltrating immune cells during the development of proliferative remodelling in PAH by integrating transcriptomic data of PAH lung tissues, followed by the validation of immune cell abundance and those genes in an experimental PAH model, which would enrich our understanding of the roles of inflammation and immune mechanisms in PAH progression and help to discover new potential biomarkers and treatment targets.

| Data collecting and processing
Gene expression profiles of human lung tissue with or without PAH were acquired from the GEO database (https://www.ncbi.nlm.nih.gov/geo/).GSE53408, GSE117261 and GSE15197, containing 96 samples of PAH and 49 normal controls, were included in this study.
To reduce the bias and variability inherent in these three different sequencing results, we integrated and processed the three data sets using the sva package in the R language, which supports using the sva function of proxy variable estimation, the ComBat function of adjusting known batch effects directly and the Fsva function of adjusting batch and latent variables in predictive problems.The limma Bioconductor package (R version 3.5.1)was used to make class comparison analysis for DEGs between PAH samples and normal samples, and DEGs with a threshold of p < 0.05 and |log2FC| >1 were selected for subsequent analysis.

| GO and KEGG pathway enrichment
The Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) was used for GO biological pathway and cellular component enrichment. 17The GO plot package was used to combine and integrate the expression data with the results of the DAVID analysis.KEGG pathway annotation and analysis of DEGs, which can determine the major metabolic and signal transduction pathways involved in these genes, were analysed by Cluster profiler in R software, performing statistical analysis and visualization of functional clustering of gene collections. 18p values <0.05 were considered statistically significant.

| Evaluation of tissue-infiltrating immune cells
CIBERSORT is a bioinformatics tool that uses the deconvolution method to characterize cell subsets of interest in high-dimensional genomic data derived from bulk tissue samples, which provides an estimate of the abundance of member cell types in a mixed cell population using gene expression data.It is used to estimate the proportion of 22 types of tissue-infiltrating immune cells (TIICs) in each sample, such as naive B cells, plasma cells, CD4+ resting memory T cells, Tregs, activated natural killer cells, monocytes, M0 macrophages, resting dendritic cells, activated mast cells, eosinophils and neutrophils. 19,20Processed GEPs were transformed into an abundance of 22 TIICs.The relative expression of 22 TIICs in each sample was determined.Significant results with a threshold of p < 0.05 were selected for further analysis.

| Construction of the PPI network
In the present study, the PPI network of DEGs was constructed by using the Retrieval of Interacting Genes (STRING) database (http:// strin g-db.org/).Analysing the functional interactions between genes may provide insights into the mechanisms of PAH development.
Herein, the PPI network of DEGs was built with a threshold interaction score of 0.4, which was visualized by the Cytoscape software (version 3.7.2).Cytoscape software is a very formidable tool for visualizing the interrelationships among a set of genes or proteins.The hub genes were screened by PPI and the CytoHubba plugin.

| Classification by machine learning
Pearson correlation analysis was used to wipe off the highautocorrelation and low correlation between the phenotypes of normal tissue and PAH tissue and gene signatures.The phenotypes of normal tissue and PAH tissue were set as '0' and '1', respectively.The correlation coefficient <0.1 was eliminated, and if the autocorrelation was more than 0.9, the gene signature with a lower correlation was removed.Then the stepwise regression approach was further used to filter the remaining gene signatures, which considered variable size, significance and contribution.Every new regression equation was employed with a significance test to assess the addition of each new gene signature.When there were no novel signatures inputted or removed, the process was terminated.After filtration, the selected descriptors served as signatures to establish the prediction models.The normal and PAH tissues were divided into training and test groups according to a ratio of 3:1.With the remaining gene signatures, logistic regression (LR), naive Bayes (NB) and support vector machine (SVM) algorithms were used to establish PAH and normal tissue classification models by Orange Canvas 3.13 (26).The performance of the three models was assessed by fivefold random sampling and test set validation by calculating the parameters of AUC, precision, specificity and sensitivity.

| Estimation correlation between infiltrating immune cells and hub genes
The normalized gene matrix was analysed by ImmuCellAI online (http://bioin fo.life.hust.edu.cn/ImmuCellAI), which is a gene set signature-based method for estimating the abundance of tissueinfiltrating immune cell populations from transcriptomic data. 21mpared with CIBERSORT mentioned above, the abundance of 18T-cell subsets could be evaluated by ImmuCellAI, such as Th17, Th1 and Treg.The immune cell abundance in all samples was estimated, and Spearman was used to evaluate the correlation between the obtained immune cell content and hub genes in each sample.
The results were visualized by the ggplot2 package within R.

| Animal models of pulmonary arterial hypertension and hemodynamic measurements
Twenty-four male Sprague-Dawley rats (220-250 g) were randomly assigned to the following two groups (n = 12 per group): The PAH group was induced with a single intraperitoneal injection of MCT (60 mg/kg) on Day 0; the parallel group of control rats was injected with an equal volume of saline (Nor).On Day 28, we anaesthetized the rats and measured their pulmonary artery pressure (PAP) through right-sided heart catheterization.After euthanizing the rats, we weighed the wall of the right ventricular (RV) and left ventricle plus the ventricular septum (LV + S), calculated the ratio of RV/LV + S to evaluate the remodelling of the right ventricle, and collected lung tissues for haematoxylin and eosin (HE) staining, immunohistochemistry (IHC) staining and immunofluorescence staining.
Ten male C57BL/6J mice were randomly divided into two groups (n = 5 per group): normoxia group and hypoxia group.These mice were exposed to normoxia (20% O 2 ) or hypoxia (10% O 2 ) for 14 days.The body weight of the mice was recorded every 3 days.
At the end, the mice were anaesthetized (200 mg/kg tribromoethanol) to measure the PAP by using the Millar catheter-transducer (Millar Instruments).Then the mice were euthanized, and the lung tissues were collected for q-PCR, Masson's trichrome staining and IHC staining.

| Haematoxylin and eosin staining, immunohistochemistry, Masson's trichrome staining and multicolor immunofluorescence staining
We followed the methods of Xiang-Guang Wu et al. 22 The left lungs of rat models were embedded in paraffin, sectioned at 5μm thickness sections and then stained with haematoxylin and eosin.For immunohistochemistry staining, briefly, the paraffin-embedded samples were sliced into 8 μm sections.Antigen retrieval was performed using the target retrieval solution Tris-EDTA buffer in a microwave oven for 10-15 min.The sections were placed in 3% H 2 O 2 for 30 min to block the endogenous peroxidase activity.3% BSA was added to evenly cover the tissues to block non-specific binding for 30 min.The tissues were incubated with antibodies against ROCK2 (1:200, 21,645-1-AP, Proteintech) overnight at 4°C, and immune detection was performed using DAB kits (G1212-200T, Servicebio).For Masson's trichrome staining, the paraffin sections were dewaxed and hydrated using a series of solutions, followed by staining with Masson A (overnight), B (1 min), C (1 min), D (6 min), E (1 min) and F (30 s) solutions (G1006, Servicebio).The slices were then rinsed with 1% glacial acetic acid and dehydrated with anhydrous ethanol before being cleared and sealed with neutral gum.Finally, microscope inspection, image acquisition and analysis were performed.
For multicolor immunofluorescence staining, briefly, the paraffinembedded samples were sectioned at 4 mm thickness.Tris-EDTA buffer (pH 8.0) was used to recover antigen with a pressure cooker for 10-15 min.3% BSA was added to evenly cover the tissues to block non-specific binding for 30 min.Subsequently, the sections were incubated with the anti-CD4 Rabbit pAb (1:3000, GB13064-2, Servicebio) in a humidified chamber overnight at 4°C, followed by washing and covering the objective tissue with secondary antibody in dark conditions for 50 min at room temperature.Then the similar experimental methods were used for subsequent staining operations.

| Quantitative real-time polymerase chain reaction
Total RNA from the lungs of the PAH model and controls was isolated using TRIzol, and the extracted RNA was reverse transcribed into complementary DNA (cDNA) using the RT First Strand cDNA Synthesis Kit (Servicebio).Quantitative real-time polymerase chain reaction (q-PCR) was carried out using SYBR Green qPCR Master Mix (Servicebio) to detect the expression of the genes.GAPDH was employed as a reference control.The primer sequences and the thermal protocol were provided in File S1.

| Molecular docking
Molecular docking is a widely used method of exploring the binding modes of small molecules within protein-ligand complexes for structure-based virtual screening.Two docking algorithms, LibDock and CDOCKER of DISCOVERY STUDIO 2018 performed in molecular docking analysis.LibDock, a site-feature docking algorithm developed by Diller and Merz, can easily provide a rapid virtual screening campaign on millions of compounds by matching the physicochemical properties of the compounds to the polar and polar features in the protein binding sites. 23CDOCKER is another powerful docking method that is CHARMm-based and has been validated for generating highly accurate docked poses. 24The crystal structure of ROCK2

| Statistical analysis
We used R software and GraphPad Prism 8 to perform statistical analyses, presented data as mean ± standard error (SD) and calculated statistical significance by the t-test.p < 0.05 were considered statistically significant.

| Screening of DEGs
A flowchart of the analysis procedures, including transcriptomic analysis and experimental validation, for this study was shown in

| GO and KEGG pathway enrichment analysis
To identify the biological roles and signalling pathways related to the 475 DEGs, GO and KEGG enrichment analyses were conducted by the Cluster profiler R package.According to the GO enrichment analysis, most DEGs were significantly enriched in protein kinase binding, extracellular exosomes, positive regulation of GTPase activity, etc. (Figure 4A).KEGG pathway analysis revealed that most DEGs were primarily enriched in the PI3K-Akt signalling pathway, lipid and atherosclerosis, protein processing in the endoplasmic reticulum, Th17 cell differentiation and the NOD-like receptor signalling pathway (Figure 4B).

| Estimating the characteristics of immune cells
Based on CIBERSORT analysis, the abundance of immune cell types in each sample was selected, and significant results (p < 0.05) from the three GEPs were presented in Figure 5A.The differences between the 22 types of immune cells and individual samples were shown in the heatmap with hierarchical clustering (Figure 5B).
The analysis presented that naïve CD4 T cells were downregulated, whereas resting memory CD4 T cells and activated memory CD4 T cells were upregulated in the lung tissues of PAH patients.

| In silico analysis of immune-related genes
The PPI network of processed DEGs determined with significant results (p < 0.05) was constructed using the STRING database and Cytoscape software.After removing the unconnected nodes, a PPI network that contained 58 genes was constructed (Figure 6A,B).CD52, HSP90AA1 and ACTR2.The six genes were used as immunerelated genes to build NB, LR and SVM classification models.Table 2 shows the performance of the four models as estimated using fivefold random sampling and test set validation.The NB and LR classification models showed better performance, and both the precise and AUC of the three evaluation methods were beyond 0.918.ROC analysis also showed that the NB and LR algorithms with the six-gene prognostic signature had high sensitivity and specificity for PAH classification.

| Estimation correlation between infiltrating immune cells and hub genes
According to CIBERSORT analysis, the differential cell type of the infiltrating immune cell population in PAH and controls was CD4 T cells.To estimate the correlation between various subtypes of CD4 T cells and the obtained genes, we analysed the normalized gene matrix by ImmuCellAI and Spearman.ROCK2, CD52, HSP90AA1 and ACTR2 were positively correlated with iTreg, nTreg, Tfh, Th1, Th17 and Th2, while ATHL1 and IL-2rb were negatively correlated with naive CD4 T cell, iTreg, Tfh, Th1, Th17, Th2 and Tr1 (Figure 6C).

| Validation of the density and distribution of infiltrating immune cells and immune-related genes in MCT-induced PAH through q-PCR
MCT-induced PAH in rats is a well-studied and published method to make preclinical models of PAH.However, there are currently no perfect animal models that replicate entirely human PAH.To investigate whether the pathological changes of PAH in human also occur in MCT-induced rats and to validate the expressions of hub genes found in our analysis, we established a PAH model.After

| Validation of the immune-related genes in hypoxia-induced PAH mice model
Following a 14-day exposure to hypoxia, increased medial thickness, collagen deposition and muscularization of the distal arterioles were observed in hypoxia-induced PAH mice model (Figure 8A), which is consistent with the model characteristics.And the increased expression of ROCK2 in the thickened and muscled microarteries of PAH mice was observed using IHC (Figure 7D,E).As shown in Figure 8C, the expressions of ROCK2, ATHL1, HSP90AA1 and ACTR2 were significantly upregulated, which was consistent with the results in the datasets of PAH patients and MCT-induced PAH models.Additionally, the expression of CD52 was also increased, while IL-2rb was not detected.These differences may be related to species differences and sample sizes.

| Molecular docking
Among the hub genes, HSP90AA1 and ROCK2 were most significantly upregulated in PAH tissues.It was reported that ROCK2 kinase inhibitors can reduce the phosphorylation of STAT3 and enhance the phosphorylation of STAT5, thereby downregulating overactivated T cells and rebuilding immune balance. 25In addition, Then we analysed the binding posture of each molecule to the docking pocket, focusing on the interaction with the aspartic acid (ASP176, ASP218, etc.) on the side of the pocket and whether the hydrophobic group can form a hydrophobic interaction with the hydrophobic region.Finally, 20 molecules were obtained by removing those poorly bounded molecules (Table 3).Figure 9. shows three randomly selected molecules docking results according to the classification of molecular structure.
TA B L E 3 Docking analysis of virtual screening hits (key hydrogen bond interactions and CDOCKER _INTERACTION_ENERGY).  Th17 differentiation via multiple mechanisms. 280][31] As a new therapeutic idea, immune microenvironment and immunotherapy have been widely studied in various diseases. 32It is of great significance to understand the progression of PAH disease from this perspective.Using the PPI network and machine learning, six hub immune-related genes were figured out.And then we demonstrated that the accumulation of CD4 T cells was infiltrated in and around vascular lesions, and the expressions of ROCK2, ATHL1, HSP90AA1 and ACTR2 were significantly upregulated in both MCT-induced PAH models and hypoxia-induced PAH mice.Among those with these validated genes, ROCK2, HSP90AA1 and ACTR2 were positively correlated with various CD4 T cells including Tregs, Tfh, Th1, Th17 and Th2, while ATHL1 were negatively correlated with them.At present, we know little about ATHL1.Although we know that PGGHG was the product of ATHL1 gene, 33 there were few studies on the function of ATHL1, which would be a novel research target.
HSP90AA1 promotes autophagy through PI3K/Akt/mTOR pathway and inhibits apoptosis through JNK/P38 pathway.Autophagic abnormalities are observed in PAECs, PASMCs and in the RV of PH from both animal models and patients. 34,35Although there had no direct evidence that HSP90AA1 plays a role in the process of PAH, HSP90AA1 is known to interact with endothelial nitric oxide synthase (eNOS), which is responsible for nitric oxide, a potent vasodilator. 36By stabilizing and promoting the activity of eNOS, HSP90AA1 may help maintain proper vascular tone and prevent excessive vasoconstriction, a key feature of PAH.8][39][40] Rho-ROCK2 2][43] Moreover, ROCK2 signalling pathway regulates the Th17/Tregs balance and controls profibrotic pathways. 44Toru The results of the molecular docking between TDI01953 and ROCK2, pemetrexed acid and ROCK2, radotinib and ROCK2.
The pyrazole ring in the drug skeleton structure forms a stable hydrogen bond with the hinge area MET172, the amino hydrogen of the indoline ring forms a hydrogen bond with ASP218, and which also forms a van der Waals interaction with ASP176.Therefore, we took ROCK2 as a potential therapeutic target to find a list of 20 small-molecule compounds that might modulate its activity through molecular docking simulations, including TDI01953, pemetrexed acid and radotinib.Pemetrexed is a folate antagonist used in the treatment of non-small cell lung cancer and malignant mesothelioma. 47Radotinib, a high-affinity BCR-ABL1 inhibitor mainly used for resistance to imatinib, was approved in Korea for second-line CML treatment in 2012. 48Notably, TDI01953 is an analogue of TDI01, which is a specific selective ROCK2 inhibitor undergoing clinical trials for pulmonary fibrosis and nonalcoholic steatohepatitis currently.Even though there is no listed ROCK2 inhibitor for the treatment of PAH, it has been evaluated in several clinical trials currently for various indications including, autoimmune diseases and cardiovascular disorders.Some of the ROCK2 inhibitors that have been tested in clinical trials include KD025 and Fasudil.KD025, a ROCK2specific inhibitor that reduces the expression of IL-17 and fibrogenesis by inhibiting the ROCK2 signalling pathway. 44,49A phase IIa clinical trial for the treatment of chronic graft-versus-host disease (cGVHD) showed quality of life improvements, corticosteroid dose reductions and limited toxicity. 44Additionally, it is worth mentioning that it can inhibit the increase of right ventricular systolic blood pressure in MCT-induced rats. 50The nonselective ROCK inhibitor Fasudil, approved for the treatment of cerebral vasospasm, can reduce the expression levels of endothelin-1 and endothelin receptor, inhibit the proliferation of endothelial cells in pulmonary arterioles, and promote the generation of nitric oxide in rats with pulmonary hypertension caused by left heart disease. 51ROCK2 inhibitors show promise as a novel therapeutic approach for PAH.However, further clinical trials are needed to evaluate their safety and efficacy in larger patient populations.
Taken together, our study identified ROCK2, ATHL1, HSP90AA1 and ACTR2 as immune-associated genes and provided insights into the reverse of the pulmonary vascular remodelling from the perspective of immunity and inflammation, which may be an effective and promising therapeutic strategy for PAH.There are some limitations to this study.First, this study employed microarray analysis, and all findings were predicated on gene expression levels.However, it should be noted that disparities at the genetic level of biomarkers do not inevitably translate to differences in protein-level function and mechanism.The detailed molecular mechanisms associated with these four genes should be further addressed elaborately in the future.Second, it is still difficult to discern whether inflammation is the cause or initiating factor of pulmonary vascular remodelling or the consequence of the pulmonary high-pressure microenvironment.
Third, studies showed that pathologic specimens from PAH patients reveal an accumulation of inflammatory cells in and around vascular lesions, including macrophages, T and B cells, dendritic cells and mast cells. 52However, in this study, we only obtained the result of the differences in T cells, which may be due to the bias of the limited samples included in our data set and the cellular heterogeneity in lung tissues.Furthermore, this study provides preliminary evidence of a correlation between inflammation and immune and PAH patho- was acquired from the Protein Data Bank (PDB ID: 6ED6).The interaction energy for every final pose of ligands with CHARMm energy was computed, and the top scoring poses are retained.The structure of ROCK2 was prepared by removing water and adding hydrogen to a clean protein module.The compounds were also prepared by adding hydrogen conversing into 3D structures and pH-based ionization.We used molecular docking DS 2016 to identify potential ROCK2 inhibitors by screening the Clinical & FDA-approved Drug Library (https://www.apexbio.cn/disco veryp robet m-clini calfda-appro ved-drug-libra ry.html).The crystal structure of ROCK2, complexed with N-[(2,3-dihydro-1,4-benzodioxin-5-yl) methyl]-4-(pyridin-4-yl) benzamide, was retrieved from the Protein Data Bank (PDB ID: 6ED6).Protein preparation was conducted by removing the water and co-crystallized N-[(2,3-dihydro-1,4-benzodioxin-5-yl) methyl]-4-(pyridin-4-yl) benzamide, adding hydrogens and missing loop regions, calculating protein ionization and protonating the protein structure.The binding site was defined by the primary ligand.After hydrogenating the eutectic ligand and minimizing the energy, CDOCKER redocks with the ROCK2 protein (6ED6) again.The top 10 positions with the lowest energy reserved were sorted in ascending order (CDOCKER_INTERACTION_ENERGY).The eutectic ligand in 6ED6 was copied into the docking results, and the RMSD values of the non-hydrogen atoms in 10 positions were calculated with it as a reference, as shown in

1
The RMSD values of redocked conformations compared with coligand.The redocked conformations are arranged in ascending order of the CDOCKER_INTERACTION_ENERGY values.topfour lowest energy positions were all less than 0.5.The eutectic ligands were superposed in the four lowest energy positions, as shown in Figure1.Eutectic ligands almost completely overlapped with the four positions, which means the CDOCKER method works well.The eutectic ligand structure was copied into the database, and 1969 matching molecules without bases were selected from the clinical approval trial database.After hydrogenation, the CHARMm Force field is used for energy minimization and Momanye Rone electrostatic charge is loaded on the molecules.

Figure 2 .
Figure 2.There were 84 PAH samples and 49 normal controls from the GSE53408, GSE117261 and GSE15197 datasets merged and processed using the sva package in R, and 2758 genes were analysed and normalized in total (Figure 3A,B).Differential expression analysis was then performed using the limma Bioconductor package, with a threshold of |log 2 fold change (FC)| ≥ 1 and p ≤ 0.05.This analysis identified 221 upregulated genes and 254 downregulated genes in PAH tissues when compared to normal samples, with a total of 475 DEGs identified, as shown in Figure 3C.
The correlation analysis of PAH and controls further revealed the clues to the correlation among 22 immunocytes (Figure5D,E).The results indicated that the correlation of immune cell subsets, including NK cells, resting dendritic cells, CD8 T cells and resting memory CD4 T cells, was strong in the normal group, while the correlation with PAH was weak.However, the correlation of CD4 T cells, including resting memory CD4 T cells, activated memory CD4 T cells and naïve CD4 T cells, remained significant compared with other cell subsets.Notably, the correlation of activated memory F I G U R E 1 Molecular superposition diagram (yellow ball-and-stick model for eutectic ligands, grey for the lowest energy values of the four conformations).| 3869 HE et al.CD4 T cells with CD8 T cells and memory B cells was amplified in PAH.

F I G U R E 3 F I G U R E 4 F I G U R E 5
Thirty-one genes were considered hub genes, with an criteria of edge count beyond three.These genes were HSP90AA1, HSP90AB1, NCL, RANBP2, ASH1L, RPS3A, CEVPE, DNAJC3, UQCRC2, CCT8, HELLS, PRKC1, SMG1, KIF20B, ANLN, CDCA2, CHUK, CYBB, CYTIP, NAP1L1, SOD2, ACTR2, CD52, DIEXF, DNAJC21, ERC1, IL2RB, MACF1, MIB1, ROCK2 and XRN1, which may play an important role in PAH progression.To further explore the immune-related gene signatures modulated during PAH development, PAH classifiers were built using machine learning algorithms.Before developing the models, Pearson F I G U R E 2 Study workflow.Screening of differentially expressed genes (DEGs) in pulmonary arterial hypertension (PAH) tissues compared with control tissues.(A, B) GSE53408, GSE117261 and GSE15197 datasets (49 normal controls and 84 PAH samples) were merged and normalized by sva package.(C) Volcano plots of all quantified genes in the gene expression profiles analysis of PAH and controls.Volcano map screened out statistically significant DEGs (p < 0.05 and |log 2 FC| > 1).Among these DEGs, there were 221 upregulated genes (represented by red dots) and 254 downregulated genes (represented by green dots).correlation analysis was applied to filter the 31 hub genes.There were 20 genes with p < 0.05, which were used as the input variables for stepwise regression.Stepwise regression basically performs multiple regressions, removing the weakest relevant variables each time.After stepwise regression, six variables were left, which can best illustrate the distribution, including IL-2rb, ROCK2, ATHL1, Functional enrichment analysis.(A) Enriched Gene Ontology (GO) terms for molecular function (MF), biological process (BP) and cellular component (CC).(B) Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the pathways for differentially expressed genes (DEGs).The x axis shows the gene numbers and the colour of bar represents the p-value.There were 16 pathways meeting the cutoff value of p < 0.05.Analysis of immune landscape in the three datasets.(A) Histogram showing the abundance of immune cell subtypes in each sample (B) Heatmap of differentially cell types and individual samples.The expression profiles greater than the mean are coloured in red, and those below the mean are coloured in green.Blue, normal lung tissues; pink, PAH specimens.(C) violin plot showing the comparisons between immune cells in normal controls and PAH patients in the three datasets.Correlation matrix of the immunocyte proportions in Normal samples (D) and PAH samples (E).The colour of red and blue represent positive and negative correlations, respectively.The value represents the correlation coefficient.
a 28-day exposure to MCT, rats had a significant elevation in PAP (Pulmonary Artery Pressure) and RV/LV + S compared with those in control rats (Figure 7A).Meanwhile, MCT treatment resulted in pulmonary vascular remodelling, as evidenced by the increased media (Figure 7B).As shown in Figure 7C, the expressions of ROCK2, ATHL1, HSP90AA1 and ACTR2 were significantly upregulated, which was consistent with the results in the datasets of PAH patients.Meanwhile, we found that the expression of CD52 and IL-2rb showed no difference between PAH models and controls.Then, we validated the density and distribution of infiltrating CD4 T cells by multicolor immunofluorescence staining and detected the expression of a representative hub gene (ROCK2) in successfully established PAH models by immunohistochemistry (Figure 7D,E).

F I G U R E 6
Immune-related genes interaction and selection.(A) Immune-related genes interaction network built by using STRING database and Cytoscape software.(B) The degree of the immune-related genes with count of node more than three.(C) the correlation between various subtypes of immune cells and obtained genes by ImmuCellAI and Spearman.TA B L E 2 Performance of logistic regression, support vector machine (SVM) and naive Bayes models estimated by 10-fold, random sampling and test validation.

F I G U R E 7
Validation by MCT-induced PAH rat model.(A) Representative tracing of PAP in rat after MCT or Nor treatment and the mean values of PAP and Fulton index (RV/LV + S) in the two groups.****p < 0.0001 versus Nor group.(B) Representative images of haematoxylin and eosin staining of lung sections from PAH rat or Nor rat.(C) The expression of ROCK2, ATHL1, HSP90AA1, ACTR2, IL-2rb and CD52.n = 12 in Nor group, n = 8 in PAH group (MCT-induced mortality in the model was about 33%).*p < 0.05, **p < 0.01, ns (no significance) versus Nor group.(D) IHC staining for ROCK2 from normal and PAH rats.(E) the density and distribution of infiltrating CD4 T cells (CD3+ CD4+)by multicolor immunofluorescence staining and detected the expression of representative hub gene (ROCK2) in successfully established PAH models by immunohistochemistry. F I G U R E 8 Validation through hypoxia-induced PAH mice model.(A) Representative images of IHC of ROCK2 from PAH mice or Nor mice.(B) Representative images of Masson's trichrome staining of lung sections from PAH mice or Nor mice.(C) The expression of ROCK2, ATHL1, HSP90AA1, ACTR2 and CD52.n = 5 in each group.*p < 0.05, **p < 0.01 versus Nor group.| 3873 HE et al.
inhibitors can also prevent tissue fibrosis and induce the regression of established fibrosis, which would relieve vasoconstriction.Here, we used ROCK2 as a potential therapeutic target and tried to screen targeted small-molecule compounds through molecular docking.We used the binding site sphere (x = 26.896090,y = 47.022057,z = 52.444471,r = 9.798772) of 6ED6 for molecular docking by LibDock.The docked pose with the highest LibDock score was obtained for each compound.Each molecule retained 10 conformations with the lowest energy.A total of 11,860 conformations were retained from the docking results.Then we analyse the hydrogen bond, hydrophobic, electrostatic and other interactions between molecules and amino acid residues by using the analysed ligand pose in the receptor-ligand interaction module.Based on the result of analysed ligand poses, 638 conformations were retained by removing all conformations that formed hinge regions.Under the condition that each molecule retaining the lowest energy conformation and CDOCKER_INTERACTION_ENERGY value less than −40 kJ/mol, 79 molecules were screened by using the filter tool.

4 |
D ISCUSS I ON AND CON CLUS I ON SIn our study, a total of 475 DEGs (221 upregulated and 254 downregulated) were identified in the transcriptional profiles in PAH tissues, mainly involved in pathways related to inflammatory responses or response to hypoxia.Using CIBERSORT analysis, we found that CD4 T cells were the differential cell type of the infiltrating immune cell population in PAH.Our findings suggest that the PI3K-Akt signalling pathway was enriched with most differentially expressed genes, which is important in the development of PAH associated with CD4 T cells, including anti-inflammatory Tregs and proinflammatory Th17 cells.The differentiation of Th cells from naive CD4 T cells is initiated by the stimulation of T-cell receptors (TCRs) and associated costimulatory molecules upon antigen presentation, which leads to the activation of PI3K and mTORC.26,27PI3K-Akt-mTORC1 signalling positively regulates The hydrophobic ring of methyl indole and form hydrophobic interactions with ALA231, GLY101, LEU123, ALA102, PHE103, etc. (A1) The surfaces (H-Bonds) of ROCK2 (6ED6) binding with compound TDI01953.(A2) Binding maps of TDI01953 with the active site of ROCK2.(A3) TDI01953 (Stick) interactions with amino residues (Line) of the active site of ROCK2.The amine group on the pyrimidine ring forms two hydrogen bonds with the hinge region MET172, which can also form hydrogen bonds with ALA102, PHE103 and ASP232.Two carboxyl groups can form a strong electrostatic interaction with the loop ring and two lysines (LYS121 and LYS216) in the catalytic region.(B1) The surfaces (H-Bonds) of ROCK2 (6ED6) binding with pemetrexed acid.(B2) Binding maps of pemetrexed acid with the active site of ROCK2.(B3) pemetrexed acid (Stick) interactions with amino residues (Line) of the active site of ROCK2.The trifluoromethyl group of radotinib forms a hydrogen bond with MET172 of ROCK2, and which can also form a hydrogen bond with ASP176, ARG100 and LYS121; the pyrazine ring can form a π-π stack with ASP232.(C1) The surfaces (H-Bonds) of ROCK2 (6ED6) binding with radotinib.(C2) Binding maps of radotinib with the active site of ROCK2.(C3) Radotinib (Stick) interactions with amino residues (Line) of the active site of ROCK2.
genesis; experiments need to be designed to further elucidate their mechanisms of action; and further in vivo data and clinical trials are warranted to investigate the role of inflammation and immune in PAH progression.

Table 1 .
The RMSD values of the