Screening and identification of key regulatory connections and immune cell infiltration characteristics for lung transplant rejection using mucosal biopsies.

This study aimed to explore key regulatory connections underlying lung transplant rejection. The differentially expressed genes (DEGs) between rejection and stable lung transplantation (LTx) samples were screened using R package limma, followed by functional enrichment analysis and protein-protein interaction network construction. Subsequently, a global triple network, including miRNAs, mRNAs, and transcription factors (TFs), was constructed. Furthermore, immune cell infiltration characteristics were analyzed to investigate the molecular immunology of lung transplant rejection. Finally, potential drug-target interactions were generated. In brief, 739 DEGs were found between rejection and stable LTx samples. PTPRC, IL-6, ITGAM, CD86, TLR8, TYROBP, CXCL10, ITGB2, and CCR5 were defined as hub genes. Eight TFs, including STAT1, SPIB, NFKB1, SPI1, STAT5A, RUNX1, VENTX, and BATF, and five miRNAs, including miR-335-5p, miR-26b-5p, miR-124-3p, miR-1-3p, and miR-155-5p, were involved in regulating hub genes. The immune cell infiltration analysis revealed higher proportions of activated memory CD4 T cells, follicular helper T cells, γδ T cells, monocytes, M1 and M2 macrophages, and eosinophils in rejection samples, besides lower proportions of resting memory CD4 T cells, regulatory T cells, activated NK cells, M0 macrophages, and resting mast cells. This study provided a comprehensive perspective of the molecular co-regulatory network underlying lung transplant rejection.

Introduction

Lung transplantation (LTx) has been considered as the ultimate rescue therapy for patients with end-stage lung diseases, such as pulmonary hypertension [1], pulmonary fibrosis [2], and recent COVID-19-related respiratory failure [3]. In the recent era, nearly 4000 LTxs are performed annually worldwide [4]. Nevertheless, the median 5-year survival rate and the median 10-year survival rate following LTx, are 54% and 32%, respectively, which are worse compared with those for other solid organ transplantations and have not changed substantially over the past two decades [5], [6]. The leading causes of death following LTx are primary graft dysfunction (PGD) and chronic lung allograft dysfunction (CLAD), contributing to early and late complications, respectively [7].

Rejection is a major complication during the development of PGD and CLAD, mainly including acute cellular rejection, antibody-mediated rejection, and chronic rejection [8]. Rejection is closely related to bronchiolitis obliterans syndrome, which represents a persistent obstructive decline in lung function and is regarded as an endpoint in most clinical studies [9]. The gold standard to detect rejection following LTx is the histopathological grading of transbronchial biopsies; however, it is associated with a considerable risk of complications and limited sensitivity, specificity, and reproducibility [10], [11], [12]. Thus, potential molecular biomarkers for lung transplant rejection need to be identified using other detection methods.

Recently, a novel and safe method was developed for diagnosing rejection in lung transplants through the molecular assessment of lung mucosal biopsies [13]. Stable and rejection LTx mucosal samples were diagnosed. However, the differences in gene expression and the immune status between these samples were not assessed. This study was novel in using bioinformatics analysis to screen genetic alterations and identify hub genes playing a critical role in lung transplant rejection [14]. The protein–protein interaction (PPI) network, an integrative miRNA–transcription factor (TF)–mRNA co-regulatory network, and drug–target interactions in lung transplant rejection were also generated. Furthermore, immune cell infiltration characteristics in stable and rejection LTx samples and the relationships between hub genes and immune cells were also revealed, further clarifying the key mechanisms in molecular immunology underlying lung transplant rejection.

Materials and methods

Microarray data

Relevant microarray data regarding lung transplant rejection was acquired from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). The mRNA microarray dataset GSE125004 deposited by Chang et al. was downloaded from the GEO database. The database comprised mucosal biopsies from the third bronchial bifurcation (3BMBs) of lung transplant recipients, including 167 stable and 24 rejection samples; the mean post-transplant biopsy time in the stable and rejection groups was 578 and 629 days, respectively (see details in ClinicalTrials.gov NCT02812290) [13]. The GPL16043 GeneChip PrimeView Human Gene Expression Array (with external spike-in RNAs) was used to obtain gene expression profiles. The flow chart of the procedures used in the present study is shown in Fig. 1 .

Fig. 1

Flow chart of the study procedures.

Screening of differentially expressed genes

The R 3.5.0 software (https://www.r-project.org/) and the limma package (version 3.30.11, http://www.bioconductor.org/packages/release/bioc/html/limma.html) from the Bioconductor project were used to screen the differentially expressed genes (DEGs) between rejection and stable lung transplant samples [15], [16]. First, data were normalized using the normalizeBetweenArray function from R package limma (Fig. 2 ). Subsequently, the t test and the Benjamini–Hochberg method were used to calculate the P value and adjusted P value (FDR), respectively. Finally, DEGs were identified under the cutoff thresholds: FDR < 0.05 and |log2FC|  > 1.

Fig. 2

Boxplots of sample data before and after normalization.

Functional enrichment analysis

Gene Ontology (GO) offers a biological model that classifies gene functions into three categories: cellular components (CC), biological processes (BP), and molecular functions. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database that can identify functional and metabolic pathways using genome sequences or high-throughput data. The “Custom Analysis” mode of Metascape was used to perform GO and KEGG pathway analyses [17], with a P value < 0.01 as the cutoff criterion.

Protein–protein interaction (PPI) network construction, module screening, and hub gene identification

The Search Tool for the Retrieval of Interacting Genes (STRING 11.0; http://string.embl.de/) [18] is a biological database and web resource that predicts comprehensive interactions of genes at the protein level. The parameter was set as medium confidence > 0.4, and the PPI network of DEGs was screened. Subsequently, the PPI network was visualized using Cytoscape software 3.6.1 [19]. The significant modules of the PPI network were selected using Molecular Complex Detection (MCODE) plug-in [20], and node score cutoff, 0.2; K-Core, 2; max depth, 100; degree cutoff, 4; and MCODE score > 10 were set as the cutoff criterion. In addition, nodes with a high degree of connectivity contribute more to the stability of the PPI network, and hence DEGs with degree connectivity of >130 were defined as hub genes using the NetworkAnalyzer [21].

miRNA–TF–DEG regulatory network analysis

The miRNA–DEG regulatory network was predicted and visualized using miRNet ( ) [22], [23], a comprehensive analytical tool that integrates multiple high-quality miRNA-target data sources from 11 databases (miRecords, miRanda, PharmacomiR, PhenomiR, miRTarBase, starBase, miR2Disease, SM2miR, TarBase, HMDD, and EpimiR). The cutoff criterion was set as follows: organism, Homo sapiens; tissue, lungs; degree cutoff for miRNA nodes, 20.

The TF–DEG regulatory network was identified using the iRegulon plug-in in Cytoscape, which included the TF-target pairs of multiple human databases such as Transfac and Encode [24]. The TF motifs with a normalized enrichment score (NES) > 3 and targets > 150 were considered as the threshold for selecting potential relationships. Ultimately, the miRNA–TF–DEG regulatory network was constructed using Cytoscape.

Immune cell infiltration evaluation and analysis

CIBERSORT [25] is an analytical tool developed by Newman to provide an estimate of the relative levels of 22 phenotypes of human hematopoietic cells in a mixed cell population using gene expression data. Normalized gene expression profiles were uploaded to the CIBERSORT web portal (), and the algorithm was run using the LM22 signature and 500 permutations. Cases that met the CIBERSORT P < 0.05 requirements, indicating that the inferred fractions of immune cell populations produced by CIBERSORT were accurate, were considered to be eligible for further analysis. In each sample, the proportion of all immunocyte types was equal to 1.

Prediction of drug–target interactions

Drugs were selected based on the hub genes that served as promising targets using the Drug–Gene Interaction Database (DGIdb; ) [26]. In this study, the final drug list included only drugs approved by the Food and Drug Administration (FDA). Drug–gene interactions were constructed and visualized using Cytoscape.

Statistical analysis

A two-tailed Student t test was used to analyze the differences between immune cell fractions of eligible rejection and stable lung transplant samples using GraphPad Prism 7.0 software. In addition, Pearson correlation analysis was used to explore the relationship between immune cell proportions in rejection and stable lung transplant samples and the relationship between the expression of hub genes and immune cell proportions in rejection lung transplant samples.

Results

Identification of DEGs

Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.intimp.2020.106827.

A total of 18,835 genes were detected in lung mucosal biopsies from 191 lung transplant recipients, of which 739 [459 (62.11%) upregulated and 280 (37.89%) downregulated] were identified as DEGs (Table S1 ). The most significantly upregulated and downregulated genes were indoleamine 2,3-dioxygenase 1 (logFC = 4.22) and glutathione S-transferase alpha 2, respectively (logFC = −2.65). The volcano plot of DEGs is shown in Fig. 3 A, and the expression levels of top 50 upregulated and top 50 downregulated DEGs are represented as a heat map in Fig. 3B.

Table S1

Fig. 3

(A) Volcano plot of all DEGs. Top five upregulated and top five downregulated DEGs are marked. The X-axis represents log2 FC, and the Y-axis represents the log-transformed adjusted P values. (B) Heatmap of the top 100 DEGs. The X-axis represents samples, and the Y-axis represents DEGs. IDO1, Indoleamine 2,3-dioxygenase 1; CXCL11, C-X-C motif chemokine ligand 11; CXCL10, C-X-C motif chemokine ligand 10; CXCL9, C-X-C motif chemokine ligand 9; CCL18, C-C motif chemokine ligand 18; GSTA2, glutathione S-transferase alpha 2; GRP, gastrin-releasing peptide; CYP26A1, cytochrome P450 family 26 subfamily A member 1; SULT1E1, sulfotransferase family 1E member 1; SLC23A1, solute carrier family 23 member 1.

Enrichment analysis of DEGs

All DEGs were subjected to the GO and KEGG pathway enrichment analyses, revealing several functional categories and pathways. The 459 upregulated DEGs were enriched mainly in the categories “cytokine-mediated signaling pathway” (GO:0019221), “lymphocyte activation” (GO:0046649), “activation of immune response” (GO:0002253), “adaptive immune response” (GO:0002250), and “leukocyte migration” (GO:0050900) (Fig. 4 A, B). The 280 downregulated DEGs were enriched mainly in the categories “cilium” (GO:0005929) and “motile cilium” (GO:0031514) (Fig. 4C, D).

Fig. 4

Functional enrichment analysis of DEGs using Metascape. (A and C) Top 20 enriched terms of up- and downregulated DEGs. The color of terms from light to dark represents the P value of terms from high to low. (B and D) Network of the top 20 enriched terms. Each term is represented by a circle node, where its size is proportional to the number of DEGs falling into that term. Nodes of the same color belong to the same cluster, and the thickness of the edge represents the similarity score.

PPI network analysis and hub gene identification

The PPI network of DEGs was constructed and visualized using Cytoscape based on the STRING database. In total, 864 nodes (genes) and 3878 edges (interactions) were screened out (Fig. S1 ). After using the MCODE plug-in to identify modules from the PPI network, the top three central modules with MCODE scores > 10 were selected. Module 1 with MCODE scores 40.857 consisted of 43 nodes and 858 edges (Fig. 5 A). Module 2 with MCODE scores 34.124 consisted of 98 nodes and 1655 edges (Fig. 5B). Module 3 with MCODE scores 15.623 consisted of 54 nodes and 414 edges (Fig. 5C). Through the GO and KEGG pathway enrichment analyses, Module 1 was found to be enriched mainly in “Chromosome segregation” (GO: 0007059) and “Mitotic cell cycle phase transition” (GO: 0044772). Module 2 was enriched mainly in “Cytokine-mediated signaling pathway” (GO: 0019221) and “Response to virus” (GO: 0009615), while Module 3 was enriched mainly in “Staphylococcus aureus infection” (hsa03010) and “Side of membrane” (GO: 0098552) (Table 1 ).

Fig. S1

PPI network analysis of DEGs. The PPI network of DEGs was constructed and visualized using Cytoscape based on the STRING database. Upregulated DEGs and downregulated DEGs are marked in red and green, respectively. The size of nodes from small to large represents the degree of connectivity of nodes from low to high.

Fig. 5

Module analysis of the PPI network. Three central modules of the PPI network were identified and visualized using the MCODE plug-in in Cytoscape. (A) Module 1. (B) Module 2. (C) Module 3. The size of nodes from small to large represents the degree of connectivity of nodes from low to high, and the hub genes are marked in yellow.

Table 1

Top enriched GO-BP and KEGG pathways of DEGs in the top three modules.

Cluster	Term	Description	Count	LogP
Module 1
GO: BP	GO: 0007059	Chromosome segregation	26	–37.271
GO: BP	GO: 0044772	Mitotic cell cycle phase transition	20	–20.554
GO: CC	GO: 0000775	Centromeric region of the chromosome	15	–20.538
Module 2
GO: BP	GO: 0019221	Cytokine-mediated signaling pathway	58	–59.697
GO: BPGO: BPKEGGKEGGKEGG	GO: 0009615GO: 0034341hsa05323hsa05164hsa05332	Response to the virusResponse to interferon-gammaRheumatoid arthritisInfluenza AGraft-versus-host disease	342513158	–37.898–29.700–16.347–15.404–11.329
Module 3
GO: CC	GO: 0098552	Side of the membrane	21	–19.483
GO: BPGO: BP	GO: 0002250GO: 0001817	Adaptive immune responseRegulation of cytokine production	2017	–16.817–12.460
KEGGKEGGKEGG	hsa03010hsa04380hsa04650	Staphylococcus aureus infectionOsteoclast differentiationNatural killer cell–mediated cytotoxicity	1385	–22.892–9.320–4.923

Seven DEGs with a degree of connectivity > 130 were defined as hub genes for lung transplant rejection. These DEGs included protein tyrosine phosphatase receptor type C (PTPRC; degree = 162), interleukin 6 (IL-6; degree = 160), integrin subunit alpha M (ITGAM; degree = 158), CD68 (degree = 149), Toll-like receptor 8 (TLR8; degree = 132), TYRO protein tyrosine kinase–binding protein (TYROBP; degree = 129), chemokine (C-X-C motif) ligand 10 (CXCL10; degree = 125), integrin, beta 2 (ITGB2; degree = 125), and chemokine (C-C motif) receptor 5 (CCR5; degree = 123) (Table 2 ). Of note, PTPRC, IL-6, CD86, TYROBP, and ITGB2 were included in Module 2, and ITGAM, TLR8, and CCR5 were in Module 3.

Table 2

Top nine hub genes with a higher degree of connectivity in the PPI network.

Gene	Degree	Stress	BetweennessCentrality	Eccentricity	ClusteringCoefficient
PTPRC	162	348,336	0.037717	6	0.268231
IL-6	160	542,604	0.072158	5	0.219811
ITGAM	158	284,234	0.028023	6	0.279529
CD86	149	221,138	0.017972	6	0.30156
TLR8	132	196,016	0.017242	6	0.312977
TYROBP	129	223,676	0.022391	6	0.276405
CXCL10	125	164,278	0.0143	6	0.330323
ITGB2	125	204,386	0.018708	6	0.312516
CCR5	123	136,128	0.009614	6	0.360123

MiRNA–TF–DEG regulatory network analyses

The miRNA–DEG regulatory network was predicted and constructed using miRNet to further understand the regulatory relationship between DEGs and their upstream miRNAs. As shown in Fig. 6 , the miRNA–DEG network consisted of 23 miRNAs and 385 DEGs. Of note, four hub genes, including IL-6, CD86, CXCL10, and ITGB2, were involved in the miRNA–DEG network. Top 7 miRNAs having more than 30 targets in the network, which may play critical roles in regulating DEGs, are shown in Table 3 . Among these, five miRNAs, including miR-335-5p, miR-26b-5p, miR-124-3p, miR-1-3p, and miR-155-5p, were predicted to target hub genes.

Fig. 6

MiRNA–DEG regulatory network analysis. MiRNAs that target DEGs were predicted using miRNet. The network was visualized using miRNet. The rectangles represent miRNAs, and the circles represent DEGs. Hub miRNAs and DEGs are marked in purple and yellow, respectively.

Table 3

Top seven hub miRNAs with a higher degree of connectivity in the miRNA–DEG regulatory network.

miRNA	Degree	Betweenness
hsa-mir-335-5p	133	36,761.29
hsa-mir-26b-5p	77	18,629.87
hsa-mir-124-3p	44	9199.693
hsa-mir-192-5p	37	5576.077
hsa-mir-1-3p	35	6479.095
hsa-mir-92a-3p	32	5412.899
hsa-mir-155-5p	31	5479.08

The iRegulon analysis revealed eight most important TFs that regulated DEGs, including STAT1, SPIB, NFKB1, SPI1, STAT5A, RUNX1, VENTX, and BATF. The TF–DEG regulation network of 611 nodes and 2306 interaction pairs was constructed using Cytoscape to determine the regulatory connections between TFs and DEGs (Fig. 7 ; Table 4 ). Finally, the miRNA–TF–DEG regulatory network was established with eight TFs, five miRNAs, and nine hub genes (Fig. 8 ).

Fig. 7

TF–DEG regulatory network analysis. TFs that regulate DEGs were identified using the iRegulon plug-in in Cytoscape. The network was visualized using Cytoscape. Brown hexagons represent TFs, red circles represent upregulated DEGs, and green circles represent downregulated DEGs.

Table 4

Eight TFs with NES higher than 3 and targets more than 150.

Transcription factor	NES	Targets
STAT1	10.871	455
SPIB	6.259	410
NFKB1	4.906	365
SPI1	4.375	336
STAT5A	3.558	199
RUNX1	3.212	196
VENTX	3.126	187
BATF	3.585	158

Fig. 8

TF–miRNA–DEG regulatory network analysis. The network was visualized using Cytoscape. Brown hexagons represent TFs, purple V shapes represent miRNAs, and blue circles represent hub DEGs.

Characterization of immune cell infiltration

The immunocyte compositions of all mucosal biopsies from lung transplant recipients were investigated according to the CIBERSORT algorithm. Twenty-four rejection samples and 64 stable samples that matched the requirements of CIBERSORT P < 0.05 were filtered out (Table S2 ). As a result, memory B cells were not detected in stable samples, while M0 macrophages were not detected in rejection samples. The immunocyte subpopulations of each sample are summarized in Fig. 9 A, and the percentages and subpopulations of immunocytes were identified and visualized using the heatmap in Fig. 9B. A weak-to-moderate correlation between various immunocyte subpopulation fractions was further revealed through the correlation analysis of rejection and stable samples (Fig. 9C, D). In stable samples, activated mast cells had the strongest positive correlation with neutrophils (r = 0.53, P < 0.001), whereas CD8 T cells had the strongest negative correlation with resting memory CD4 T cells (r = −0.63, P < 0.001). In rejection samples, naïve CD4 T cells had the strongest positive correlation with activated dendritic cells (r = 0.65, P < 0.001), whereas activated memory CD4 T cells had the strongest negative correlation with activated NK cells (r = −0.51, P = 0.01).

Table S2

Fig. 9

Immune infiltrate landscape of mucosal tissues from lung transplant recipients. (A) Stacked bar chart representing deviations in immune infiltration in each sample. (B) Heatmap showing the 22 fractions of immunocytes. The X-axis represents samples, and the Y-axis represents immune cell infiltration. (C and D) Correlation matrix of the immunocyte proportions in stable (C) and rejection (D) samples, performed by Pearson correlation analysis. The numbers represent correlation coefficients.

In addition, compared with stable samples, higher proportions of activated memory CD4 T cells, follicular helper T cells, γδ T cells, monocytes, M1 and M2 macrophages, and eosinophils were detected in rejection samples, along with lower proportions of resting memory CD4 T cells, regulatory T cells (Tregs), activated NK cells, M0 macrophages, and resting mast cells (Fig. 10 ).

Fig. 10

Differences in proportions of each immune cell type in stable and rejection samples, analyzed by the two-tailed Student t test. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Analysis of the relationship between hub gene expression and immune infiltration level

Infiltration of immune cell populations has a significant impact on lung transplant rejection. Therefore, this study investigated whether the expression of hub genes correlated with the proportions of immunocytes, especially of immunocytes significantly different between stable and rejection lung transplant samples. As a result, all nine hub genes showed significant correlations with the infiltrating levels of immune cells in rejection lung transplant samples (Fig. 11 ). For example, IL-6 expression positively correlated with the proportion of monocytes (r = 0.82, P < 0.001). ITGAM expression positively correlated with the proportion of eosinophils (r = 0.76, P < 0.001) and negatively correlated with the proportion of Tregs (r = −0.52, P < 0.001). CD86 expression negatively correlated with the proportion of resting mast cells (r = −0.51, P = 0.01). CXCL10 expression positively correlated with the proportion of M1 macrophages (r = 0.59, P < 0.001). CCR5 expression positively correlated with the proportion of activated memory CD4 T cells (r = 0.53, P < 0.001) and γδ T cells (r = 0.50, P = 0.01).

Fig. 11

Correlation index between the expression of hub genes and the immune infiltration level, performed by Pearson correlation analysis. The numbers represent correlation coefficients.

Drug–gene interaction analysis

A total of 66 drug–gene interactions were identified according to DGIdb prediction, including 8 hub genes (PTPRC, IL-6, ITGAM, CD86, TLR8, CXCL10, ITGB2, and CCR5) and 60 kinds of drugs (Fig. 12 ). However, no FDA-approved drugs targeting the hub gene TYROBP were identified. Three hub genes, including ITGB2, IL-6, and CXCL10, were found to be theoretically important in the inhibition of rejection following, and they were targeted by 16, 13, and 9 FDA-approved drugs, respectively. Among these drugs, six drugs (siltuximab, belatacept, abatacept, lifitegrast, phenprocoumon, and maraviroc) were identified as antagonists or inhibitors, and three drugs (ethynodiol diacetate, imiquimod, and glycerin) were identified as agonists. Of note, antagonists or inhibitors might have more significance in lung transplant rejection because all hub genes were upregulated in rejection groups.

Fig. 12

Prediction of drug–gene interactions by DGIdb. The interactions were visualized using Cytoscape. Blue hexagons represent hub genes, and triangles represent drugs. Antagonists or inhibitors are marked in green, and agonists are marked in red.

Discussion

Rejection is one of the main risk factors for limiting lung allograft survival. However, the genetic differences and characteristics of rejection and stable lung transplant samples are still unknown. In the present study, a genomic analysis was used to screen 739 DEGs from the GSE125004 microarray dataset in the GEO. As suggested by the GO and KEGG enrichment analyses, 459 upregulated DEGs were manifested mainly in items related to immune response, suggesting the abnormal immune regulation following LTx [27]. Moreover, the 280 downregulated DEGs were manifested mainly in items related to the cilium. The function of the mucociliary apparatus is impaired due to the ischemia of the airway mucosa following LTx, predisposing patients to repeated pulmonary infections [28], [29].

According to the PPI network, nine DEGs with high degrees were identified as hub genes and upregulated during the rejection of LTx: PTPRC, IL-6, ITGAM, CD86, TLR8, TYROBP, CXCL10, ITGB2, and CCR5. The PTPRC gene encodes the protein tyrosine phosphatase CD45, which acts as a hematopoietic JAK phosphatase required for lymphocyte activation [30]. In a rat model of LTx, an increased proportion of CD4+/CD45RC–T cells is a marker of acute rejection of lung grafts [31]. Instead, CD4+/CD45RC + T cells are able to abrogate lung allograft rejection, and their proportion significantly increases in tolerant rats following LTx [32]. An elevated IL-6 secretion is involved in the acute and chronic rejection of LTx [33], [34] and is linked to acute complications such as PGD [35], [36]. Complement receptor CR3/Mac-1 is composed of CD11b (encoded by ITGAM) and CD18 (encoded by ITGB2), regulating the adhesion of leukocytes and the phagocytosis of complement-coated particles [37]. The inhibition of Mac-1 reduced the recruitment of neutrophils and improved lung graft function in a rat model of severe lung allograft reperfusion injury [38]. CD86 molecules are members of the immunoglobulin superfamily and serve as co-stimulatory molecules on the surface of antigen-presenting cells (APCs). CD86 molecules interact with CD28 on T cells and activate T-cell responses, contributing to lung allograft rejection [39], [40]. TLR8, a Toll-like receptor (TLR), acts as an important sensor for bacteria and viruses and can activate an innate immune response against infection [41], [42]. The upregulation of TLR8 during lung graft rejection may be associated with infection, while innate immunity activated by TLR signaling is involved in the initiation and maintenance of graft rejection [43], [44]. TYROBP (DAP12) encodes a signaling adapter for multiple pattern recognition receptors in myeloid and NK cells [45]. DAP12 can promote lung transplant–mediated ischemia/reperfusion injury through facilitating the extravasation of neutrophils into the pulmonary tissue [46]. The levels of CXCL10, an IFN-γ-induced small pro-inflammatory cytokine, have been found to be elevated in bronchoalveolar lavage (BAL) from patients with restrictive allograft syndrome. The high level of CXCL10 may act as a biomarker that predicts poor lung allograft survival [47]. CCR5, a pro-inflammatory chemokine receptor, is upregulated in BALs and promotes T-cell infiltration during episodes of acute lung rejection [48], [49], [50], [51], [52]. Treatment with antibodies against CCL5 markedly reduces inflammatory events and attenuates acute lung allograft rejection [52].

Significant modules of the PPI network were screened out in this study. Modules 2 and 3 contained five and three hub genes, respectively. The enrichment analysis of Modules 2 and 3 showed that DEGs in these modules were associated mainly with items related to immune responses (cytokine-mediated signaling pathway, adaptive immune response, regulation of cytokine production, and natural killer cell–mediated cytotoxicity), immune diseases (rheumatoid arthritis and graft-versus-host disease), and infections (response to virus, response to interferon-gamma, influenza A, and Staphylococcus aureus infection). Infections have immunological interactions that promote acute rejection and CLAD following LTx, which negatively impacts the function and survival of lung grafts [53], [54]. The results of enrichment analysis were consistent with biopsy characteristics indicating that parts of samples were infected by bacteria, virus, or fungi [13].

MiRNAs have emerged not only as biomarkers for predicting lung allograft rejection but also as potential therapy for promoting long-term airway integrity and graft survival [55]. A total of five miRNAs participating in the regulation of hub genes were predicted using miRNet, including miR-335-5p, miR-26b-5p, miR-124-3p, miR-1-3p, and miR-155-5p. MiR-335-5p sustains the repair function of lung fibroblasts. In cigarette smokers, miR-335-5p is downregulated in lung fibroblasts and tissues and causes lung injury [56]. MiR-26b is a positive regulator of lung inflammatory response by activating alveolar macrophages (AMs) [57]. MiR-124a is highly expressed in the lungs. It exerts a protective effect against acute lung injury by promoting M2 macrophage polarization and suppressing pro-inflammatory responses [58], [59], [60]. In a study investigating the microRNA profiling of airway epithelium from lung transplant recipients, miR-124 was significantly downregulated in acute rejection groups than in stable groups. Consistent with the results of the present study, miR-124 influences a wider repertoire of targets and may play critical roles in lung transplant rejection [61]. MiR-1 is important for the normal functioning of airway smooth muscle cells; its expression is downregulated in response to acute lung injury (ALI) [62]. A meta-analysis of 275 solid organ transplant recipients revealed that the expression of miR-155 was consistent with the dynamic change in acute rejection degree, suggesting it as a biomarker for monitoring the abnormal allograft status in solid organ transplantation [63].

To further investigate the transcriptional regulatory network of lung transplant rejection, an iRegulon analysis was performed, and eight most important TFs were identified: STAT1, SPIB, NFKB1, SPI1, STAT5A, RUNX1, VENTX, and BATF. STAT1 and STAT5A are members of the signal transducers and activators of transcription (STAT) family; their activation can initiate the inflammatory response in ALI [64], [65]. NFKB1 activation is responsible for innate inflammation and remodeling of bronchiolar mucosal epithelium in chronic lung diseases [66]. SPI1 (PU.1) belongs to E-twenty-six family and functions to initiate pulmonary inflammatory cascade through activating AMs [67], [68]. The inhibition of PU.1 blocks the innate immune function of AMs, thus attenuating airway inflammatory disease such as asthma [69]. RUNX1 is a member of the Runt-related transcription factor (RUNX) family and is highly expressed in pulmonary tissues. RUNX1 deficiency caused by ALI or acute respiratory distress syndrome (ARDS) can lead to intense lung inflammation [70], [71]. VENTX is a human homeobox transcriptional factor that accelerates the maturation and pro-inflammatory activation of macrophages and dendritic cells in autoimmune diseases; however, its role in the lungs has not been reported [72], [73]. BATF, a basic leucine zipper TF, has been considered as a target for preventing rejection following LTx because it activates Th17-induced autoimmunity against lung self-antigens (SAgs) and leads to anti-MHC induced rejection [74].

The CIBERSORT algorithm revealed differences in the immune cell infiltration characteristics between stable and rejection LTx samples. Significant differences in the proportions of different subtypes of T cells were consistent with the observations in solid organ transplantation: (1) CD4 + T-cell memory responses were activated during chronic infection following LTx [75]. (2) Activated pro-inflammatory follicular helper T cells promoted the secretion of donor-specific antibodies, resulting in antibody-mediated rejection following kidney transplantation [76]. (3) The levels of anti-inflammatory Tregs decreased in the lung allografts following LTx. Once recovered, these cells prolonged graft survival by suppressing IL-17 production [77], [78]. (4) The levels of pro-inflammatory γδ T cells increased in the lung allografts following LTx, thus promoting IL-17 secretion and contributing to acute rejection and obliterative bronchiolitis (OB) [77], [79].

NK cells can promote lung allograft tolerance through the depletion of donor APCs [80], [81], [82]. Thus, a low percentage of NK cells may contribute to lung allograft rejection. In addition, higher proportions of monocytes, macrophages, and eosinophils were found in rejection groups. Evidence shows that monocytes and macrophages initiate alloimmune responses following LTx, and are involved in complications, including ventilator-induced lung injury, ischemia/reperfusion injury, PGD, and rejection [83], [84]. In mouse models of LTx, eosinophils were found to induce lung allograft acceptance by downregulating T-cell immune responses [85], [86]. However, the role of high proportions of eosinophils during the lung graft rejection of human recipients is unclear.

The correlation analysis revealed the relationships between hub gene expression and immune cell proportions in rejection samples following LTx; some of these were also found in other situations. For example, high levels of IL-6 secreted by monocytes in kidney allograft biopsies following kidney transplantation predict progressive allograft damage [87]. CXCL10 can induce the chemotaxis of pro-inflammatory M1 macrophages and act as a marker of M1 macrophages [88], [89], [90]. CCR5 and CD4 constitutively interact on the plasma membrane of central memory CD4 + T cells, which are highly associated with the rapid progression of acute HIV-1 infection [91], [92].

Finally, FDA-approved drugs that potentially target hub genes were predicted, and six antagonists/inhibitors were found: lifitegrast, siltuximab, belatacept, abatacept, phenprocoumon, and maraviroc. Lifitegrast targets ITGB2 and is applied in the treatment of ocular chronic graft-versus-host disease [93]. Siltuximab is a human–mouse chimeric monoclonal antibody that targets human IL-6 and inhibits the IL-6 pro-inflammatory pathway in pig-to-baboon organ xenotransplantation models [94]. Belatacept is a selective co-stimulation blocker that binds to membrane CD80 and CD86 of APCs; it has been applied in immunosuppression during lung, kidney, and heart transplantation [95], [96], [97]. However, renal insufficiency was found in lung transplant recipients receiving belatacept; hence, multicenter prospective studies are required to evaluate the safety and efficacy of this drug [98], [99]. Abatacept is closely related to belatacept and acts as an immunosuppressant marketed for rheumatic diseases [100]. Phenprocoumon is predicted to target CCR5; however, its overdose may induce diffuse alveolar hemorrhage [101]. Maraviroc, a specific antagonist of CCR5, can prolong cardiac allograft survival in a rhesus monkey cardiac allograft model, suggesting its potential value for lung graft protection [102].

In summary, an integrated analysis based on the mRNA microarray dataset GSE125004 was performed in this study to identify specific mRNAs playing a pivotal role in lung transplant rejection. Besides, a miRNA–TF–DEG regulatory network in lung transplant rejection was predicted and established. Moreover, the immune cell infiltration characteristics of lung transplant samples and the relationship between hub gene expression and immune infiltration levels in lung transplant rejection samples were revealed. Finally, drugs that potentially targeted hub genes were investigated.

Inevitably, the present study had some innate limitations. The sample size was relatively small, reducing the credibility of the miRNA–TF–DEG co-regulatory network and immune cell infiltration analyses. In addition, information on specific patient/biopsy characteristics for the stable and rejection groups is not available. Moreover, the present study lacked further experiments to verify the analysis results as a solid foundation. All in all, this study provided a comprehensive perspective of regulatory mechanism networks underlying lung transplant rejection and also identified potential molecule targets of early diagnosis and treatment for prolonging lung allograft survival.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.