Identifying 4 Novel lncRNAs as Potential Biomarkers for Acute Rejection and Graft Loss of Renal Allograft

Abstract

Acute rejection (AR) after kidney transplant is one of the major obstacles to obtain ideal graft survival. Reliable molecular biomarkers for AR and renal allograft loss are lacking. This study was performed to identify novel long noncoding RNAs (lncRNAs) for diagnosing AR and predicting the risk of graft loss. The several microarray datasets with AR and nonrejection specimens of renal allograft downloaded from Gene Expression Omnibus database were analyzed to screen differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs). Univariate and multivariate Cox regression analyses were used to identify optimal prognosis-related DElncRNAs for constructing a risk score model. 39 common DElncRNAs and 185 common DEmRNAs were identified to construct a lncRNA-mRNA regulatory relationship network. DElncRNAs were revealed to regulate immune cell activation and proliferation. Then, 4 optimal DElncRNAs, ATP1A1-AS1, CTD-3080P12.3, EMX2OS, and LINC00645, were selected from 17 prognostic DElncRNAs to establish the 4-lncRNA risk score model. In the training set, the high-risk patients were more inclined to graft loss than the low-risk patients. Time-dependent receiver operating characteristics analysis revealed the model had good sensitivity and specificity in prediction of 1-, 2-, and 3-year graft survival after biopsy ($\text{AUC}=0.891,0.836,\text{and}0.733$, respectively). The internal testing set verified the result well. Gene set enrichment analysis which expounded NOD-like receptor, the Toll-like receptor signaling pathways, and other else playing important role in immune response was enriched by the 4 lncRNAs. Allograft-infiltrating immune cells analysis elucidated the expression of 4 lncRNAs correlated with gamma delta T cells and eosinophils, etc. Our study identified 4 novel lncRNAs as potential biomarkers for AR of renal allograft and constructed a lncRNA-based model for predicting the risk of graft loss, which would provide new insights into mechanisms of AR.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Kidney transplant is increasingly and widely acceptable treatment for end-stage renal disease. However, acute rejection (AR) that occurs days to months after renal transplantation is one of the major obstacles to obtain ideal long-term graft survival, while plenty of patients with end-stage renal disease are still waiting for renal grafts [1–3]. Therefore, it is important to routinely monitor the function of allograft and diagnose AR promptly. The serum creatinine as a regular biomarker for inferring AR can be late for the change of allograft condition, while kidney biopsy, regarded as the gold standard for diagnosing AR, still has limitations, such as potentially variable artificial diagnosis and others, which reflect that it requires additional molecular biomarkers to become a unified diagnostic criteria [4, 5]. Although the rate of acute rejection has drastically decreased over the past five decades due to the advances in immunosuppressive therapy, some evidence revealed that improvements of traditional immunosuppression in transplant recipients for preventing AR could improve the early outcome of renal allograft but not the long-term graft survival [1, 6]. Moreover, high-dose immunosuppressants such as steroids and/or lymphocyte-depleting antibodies enhance the risk of infection, sepsis, and even cancer, potentially bringing about death of patients with functioning allograft [7–9]. Hence, it is necessary to identify novel biomarkers for AR of renal allograft, especially prognosis-related biomarkers, which probably make available for speculating the underling mechanisms of AR and provide druggable targets in the future.

Long noncoding RNAs (lncRNAs) defined as more than 200 nucleotides RNA without the ability of protein coding were widely investigated and found to affect development of human diseases through regulating neighbouring genes and chromatin topology, scaffolding, and decoying proteins [10, 11]. Their effective roles have been recognized in renal transplantation, renal ischemia, fibrosis, injury, inflammation, glomerular diseases, and renal cell carcinoma [12]. Chen et al. [13] and Zou et al. [14] identified dysregulated lncRNAs in AR specimens when compared with nonrejection specimens. Xu et al. [15] also screened the lncRNAs associated with chronic damage and graft loss after renal transplantation. Besides, lncRNA RP11-354P17.15-001 in urine was reported that it might serve as a novel biomarker for AR of kidney [16]. However, the underlying mechanism through which lncRNAs contribute the rejection or injury is rarely known. Atianand et al. has discovered that lncRNAs take part in development and differentiation of innate and adaptive immune cells [17]. Meanwhile, some studies revealed that allograft-infiltrating immune cells own potential diagnostic and prognostic values in AR after renal transplantation [18, 19], which turns out to be a major target for immunosuppression therapy [20]. To the best of our knowledge, no exact lncRNA has been regarded as the known criteria to help diagnose AR and predict the prognosis of graft survival so far. Thus, lncRNAs related to AR of kidney transplant, graft survival, and allograft-infiltrating immune cells were screened by our study based on the public datasets. This will probably help develop novel biomarkers for AR of renal allograft, directly or indirectly impact current clinical practices, and provide new insights into the mechanism of AR.

2. Material and Methods

2.1. GEO Dataset Acquisition and Preprocessing

Gene Expression Omnibus database (GEO, https://www.ncbi.nlm.nih.gov/geo/) is a public repository of high-throughput gene expression data at the National Center of Biotechnology Information, which is submitted by researchers in various fields, and the data are freely accessible. Five datasets, GSE34437 [21], GSE75693 [22], GSE50058 [23], GSE76882 [24], and GSE21374 [25], containing RNA sequencing data extracted from kidney biopsies from kidney transplant patients with different graft state were selected. Their series matrix and MINiML formatted family files were downloaded. Two platforms were identified: GSE76882 dataset was from Affymetrix HT HG-U133+ PM Array Plate, while others were all from Affymetrix Human Genome U133 Plus 2.0 Array. The data tables of the two platforms were also downloaded to obtain the gene symbol for reannotating the Affymetrix Probe Set ID in the series matrix files. Besides, we extracted the corresponding clinical information from the MINiML formatted family files using SangerBox tool (Version 1.0.9, http://sangerbox.com/). Among these datasets, the diagnosis criteria of AR for allograft were biopsy-proven based on the Banff criteria scored by blinded pathologists, which accompanied rising serum creatinine. Nonrejection or stable biopsies were regarded as stable graft function and absence of significant pathological injury based on the Banff criteria. As for the protocol surveillance or for-cause biopsies, the GSE34437, GSE50058, and GSE76882 datasets contained mixed samples, while the GSE75693 and GSE21374 datasets embrace for-cause biopsies. In our study, biopsy-proven acute rejection, borderline acute rejection, and acute rejection with interstitial fibrosis and tubular atrophy were incorporated into the AR group in our study. Specimens diagnosed with BK virus nephritis in the datasets were excluded. The number of patients with different graft state and outcome in the five datasets is shown in Table 1. The more detailed clinical information of the kidney biopsies is given in Tables S1 and S2.

Table 1

Clinical characteristics of the GEO datasets.

Characteristics	Merged			GSE50058	GSE76882	GSE21374
	GSE34437	GSE75693	Total
Graft state
Acute rejection	17	15	32	43	83	—
Stable	16	30	46	58	99	—
Graft outcome
Failed	—	—	—	—	—	51
Nonfailed	—	—	—	—	—	231
Total	33	45	78	101	182	282

Therefore, GSE34437 and GSE75693 datasets were merged into a single cohort, in the following called the merged dataset, due to the small number of specimens and followed by batch normalization using “sva” R package to eliminate the batch effect. Then in the merged dataset, GSE50058 and GSE76882 datasets were used for differential expression analysis. The time from biopsy to the graft failure or patient censored in GSE21374 dataset was defined as graft survival, which was used for subsequent graft survival analysis of the screened genes. The 282 patients with information of graft survival in GSE21374 dataset were randomly divided into the training set ($n=198$) and the testing set ($n=84$) by a ratio of 7 to 3 using “caret” R package.

2.2. Differential Expression Analyses of lncRNAs and mRNAs

The expression data of GEO datasets were checked whether they were standardized and were quantile normalized using the normalizeBetweenArrays function from “limma” R package (Figure S1). Then, lncRNAs and protein-coding genes were extracted from the matrix based on the human genome (CRCh38) from database Ensembl (http://asia.ensembl.org/index.html). Moreover, the gene was excluded if it is an unrecognized gene or the sum of expression level for each specimen is less than 1. Finally, 1236 lncRNAs and 1234 lncRNAs were obtained from the merged dataset and GSE50058 datasets, respectively. 16413, 16432, and 14064 protein-coding genes were obtained from the merged datasets GSE50058 and GSE76882, respectively.

The “limma” R package was used to screen differentially expressed lncRNAs (DElncRNAs) and differentially expressed mRNAs (DEmRNAs) between kidney transplant patients with AR and stable function. Then, the DElncRNAs from the merged dataset and GSE50058 dataset were intersected to obtain common DElncRNAs, while the DEmRNAs from the merged dataset, GSE50058 and GSE76882 were intersected to obtain common DEmRNAs for subsequent analyses. Adjusted $p$ value < 0.05 and absolute ${\log }_2\text{fold}\text{change}\text{FC}>0.5$ were set to screen DElncRNAs. Adjusted $p$ value < 0.05 and absolute ${\log }_2\text{fold}\text{change}\text{FC}>1$ were used to screen DEmRNAs.

2.3. Construction of lncRNA-mRNA Regulatory Relationship Network

To elucidate the expressional relationship between dysregulated lncRNAs and mRNAs in AR of kidney allograft, the Pearson correlation analysis was performed with the cutoff criteria of absolute coefficient of correlation $r>0.7$ and $p$ value < 0.01. The visualization of the regulatory network was constructed using Cytoscape (V3.7.2, https://cytoscape.org/).

2.4. Functional Annotation and Enrichment Analyses of DElncRNAs

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using “clusterProfiler” R package for revealing the potential gene function terms and enriched signaling pathways of the DEmRNAs in the lncRNA-mRNA regulatory relationship network, which could indirectly indicate the potential biological function of DElncRNAs. GO enrichment analysis contains the aspects of molecular function (MF), cellular components (CC), and biological processes (BP). Before doing the enrichment analyses, gene symbol was converted to Entrez Gene IDs using “http://org.hs.eg.db” R package in order to the mapping between GO or KEGG terms and Entrez Gene IDs. The cutoff value for the enrichment analyses was adjusted $p$ value < 0.05.

2.5. Identification of Prognostic lncRNAs

The normalized expression data of DElncRNAs and corresponding survival time of grafts from the training set ($n=198$) in GSE21374 dataset were used to identify the prognosis-related DElncRNAs using univariate Cox proportional hazards regression analysis. Then, the patients were divided into two groups, high- and low-expression groups by the median expression level as a cutpoint of each prognostic DElncRNA, which was determined by “survminer” R package. Then, “survival” R package was performed to implement a log-rank test and draw Kaplan-Meier curves to compare the graft survival rate between high- and low-expression levels of the prognostic DElncRNAs [26]. Afterwards, the prognostic DElncRNAs screened from the univariate Cox proportional hazard regression were enrolled in stepwise multivariate Cox proportional hazard regression method to select the optimal gene model with the minimum AIC value. A $p$ value < 0.05 was considered as significant.

2.6. Establishment and Assessment of Multi-lncRNA Model

The final optimal prognosis-related DElncRNAs screened from multivariate Cox regression were used to construct a risk score model to evaluate graft survival in renal transplant patients who were diagnosed with acute rejection. The risk score of each patient was calculated by multiplication of the regression coefficient of the lncRNAs obtained from multivariate Cox regression and their expression level and finally summing, which was the following: $$\begin{array}{}\text{(1)}& \text{Risk}\text{score}=\sum _{i=1}^n{\beta }_i\ast {\text{Exp}}_i,\end{array}$$where $n$, $\beta$, and $\text{Exp}$ are the number, the regression coefficient, and the expression value of the prognosis-related DElncRNAs, respectively.

Then, the median risk score was regarded as the cutoff point to stratify the patients into high-risk and low-risk groups in the training set. The Kaplan-Meier (KM) survival analysis was performed to compare the graft survival between these two groups. Moreover, the time-dependent receiver operating characteristics (tROC) curve was used to evaluate the specificity and sensitivity of the graft survival prediction using “survivalROC” R package.

In order to assess the prediction value of the model, the testing set ($n=84$) in GSE21374 dataset was used to perform validation. The risk score of each renal transplant patient was calculated through the above formula with the same coefficients used in the training set. Then, the patients were stratified into high-risk and low-risk groups using the same cutoff point in the training set. Besides, we used KM survival analysis and tROC analysis to validate the model in the testing set.

2.7. Immune Cell Infiltration Analysis and Gene Set Enrichment Analysis of Final Optimal lncRNAs

CIBERSORT (cell type identification by estimating relative subsets of RNA transcripts, https://cibersort.stanford.edu/) [27] is a computational approach for characterizing cell composition of complex tissues including fresh, frozen, and fixed tissues from the gene expression profiles. It outperformed many other estimating methods with respect to noise and unknown mixture content [28]. A signature matrix of the collated 547-gene expression datasets associated with 22 immune cell types was downloaded from CIBERSORT web portal, which encompasses T cells, B cells, natural killer cells, dendritic cells, eosinophils, etc. Then, the significantly different proportions of infiltrating immune cells in the merged dataset and GSE76882 dataset between acute rejection and nonrejection specimens were calculated using the CIBERSORT algorithm and Wilcoxon method for variance by R software with the criteria of $p$ value < 0.05. The correlation between common infiltrating immune cells and expression of the final optimal lncRNAs was further analyzed using Pearson correlation analysis, which considered $\mid r\mid >0.3$ and $p$ value < 0.01 as significant.

For identifying the biological pathways of the final optimal lncRNAs, gene set enrichment analysis (GSEA) [29] was performed on JAVA 8.0 platform using GSEA software (V4.0.4, http://software.broadinstitute.org/gsea/). The median expression of the lncRNAs was regarded as a cutoff point to divide all the samples into high- and low-expression groups. Then, the annotated gene sets “c2.cp.kegg.v7.0.symbols” and “c5.all.v7.0.symbols” downloaded from the Molecular Signatures Database (MSigDB) were adopted as the reference gene set to calculate enrichment score (ES). The number of permutations was 1000. $\text{Gene}\text{size}<15$ or >500 was excluded. The enrichment pathways with $\mid \text{ES}\mid >1.5$, normalized $p$ value < 0.05, and $\text{FDR}<0.25$were expected as significant [29].

2.8. Statistical Analysis

R software (V3.6.1, The R Foundation for Statistical Computing, 2019) was used to perform all statistical analysis. Volcano plots of DElncRNAs and DEmRNAs were plotted using “ggrepel” R package, while heat maps of DElncRNAs and DEmRNAs were plotted using “pheatmap” R package with zero-mean normalization. The expressional boxplots between the AR and NR groups were analyzed using the Mann–Whitney $U$ test. The boxplots and correlation graphs were drawn by Prism 8 (GraphPad). To validate the association between the final optimal lncRNAs and AR of renal allograft, logistic regression analysis was conducted and area under the ROC curve (AUC) was also calculated to illustrate their diagnostic accuracy. For Kaplan-Meier curves, $p$ values and hazard ratio (HR) with 95% confidence interval (CI) were generated by the log-rank tests and Cox regression methods. All statistical tests were two-sided. $p$ value < 0.05 was considered as statistically significant.

3. Results

3.1. Identification of DElncRNAs and DEmRNAs in AR

The flow chart of our whole study is shown in Figure 1.

[figure omitted; refer to PDF]

After differential expression analysis via limma R package, there were 134 DElncRNAs and 553 DEmRNAs in the merged dataset; 63 DElncRNAs and 439 DEmRNAs in GSE50058 dataset; and 620 DEmRNAs in GSE76882. The cluster heat maps of top 20 DElncRNAs and volcano plots of DElncRNAs in the merged dataset and GSE50058 dataset are shown in Figures 2(a)–2(d). Then, the common 39 DElncRNAs and 185 DEmRNAs were obtained via integration using “vennDiagram” R package, which are shown in Figures 2(e) and 2(f).

[figures omitted; refer to PDF]

3.2. lncRNA-mRNA Regulatory Relationship Network Construction and Functional Enrichment Analysis

For identifying the coexpression relationship between DElncRNAs and DEmRNAs in AR of renal allograft, the Pearson correlation analysis was performed based on GSE50058 dataset and totally there were 101 lncRNA-mRNA pairs with positive correlation and 16 lncRNA-mRNA pairs with negative correlation ($\mid r\mid >0.7$, $p$ value < 0.05). The relevant positive and negative correlation networks were constructed using Cytoscape and they are shown in Figures 3(a) and 3(b). Besides, in order to figure out the potential biological function of common DElncRNAs, the coexpression DEmRNAs were used for GO and KEGG analyses. The significantly enriched GO terms are shown in Figure 3(c). Top 10 enriched terms in BP, CC, and MF are listed in the figure. In BP category, most of the enriched terms were associated with regulation of immune cell activation and proliferation. In CC and MF category, “external side of plasma membrane” and “G protein-coupled receptor binding” were the most enriched terms, respectively. The significantly enriched KEGG pathways are shown in Figure 3(d). “Cytokine-cytokine receptor interaction” was the most enriched pathway. “Allograft rejection” and “MHC protein complex binding” were also enriched which demonstrated the good reliability of the differential expression analysis screened by our study. Moreover, signaling pathways like “NOD-like receptor signaling pathway” and “Toll-like receptor signaling pathway” were enriched. Pathways involved in “Allograft rejection” and “Graft-versus-host disease” were enriched as well.

[figures omitted; refer to PDF]

3.3. Screening of Prognosis-Related DElncRNAs

The graft survival data after biopsy of renal allograft and the expression data of 39 common DElncRNAs in the training set were subjected to univariate Cox proportional hazards regression analysis with the significant threshold of $p$ value < 0.05. Therefore, 15 DElncRNAs associated with the graft survival were identified, among which there were 2 positive DElncRNAs ($\text{HR}>1$) and 13 negative DElncRNAs ($\text{HR}<1$). Then, these 17 DElncRNAs were used in subsequent stepwise multivariate Cox proportional hazards regression analysis. Finally, the optimal 4-lncRNA combination was obtained with the minimum Akaike’s information criterion value ($\text{AIC}=317.29$), which contained ATP1A1-AS1, CTD-3080P12.3, EMX2OS, and LINC00645. The results of Cox regression analysis are shown in Table 2.

Table 2

Univariate and multivariate Cox regression of 39 DElncRNAs in the training set of GSE21374.

Variable	Univariate			Multivariate
	HR	95% CI	$p$ value	HR	95% CI	$p$ value
TRG-AS1	1.951	1.426-2.688	$$<0.001$$	—	—	—
LINC00645	0.58	0.443-0.758	$$<0.001$$	0.549	0.396-0.761	$$<0.001$$
LINC01187	0.461	0.305-0.697	$$<0.001$$	—	—	—
TCL6	0.202	0.085-0.480	$$<0.001$$	—	—	—
DANCR	0.25	0.116-0.540	$$<0.001$$	—	—	—
LINC00982	0.437	0.270-0.707	$$0.001$$	—	—	—
CTD-3080P12.3	0.159	0.050-0.504	$$0.002$$	0.102	0.021-0.505	$$0.005$$
EMX2OS	0.212	0.078-0.576	$$0.002$$	0.324	0.091-0.916	0.082
TRAM2-AS1	0.236	0.090-0.620	$$0.003$$	—	—	—
LINC00671	0.502	0.309-0.817	$$0.005$$	—	—	—
WAC-AS1	3.154	1.414-7.036	$$0.005$$	—	—	—
ATP1A1-AS1	0.106	0.021-0.537	$$0.007$$	0.14	0.020-0.968	$$0.046$$
AC112198.1	0.367	0.174-0.773	$$0.008$$	—	—	—
WDFY3-AS2	0.682	0.507-0.916	$$0.011$$	—	—	—
LINC00886	0.403	0.183-0.892	$$0.025$$	—	—	—
AL022344.5	3.228	1.129-9.231	$$0.029$$	—	—	—
RPARP-AS1	0.242	0.066-0.890	$$0.033$$	—	—	—
SMIM2-AS1	0.468	0.205-1.071	$$0.072$$	—	—	—
RUSCI-AS1	2.547	0.879-7.385	$$0.085$$	—	—	—
C12orf77	7.544	0.712-79.89	$$0.093$$	—	—	—
ITGB2-AS1	1.338	0.952-1.880	$$0.094$$	—	—	—
PCED1B-AS1	2.276	0.866-5.986	$$0.095$$	—	—	—
LINC00592	3.047	0.810-11.463	$$0.099$$	—	—	—
ADIRF-AS1	0.477	0.197-1.154	$$0.101$$	—	—	—
AC005523.2	0.391	0.099-1.533	$$0.177$$	—	—	—
ZNF213-AS1	0.416	0.110-1.570	$$0.196$$	—	—	—
TRHDE-AS1	0.769	0.479-1.235	$$0.277$$	—	—	—
LINC00911	2.155	0.377-12.30	$$0.387$$	—	—	—
AC092192.1	1.494	0.534-4.418	0.446	—	—	—
ELOVL2-AS1	2.084	0.316-13.740	$$0.446$$	—	—	—
GS1-124K5	2.14	0.247-18.58	$$0.489$$	—	—	—
EPB41L4A-AS1	0.651	0.182-2.325	$$0.508$$	—	—	—
APTR	0.761	0.221-2.623	$$0.666$$	—	—	—
LINC01410	0.601	0.053-6.847	$$0.682$$	—	—	—
HMMR-AS1	0.698	0.104-4.654	$$0.711$$	—	—	—
LINC00472	0.945	0.639-1.398	$$0.778$$	—	—	—
LINC01222	0.776	0.114-5.259	$$0.795$$	—	—	—
FLJ37453	0.934	0.145-6.001	$$0.943$$	—	—	—
SEPSES-AS1	1.033	0.335-3.179	$$0.955$$	—	—	—

HR: hazard ratio; CI: confidence interval.

3.4. Expression Profiles and Survival Analysis of the Optimal 4 lncRNAs

The expression profiles of the optimal 4 DElncRNAs between AR and NR specimens in the merged dataset and GSE50058 dataset are presented in Figure 4(a), which indicated that the 4 DElncRNAs were significantly downregulated in AR of renal allograft ($p<0.001$). Moreover, the result of logistic regression analysis based on the merged dataset and GSE50058 dataset verified that the 4 DElncRNAs were significantly associated with the AR of renal allograft ($p<0.01$) (Table 3).

[figures omitted; refer to PDF]

Table 3

Logistic regression analyses and ROC analyses of the 4 optimal lncRNAs.

Predictors	Merged GSEa ($n=78$)		GSE50058 ($n=101$)
	OR	AUC (SE)	OR	AUC (SE)
ATP1A1-AS1	0.032${}^{\ast \ast \ast }$	0.863 (0.041)${}^{\ast \ast \ast }$	0.123${}^{\ast \ast \ast }$	0.780 (0.048)${}^{\ast \ast \ast }$
CTD-3080P12.3	0.147${}^{\ast \ast \ast }$	0.852 (0.044)${}^{\ast \ast \ast }$	0.160${}^{\ast \ast \ast }$	0.744 (0.049)${}^{\ast \ast \ast }$
EMX20S	0.147${}^{\ast \ast \ast }$	0.750 (0.060)${}^{\ast \ast \ast }$	0.501${}^{\ast \ast }$	0.666 (0.057)${}^{\ast \ast }$
LINC00645	0.232${}^{\ast \ast \ast }$	0.756 (0.059)${}^{\ast \ast \ast }$	0.455${}^{\ast \ast \ast }$	0.718 (0.052)${}^{\ast \ast \ast }$

^aMerged $\text{GSE}=\text{GSE}34437+\text{GSE}75693$. ${}^{\ast \ast \ast }p$ value < 0.001. ${}^{\ast \ast }p$ value < 0.01. ROC: receiver operating characteristics; OR: odds ratios; AUC: area under the curve; SE: standard error.

In the meanwhile, ROC analysis was performed as well to obtain their AUC values and standard error. In the merged dataset, $\text{AU}{\text{C}}_{\text{ATP}1\text{A}1‐\text{AS}1}=0.863$ (0.041), $\text{AU}{\text{C}}_{\text{CTD}‐3080\text{P}12.3}=0.852$ (0.044), $\text{AU}{\text{C}}_{\text{EMX}20\text{S}}=0.750$ (0.060), and $\text{AU}{\text{C}}_{\text{LINC}00645}=0.756$ (0.059). In GSE50058 dataset, $\text{AU}{\text{C}}_{\text{ATP}1\text{A}1‐\text{AS}1}=0.780$ (0.048), $\text{AU}{\text{C}}_{\text{CTD}‐3080\text{P}12.3}=0.744$ (0.049), $\text{AU}{\text{C}}_{\text{EMX}20\text{S}}=0.6660.057$, and $\text{AU}{\text{C}}_{\text{LINC}00645}=0.718$ (0.052).

KM survival curves were plotted to evaluate the prognostic value of the 4 DElncRNAs (Figure 4(b)), in which the median expression value of each DElncRNA was regarded as a cutoff point to partition the patients into the high-expression and low-expression groups. Low expression of the 4 DElncRNAs was associated with the poor prognosis of the renal allograft. The KM survival curves of other 11 prognosis-related DElncRNAs screened from univariate Cox regression analysis are also shown in Figure S3.

3.5. Construction and Validation of the 4-lncRNA Model

The Cox coefficients of the 4 lncRNAs obtained from the multivariate Cox proportional hazards regression analysis were used to multiply their expression values for calculating the risk score of each patient in GSE21374 dataset. $\text{Risk}\text{score}=-1.9663\ast \text{Ex}{\text{p}}_{\text{ATP}1\text{A}1‐\text{AS}1}+-2.28396\ast \text{Ex}{\text{p}}_{\text{CTD}‐3080\text{P}12.3}+-1.12712\ast \text{Ex}{\text{p}}_{\text{EMX}20\text{S}}+-0.59929\ast \text{Ex}{\text{p}}_{\text{LINC}00645}.$

Then, the patients in the training set were divided into high- and low-risk groups based on the median risk score determined by “survminer” R package, in which the risk score curve is shown in Figure 5(a). The corresponding graft survival status of the patients in the training set is demonstrated in Figure 5(b), which suggests that there were more patients who were classified as high risk tend to get graft failure. The expression profiles of the 4 DElncRNAs in AR and stable specimens are visualized in the cluster heat map (Figure 5(c)). The KM survival curve revealed that the high-risk group patients have worse graft survival outcome than that of the low-risk group patients ($p<0.001$) (Figure 5(d)). Moreover, tROC analysis was performed and AUC values of 1-, 2-, and 3-year graft survival after biopsy are 0.891, 0.836, and 0.733 (Figure 5(e)).

[figures omitted; refer to PDF]

In order to verify the results obtained from the training set, the testing set was used for the same analysis. The cutoff point of the risk score in the training set was used to divide the patients into high- and low-risk groups in the testing set (Figure 6(a)), and more patients got graft failure in the high-risk groups (Figure 6(b)). The cluster heat map of the 4 DElncRNAs is also presented in Figure 6(c). The KM survival curves verified the significant difference of graft survival between the high- and low-risk groups ($p<0.05$) (Figure 6(d)). AUC values of 1-, 2-, and 3-year graft survival after biopsy in the testing set are 0.805, 0.781, and 0.763 (Figure 6(e)).

[figures omitted; refer to PDF]

3.6. Correlation between the Optimal 4 lncRNAs and Immune Cell Infiltration

For determining the types of infiltrating immune cells involved in AR of renal allograft, CIBERSORT algorithm was performed to speculate the percentage of 22 types of immune cells in allograft specimens. After the calculation and filter with the criteria of $p$ value < 0.05, the merged dataset and GSE76882 contained 68 and 108 AR and stable specimens, while GSE50058 contained only 15 AR and stable specimens. Therefore, the GSE50058 dataset was excluded in the comparison analysis for revealing significant different infiltrated immune cells between AR and stable specimens.

The comparison of immune cell infiltration between AR and stable specimens in the merged dataset and GSE76882 dataset is shown in Figures 7(a) and 7(b), respectively. In these two datasets, there were 5 common significantly different immune cell types: CD8⁺ T cells ($p=0.018$), activated memory CD4⁺ T cells ($p=0.001$), gamma delta T cells ($p=0.028$), and eosinophils ($p=0.011$) were more infiltrated in AR specimens, while resting dendritic cells ($p=0.005$) was less infiltrated in AR specimens compared with stable specimens. Then, the correlation between the optimal 4 lncRNAs and the 5 common significantly different immune cell types is explored and presented in Figures 7(c)–7(f) based on the merged dataset. The expressions of ATP1A1-AS1, CTD-3080P12.3, and EMX2OS were all negatively correlated with the gamma delta T cells and eosinophils. The expression of LINC00645 was negatively correlated with CD8⁺ T cells and positively correlated with resting dendritic cells.

[figures omitted; refer to PDF]

3.7. Gene Set Enrichment Analysis

The potential biological functions of ATP1A1-AS1, CTD-3080P12.3, EMX2OS, and LINC00645 were mined using GSEA. As the results shown in Figure 8(a), lymphocyte differentiation and mitotic cell cycle checkpoint gene sets were found to be enriched in the low expression of ATP1A1-AS1. G1 DNA damage checkpoint and mitotic cell cycle checkpoint gene sets were enriched in the low expression of CTD-3080P12.3 (Figure 8(b)). For the low expression of EMX2OS, DNA-binding transcription factor and RNA splicing gene sets were enriched (Figure 8(c)). Fucosylation and nucleobase biosynthetic process gene sets were enriched in the low expression of LINC00645 (Figure 8(d)). Moreover, several important signaling pathways which play important roles in allograft rejection were found to enriched by the low expression of these 4 lncRNAs, such as the MAPK signaling pathway, MTOR signaling pathway, Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, p53 signaling pathway, TGF-β signaling pathway, JAK-STAT signaling pathway, and VEGF signaling pathway, where the findings provided good clues for exploring the potential specific function of these four lncRNAs in AR of renal allografts.

[figures omitted; refer to PDF]

4. Discussion

Mining high-throughput microarray datasets which contain thousands of genes in small sample sizes are the cost-effective manner in exploring the fields of genetic characteristics and functional genomics [30]. lncRNAs as noncoding RNAs have been reported that they could be the main regulators of immune response and exert their functional roles in immune-mediated tissue rejection [31]. Therefore, we used the public microarray datasets to investigate the potential lncRNAs as biomarkers and therapeutic targets for AR after renal transplantation.

In our study, several GEO datasets about AR after renal transplantation were used to identify differentially expressed lncRNAs and mRNAs to construct the DElncRNA-DEmRNA regulatory relationship network based on Pearson correlation analysis. GO and KEGG functional enrichment analyses were performed to reveal the potential biological function of DElncRNAs, showing that T cell activation was mostly enriched. Proliferation and regulation of immune-related cells including leukocytes, lymphocytes, and mononuclear cells were also enriched. These biological processes were achieved probably by the means of interfering cytokine-cytokine receptor interaction pathway and chemokine signaling pathway, which were enriched in KEGG analysis. This suggests that the identified DElncRNAs are deserved to be further explored for the roles in immune response in the future. Then, the DElncRNAs were used for univariate and stepwise multivariate Cox regression analysis and 4 optimal DElncRNAs associated with graft survival were obtained, which included ATP1A1-AS1, EMX2OS, CTD-3080P12.3, and LINC00645. ATP1A1-AS1 is antisense RNA 1 of Na/K-ATPase α1 (ATP1A1). In human kidney cells, ATP1A1-AS1 was discovered to negatively regulate its sense gene, ATP1A1, which regulates renal cell survival as a signal transducer [32]. Fan also claimed that ATP1A1-AS1 is a moderate negative regulator of ATP1A1 and modulates Na/K-ATPase-related signaling pathways, which play an important role in cardiac fibrosis [33]. It was also found to be dysregulated expression in cutaneous melanoma [34] and related to the prognosis of thymoma [35]. However, no more studies about the molecular mechanisms and signaling pathways of ATP1A1-AS1 can be found. According to our allograft-infiltrating immune cells analysis, we found that ATP1A1-AS1 is correlated negatively with the infiltration of gamma delta T cells and eosinophils. The role of gamma delta T cells in transplantation is underresearched. Some studies have proven that they frequently present in acute rejection of renal allografts and have direct cytolytic activity against renal epithelium, in which the mechanism for killing allogeneic renal cells is a natural killer-like way [36, 37]. Besides, in small animal models, it was observed that gamma delta T cells could produce interleukin-17 (IL-17) to contribute acute and chronic allograft dysfunction in skin [38] and lung transplantation [39]. In contrast, some evidence showed they could produce IL-4 and IL-10 to decrease Th1 responses to achieve allograft protection [40]. For eosinophils, they are regarded as the promoting factor for inducing acute allograft rejection and the increased presence of eosinophils in peripheral blood and/or renal allograft biopsy specimen would be risky factors for outcome of acute rejection [41, 42]. Therefore, according to what we found, low expression of ATP1A1-AS1 was associated with the high infiltration of gamma delta T cells and eosinophils, which probably provide some hints to reveal the causes for acute rejection of the renal allograft. Besides, our GSEA revealed that ATP1A1-AS1 is involved in allograft rejection and other several signaling pathways including the MAPK signaling pathway [43], MTOR signaling pathway [44], and Toll-like receptor signaling pathway [45], which were discovered to play a role in rejection of renal allograft. EMX2OS is an antisense transcript of homeodomain-containing transcription factor EMX2, which plays a vital role in brain development [46] and urogenital development [47]. Gu et al. revealed the relationship between the downregulation of EMX2OS and poor prognosis of classical papillary thyroid cancer [48], while Tang et al. revealed its participation in the molecular mechanisms regulating recurrent laryngeal cancer [49]. However, there have been little related researches investigating the roles of EMX2OS in kidney cells and even renal allograft with AR. As our study results indicated, downregulation of EMX2OS can predict the poor renal allograft survival after biopsy. It also correlates negatively with gamma delta T cells and eosinophils, which can be risky factors in AR of renal allograft. GSEA suggests that EMX2OS can probably exert its roles through regulation of protein binding, RNA splicing, and transcription factor binding. Besides, low expression of EMX2OS also promotes leukocyte transendothelial migration and the TGF-β signaling pathway which regulates cell development, differentiation, apoptosis, and other functions associated with cell homeostasis [50]. As for CTD-3080P12.3, there are not much research introducing its basic information and studying its biological function. Kim et al. regarded CTD-3080P12.3 as a candidate biomarker for thyroid cancer by analyzing their collected samples and The Cancer Genome Atlas dataset [51]. Esposti et al. found it was a differentially expressed lncRNA in hepatocellular carcinoma compared with adjacent cirrhotic tissues [52]. However, the role of CTD-3080P12.3 in rejection or immune response remains unknown. Our results provided some potential clues that CTD-3080P12.3 can be an independent prognostic factor for predicting the renal graft survival as shown in Table 2. The low expression of CTD-3080P12.3 is also associated with the high infiltration of gamma delta T cells and eosinophils in renal allograft with AR. Then on the basis of the functional enrichment analysis, the downregulation of CTD-3080P12.3 possibly gives positively impact on the JAK-STAT signaling pathway, NOD-like receptor signaling pathway, p53 signaling pathway, and Toll-like receptor signaling pathway. The JAK-STAT signaling pathway has the cardinal role in the development and/or function of immune cells [53], while the NOD-like receptor and Toll-like receptor signaling pathways play critical role in activating innate and adaptive immune response [54]. These findings suggest that CTD-3080P12.3 can exert its important role in immune-related aspects. LINC00645, a long intergenic noncoding RNA 645 located in human chromosome 14, is firstly discovered and found to be oncogenic in endometrial cancer by the team of Chen et al. [55]. They found that it was upregulated in endometrial cancer through sequencing the lncRNA transcriptome. Besides, it was studied thoroughly by Li et al., who explored that LINC00645 plays an oncogenic role in glioma through LINC00645/miR-205-3p/ZEB1 signaling axis triggered by the TGF-β signaling pathway [56]. However, the action of LINC00645 in renal allograft rejection or immune response remains unclear. Based on our findings, the low expression of LINC00645 is associated with the high infiltration of CD8⁺ T cells and low infiltration of resting dendritic cells. In grafts undergoing acute rejection, invading CD8⁺ T cells with immunologic specificity for the allograft can release perforin and granzymes A and B to perforate target cell membrane and induce caspase-mediated apoptosis of tubular cells [57]. Dendritic cells are critical in the induction of T cell immunity and in peripheral T cell tolerance. It has been reported that resting dendritic cells can induce peripheral CD8⁺ T cells tolerance through PD-1 and CTLA-4 [58, 59]. These discoveries verify what we found for LINC00645, whose low expression will probably whittle the tolerance of CD8⁺ T cells induced by resting dendritic cells and facilitate AR of renal allograft. Besides, our GSEA revealed that “allograft rejection” was enriched by the low expression of LINC00645, while the JAK-STAT signaling pathway and p53 signaling pathway are also enriched, which suggests that it is quite possible to play an important role in AR through these important signaling pathways. From the above, the researches about specific function and mechanism of these four novel lncRNAs remain limited. Our findings provide some important insights for the lncRNAs and more basic researches are recommended to be achieved to investigate their precise role in acute rejection of renal allograft.

After identifying the optimal 4 lncRNAs, logistic regression analysis showed their significant association with acute kidney transplant rejection and their AUC values of ROC analysis interpreted that they have good sensitivity and specificity for differentiating AR of renal allograft (Table 3). Besides, we also constructed a prognosis-related risk score system on the foundation of the 4 lncRNAs to predict the risk of renal graft loss after biopsy. Through calculating the risk score of transplant patients using the formula we provided in the result, the patients in the training set who were classified into the high-risk group were easier to get renal allograft loss than those into the low-risk group. The same significant result was revealed in the testing set. Time-dependent ROC analysis of the risk score system in the training set and the testing set both illustrated the good performance in prediction of 1-, 2-, and 3-year graft survival after biopsy. Therefore, our 4-lncRNA based model may help predict the graft survival outcome of kidney transplant patients and provide the reference for therapeutic guidance in AR of kidney transplant.

However, several limitations should be acknowledged in our studies. Firstly, to the best of our knowledge, there is only one GEO dataset (GSE21374 dataset) which contains survival time of renal graft after biopsy. Thus, in our study, we divided this dataset randomly into the training set and the testing set. The testing set was used to verify the performance of the model screened from the training set. More independent datasets are needed for validation. Secondly, the datasets we used in our study are mainly from North America, which suggests that more public datasets should get involved to eliminate the geographical difference. Thirdly, the optimal biomarkers for AR and graft loss in our study were screened through mining the public dataset merely. The essential clinical parameters supplied by these public datasets are restricted and it is not able to investigate the potential influence of the clinical parameters on the four lncRNAs and the model. Probably in the future we will plan to collect our own data exploring the effects to elaborate the model. According to what we have found, it is still valuable and crucial to conduct the functional experiment in the future to explore the potential mechanism of these ideal four lncRNAs in acute rejection, which probably assists its transition from experimental research to clinical implementation.

5. Conclusion

In summary, we used high-throughput microarray datasets to identify 4 novel lncRNAs as biomarkers for AR of renal allograft and construct a lncRNA-based risk score model to predict the risk of graft loss after renal transplantation, which probably aids clinicians in choosing or adjusting therapy of immunosuppression after renal transplantation. Besides, functional enrichment analysis and allograft-infiltrating immune cell analysis revealed the value and importance of the exploring experiment in these 4 novel lncRNAs.

Authors’ Contributions

Z.Z. and Y.T. initiated and organized the study. Z.Z., Y.T., and H.Z. selected the datasets, carried out the bioinformatics analyses and statistical analyses, and drew the figures. E.L., L.X., and X.F. helped embellish the figures and tables. Z.Z. and Y.T. wrote the manuscript. J.Z., Y.L., J.L., and Y.Y. reviewed and edited the manuscript. All authors read and approved the final version of the manuscript.

Glossary

Abbreviations

AR:Acute rejection

lncRNAs:Long noncoding RNAs

DElncRNAs:Differentially expressed long noncoding RNAs

DEmRNAs:Differentially expressed mRNAs

GEO:Gene Expression Omnibus

KM:Kaplan-Meier

FDR:False discovery rate

FC:Fold change

HR:Hazard ratio

CI:Confidence interval

ROC:Receiver operating characteristic

GSEA:Gene set enrichment analysis

CIBERSORT:Cell type identification by estimating relative subsets of RNA transcripts

GO:Gene Ontology

KEGG:Kyoto Encyclopedia of Genes and Genomes.

References

[1] A. Hart, J. Smith, M. Skeans, S. K. Gustafson, A. R. Wilk, S. Castro, J. Foutz, J. L. Wainright, J. J. Snyder, B. L. Kasiske, A. K. Israni, "OPTN/SRTR 2018 annual data report: kidney," American Journal of Transplantation, vol. 20 no. s1, pp. 20-130, DOI: 10.1111/ajt.15672, 2020.

[2] P. R. Bicalho, L. R. Requião-Moura, É. F. Arruda, R. Chinen, L. Mello, A. P. F. Bertocchi, E. Lamkowski Naka, E. J. Tonato, A. Pacheco-Silva, "Long-term outcomes among kidney transplant recipients and after graft failure: a single-center cohort study in Brazil," BioMed Research International, vol. 2019,DOI: 10.1155/2019/7105084, 2019.

[3] Z. M. el-Zoghby, M. D. Stegall, D. Lager, W. K. Kremers, H. Amer, J. M. Gloor, F. G. Cosio, "Identifying specific causes of kidney allograft loss," American Journal of Transplantation, vol. 9 no. 3, pp. 527-535, DOI: 10.1111/j.1600-6143.2008.02519.x, 2009.

[4] M. Haas, A. Loupy, C. Lefaucheur, C. Roufosse, D. Glotz, D. Seron, B. J. Nankivell, P. F. Halloran, R. B. Colvin, E. Akalin, N. Alachkar, S. Bagnasco, Y. Bouatou, J. U. Becker, L. D. Cornell, J. P. Duong van Huyen, I. W. Gibson, E. S. Kraus, R. B. Mannon, M. Naesens, V. Nickeleit, P. Nickerson, D. L. Segev, H. K. Singh, M. Stegall, P. Randhawa, L. Racusen, K. Solez, M. Mengel, "The Banff 2017 Kidney Meeting Report: revised diagnostic criteria for chronic active T cell–mediated rejection, antibody-mediated rejection, and prospects for integrative endpoints for next-generation clinical trials," American Journal of Transplantation, vol. 18 no. 2, pp. 293-307, DOI: 10.1111/ajt.14625, 2018.

[5] D. Perkins, M. Verma, K. J. Park, "Advances of genomic science and systems biology in renal transplantation: a review," Seminars in Immunopathology, vol. 33 no. 2, pp. 211-218, DOI: 10.1007/s00281-011-0243-2, 2011.

[6] H. U. Meier-Kriesche, J. D. Schold, T. R. Srinivas, B. Kaplan, "Lack of improvement in renal allograft survival despite a marked decrease in acute rejection rates over the most recent era," American Journal of Transplantation, vol. 4 no. 3, pp. 378-383, DOI: 10.1111/j.1600-6143.2004.00332.x, 2004.

[7] B. Jamil, K. Nicholls, G. J. Becker, R. G. Walker, "Impact of acute rejection therapy on infections and malignancies in renal transplant recipients," Transplantation, vol. 68 no. 10, pp. 1597-1603, DOI: 10.1097/00007890-199911270-00027, 1999.

[8] E. Yanik, S. Gustafson, B. Kasiske, A. K. Israni, J. J. Snyder, G. P. Hess, E. A. Engels, D. L. Segev, "Sirolimus use and cancer incidence among US kidney transplant recipients," American Journal of Transplantation, vol. 15 no. 1, pp. 129-136, DOI: 10.1111/ajt.12969, 2015.

[9] M. Hortlund, L. S. Arroyo Mühr, H. Storm, G. Engholm, J. Dillner, D. Bzhalava, "Cancer risks after solid organ transplantation and after long-term dialysis," International Journal of Cancer, vol. 140 no. 5, pp. 1091-1101, DOI: 10.1002/ijc.30531, 2017.

[10] J. D. Ransohoff, Y. Wei, P. A. Khavari, "The functions and unique features of long intergenic non-coding RNA," Nature Reviews Molecular Cell Biology, vol. 19 no. 3, pp. 143-157, DOI: 10.1038/nrm.2017.104, 2018.

[11] I. Ulitsky, D. P. Bartel, "lincRNAs: genomics, evolution, and mechanisms," Cell, vol. 154 no. 1, pp. 26-46, DOI: 10.1016/j.cell.2013.06.020, 2013.

[12] H. Moghaddas Sani, M. Hejazian, S. M. Hosseinian Khatibi, M. Ardalan, S. Zununi Vahed, "Long non-coding RNAs: an essential emerging field in kidney pathogenesis," Biomedicine & Pharmacotherapy, vol. 99, pp. 755-765, DOI: 10.1016/j.biopha.2018.01.122, 2018.

[13] W. Chen, W. Peng, J. Huang, X. Yu, K. Tan, Y. Chen, X. Lin, D. Chen, Y. Dai, "Microarray analysis of long non-coding RNA expression in human acute rejection biopsy samples following renal transplantation," Molecular Medicine Reports, vol. 10 no. 4, pp. 2210-2216, DOI: 10.3892/mmr.2014.2420, 2014.

[14] Y. Zou, W. Zhang, H. H. Zhou, R. Liu, "Analysis of long noncoding RNAs for acute rejection and graft outcome in kidney transplant biopsies," Biomarkers in Medicine, vol. 13 no. 3, pp. 185-195, DOI: 10.2217/bmm-2018-0272, 2019.

[15] J. Xu, J. Hu, H. Xu, H. Zhou, Z. Liu, Y. Zhou, R. Liu, W. Zhang, "Long non-coding RNA expression profiling in biopsy to identify renal allograft at risk of chronic damage and future graft loss," Applied Biochemistry and Biotechnology, vol. 190 no. 2, pp. 660-673, DOI: 10.1007/s12010-019-03082-2, 2020.

[16] J. M. Lorenzen, C. Schauerte, M. Kölling, A. Hübner, M. Knapp, H. Haller, T. Thum, "Long noncoding RNAs in urine are detectable and may enable early detection of acute T cell–mediated rejection of renal allografts," Clinical Chemistry, vol. 61 no. 12, pp. 1505-1514, DOI: 10.1373/clinchem.2015.243600, 2015.

[17] M. K. Atianand, D. R. Caffrey, K. A. Fitzgerald, "Immunobiology of long noncoding RNAs," Annual Review of Immunology, vol. 35 no. 1, pp. 177-198, DOI: 10.1146/annurev-immunol-041015-055459, 2017.

[18] Q. Zhuang, Q. Liu, S. J. Divito, Q. Zeng, K. M. Yatim, A. D. Hughes, D. M. Rojas-Canales, A. Nakao, W. J. Shufesky, A. L. Williams, R. Humar, R. A. Hoffman, W. D. Shlomchik, M. H. Oberbarnscheidt, F. G. Lakkis, A. E. Morelli, "Graft-infiltrating host dendritic cells play a key role in organ transplant rejection," Nature Communications, vol. 7 no. 1,DOI: 10.1038/ncomms12623, 2016.

[19] S. Yazdani, J. Callemeyn, S. Gazut, E. Lerut, H. de Loor, M. Wevers, L. Heylen, C. Saison, A. Koenig, O. Thaunat, L. Thorrez, D. Kuypers, B. Sprangers, L. H. Noël, L. van Lommel, F. Schuit, M. Essig, W. Gwinner, D. Anglicheau, P. Marquet, M. Naesens, "Natural killer cell infiltration is discriminative for antibody-mediated rejection and predicts outcome after kidney transplantation," Kidney International, vol. 95 no. 1, pp. 188-198, DOI: 10.1016/j.kint.2018.08.027, 2019.

[20] M. A. Lim, J. Kohli, R. D. Bloom, "Immunosuppression for kidney transplantation: where are we now and where are we going?," Transplantation Reviews, vol. 31 no. 1, pp. 10-17, DOI: 10.1016/j.trre.2016.10.006, 2017.

[21] S. Roedder, N. Kimura, H. Okamura, S. C. Hsieh, Y. Gong, M. M. Sarwal, "Significance and suppression of redundant IL17 responses in acute allograft rejection by bioinformatics based drug repositioning of fenofibrate," PLoS One, vol. 8 no. 2,DOI: 10.1371/journal.pone.0056657, 2013.

[22] T. K. Sigdel, Y. Gao, J. He, A. Wang, C. D. Nicora, T. L. Fillmore, T. Shi, B. J. Webb-Robertson, R. D. Smith, W. J. Qian, O. Salvatierra, D. G. Camp, M. M. Sarwal, "Mining the human urine proteome for monitoring renal transplant injury," Kidney International, vol. 89 no. 6, pp. 1244-1252, DOI: 10.1016/j.kint.2015.12.049, 2016.

[23] P. Khatri, S. Roedder, N. Kimura, K. de Vusser, A. A. Morgan, Y. Gong, M. P. Fischbein, R. C. Robbins, M. Naesens, A. J. Butte, M. M. Sarwal, "A common rejection module (CRM) for acute rejection across multiple organs identifies novel therapeutics for organ transplantation," Journal of Experimental Medicine, vol. 210 no. 11, pp. 2205-2221, DOI: 10.1084/jem.20122709, 2013.

[24] B. D. Modena, S. M. Kurian, L. W. Gaber, J. Waalen, A. I. Su, T. Gelbart, T. S. Mondala, S. R. Head, S. Papp, R. Heilman, J. J. Friedewald, S. M. Flechner, C. L. Marsh, R. S. Sung, H. Shidban, L. Chan, M. M. Abecassis, D. R. Salomon, "Gene expression in biopsies of acute rejection and interstitial fibrosis/tubular atrophy reveals highly shared mechanisms that correlate with worse long-term outcomes," American Journal of Transplantation, vol. 16 no. 7, pp. 1982-1998, DOI: 10.1111/ajt.13728, 2016.

[25] G. Einecke, J. Reeve, B. Sis, M. Mengel, L. Hidalgo, K. S. Famulski, A. Matas, B. Kasiske, B. Kaplan, P. F. Halloran, "A molecular classifier for predicting future graft loss in late kidney transplant biopsies," The Journal of Clinical Investigation, vol. 120 no. 6, pp. 1862-1872, DOI: 10.1172/JCI41789, 2010.

[26] I. Diboun, L. Wernisch, C. A. Orengo, M. Koltzenburg, "Microarray analysis after RNA amplification can detect pronounced differences in gene expression using limma," BMC Genomics, vol. 7 no. 1,DOI: 10.1186/1471-2164-7-252, 2006.

[27] B. Chen, M. S. Khodadoust, C. L. Liu, A. M. Newman, A. A. Alizadeh, "Profiling tumor infiltrating immune cells with CIBERSORT," Cancer Systems Biology, pp. 243-259, DOI: 10.1007/978-1-4939-7493-1_12, 2018.

[28] A. M. Newman, C. L. Liu, M. R. Green, A. J. Gentles, W. Feng, Y. Xu, C. D. Hoang, M. Diehn, A. A. Alizadeh, "Robust enumeration of cell subsets from tissue expression profiles," Nature Methods, vol. 12 no. 5, pp. 453-457, DOI: 10.1038/nmeth.3337, 2015.

[29] A. Subramanian, P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, A. Paulovich, S. L. Pomeroy, T. R. Golub, E. S. Lander, J. P. Mesirov, "Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles," Proceedings of the National Academy of Sciences, vol. 102 no. 43, pp. 15545-15550, DOI: 10.1073/pnas.0506580102, 2005.

[30] A. Ramasamy, A. Mondry, C. C. Holmes, D. G. Altman, "Key issues in conducting a meta-analysis of gene expression microarray datasets," PLoS Medicine, vol. 5 no. 9,DOI: 10.1371/journal.pmed.0050184, 2008.

[31] Y.-Z. Ge, T. Xu, W.-J. Cao, R. Wu, W. T. Yao, C. C. Zhou, M. Wang, L. W. Xu, T. Z. Lu, Y. C. Zhao, Z. K. Hu, Z. L. Xu, X. B. Yang, L. H. Zhou, S. L. Zhang, B. Zhong, Z. Xu, W. C. Li, J. G. Zhu, R. P. Jia, "A molecular signature of two long non-coding RNAs in peripheral blood predicts acute renal allograft rejection," Cellular Physiology and Biochemistry, vol. 44 no. 3, pp. 1213-1223, DOI: 10.1159/000485451, 2017.

[32] X. Fan, U. M. Ashraf, C. A. Drummond, H. Shi, X. Zhang, S. Kumarasamy, J. Tian, "Characterization of a long non-coding RNA, the antisense RNA of Na/K-ATPase α 1 in human kidney cells," International Journal of Molecular Sciences, vol. 19 no. 7,DOI: 10.3390/ijms19072123, 2018.

[33] X. Fan, Regulation of Na/K-ATPase and Its Role in Cardiac Disease, 2018.

[34] S. Xu, J. Sui, S. Yang, Y. Liu, Y. Wang, G. Liang, "Integrative analysis of competing endogenous RNA network focusing on long noncoding RNA associated with progression of cutaneous melanoma," Cancer Medicine, vol. 7 no. 4, pp. 1019-1029, DOI: 10.1002/cam4.1315, 2018.

[35] J. Gong, S. Jin, X. Pan, G. Wang, L. Ye, H. Tao, H. Wen, Y. Liu, Q. Xie, "Identification of long non-coding RNAs for predicting prognosis among patients with thymoma," Clinical Laboratory, vol. 64 no. 7, pp. 1193-1198, DOI: 10.7754/Clin.Lab.2018.180136, 2018.

[36] A. D. Kirk, S. Ibrahim, D. V. Dawson, F. Sanfilippo, O. J. Finn, "Characterization of T cells expressing the γ / δ antigen receptor in human renal allografts," Human Immunology, vol. 36 no. 1, pp. 11-19, DOI: 10.1016/0198-8859(93)90003-J, 1993.

[37] W. Moolenaar, J. Bruijn, E. Schrama, S. Ferrone, M. R. Daha, A. H. Zwinderman, P. J. Hoedemaeker, L. van Es, F. van der Woude, "T-cell receptors and ICAM-1 expression in renal allografts during rejection," Transplant International, vol. 4 no. 3, pp. 140-145, DOI: 10.1007/BF00335334, 1991.

[38] Y. Li, Z. Huang, R. Yan, M. Liu, Y. Bai, G. Liang, X. Zhang, X. Hu, J. Chen, C. Huang, B. Liu, G. Luo, J. Wu, W. He, "V γ 4 γδ T cells provide an early source of IL-17A and accelerate skin graft rejection," Journal of Investigative Dermatology, vol. 137 no. 12, pp. 2513-2522, DOI: 10.1016/j.jid.2017.03.043, 2017.

[39] Q. Wu, P. K. Gupta, H. Suzuki, S. R. Wagner, C. Zhang, O. W. Cummings, L. Fan, M. H. Kaplan, D. S. Wilkes, R. A. Shilling, "CD4 T cells but not Th17 cells are required for mouse lung transplant obliterative bronchiolitis," American Journal of Transplantation, vol. 15 no. 7, pp. 1793-1804, DOI: 10.1111/ajt.13215, 2015.

[40] L. C. Sullivan, E. M. Shaw, S. Stankovic, G. I. Snell, A. G. Brooks, G. P. Westall, "The complex existence of γδ T cells following transplantation: the good, the bad and the simply confusing," Clinical & Translational Immunology, vol. 8 no. 9, article e1078,DOI: 10.1002/cti2.1078, 2019.

[41] M. R. Weir, M. Hall-Craggs, S. Y. Shen, J. N. Posner, S. V. Alongi, F. J. Dagher, J. H. Sadler, "The prognostic value of the eosinophil in acute renal allograft rejection," Transplantation, vol. 41 no. 6, pp. 709-712, DOI: 10.1097/00007890-198606000-00008, 1986.

[42] W. Hongwei, R. S. Nanra, A. Stein, L. Avis, A. Price, A. D. Hibberd, "Eosinophils in acute renal allograft rejection," Transplant Immunology, vol. 2 no. 1, pp. 41-46, DOI: 10.1016/0966-3274(94)90076-0, 1994.

[43] Z. Zhao, L. Wang, C. Yang, T. Zhao, L. Li, L. Hu, D. Wu, R. Rong, M. Xu, T. Zhu, "Soluble FGL2 induced by tumor necrosis factor- α and interferon- γ in CD4+ T cells through MAPK pathway in human renal allograft acute rejection," Journal of Surgical Research, vol. 184 no. 2, pp. 1114-1122, DOI: 10.1016/j.jss.2013.04.011, 2013.

[44] P. Grimbert, O. Thaunat, "mTOR inhibitors and risk of chronic antibody-mediated rejection after kidney transplantation: where are we now?," Transplant International, vol. 30 no. 7, pp. 647-657, DOI: 10.1111/tri.12975, 2017.

[45] M.-L. Alegre, D. R. Goldstein, A. S. Chong, "Toll-like receptor signaling in transplantation," Current Opinion in Organ Transplantation, vol. 13 no. 4, pp. 358-365, DOI: 10.1097/MOT.0b013e3283061149, 2008.

[46] R. M. Gangemi, A. Daga, L. Muzio, D. Marubbi, S. Cocozza, M. Perera, S. Verardo, D. Bordo, F. Griffero, M. C. Capra, A. Mallamaci, G. Corte, "Effects of Emx2 inactivation on the gene expression profile of neural precursors," European Journal of Neuroscience, vol. 23 no. 2, pp. 325-334, DOI: 10.1111/j.1460-9568.2005.04559.x, 2006.

[47] G. Spigoni, C. Gedressi, A. Mallamaci, "Regulation of Emx2 expression by antisense transcripts in murine cortico-cerebral precursors," PLoS One, vol. 5 no. 1, article e8658,DOI: 10.1371/journal.pone.0008658, 2010.

[48] Y. Gu, C. Feng, T. Liu, B. Zhang, L. Yang, "The downregulation of lncRNA EMX2OS might independently predict shorter recurrence-free survival of classical papillary thyroid cancer," PLoS One, vol. 13 no. 12, article e0209338,DOI: 10.1371/journal.pone.0209338, 2018.

[49] Z. Tang, G. Wei, L. Zhang, Z. Xu, "Signature microRNAs and long noncoding RNAs in laryngeal cancer recurrence identified using a competing endogenous RNA network," Molecular Medicine Reports, vol. 19 no. 6, pp. 4806-4818, DOI: 10.3892/mmr.2019.10143, 2019.

[50] J. Massagué, "TGF- β signal transduction," Annual Review of Biochemistry, vol. 67 no. 1, pp. 753-791, DOI: 10.1146/annurev.biochem.67.1.753, 1998.

[51] D. Kim, W. K. Lee, S. Jeong, M. Y. Seol, H. Kim, K. S. Kim, E. J. Lee, J. Lee, Y. S. Jo, "Upregulation of long noncoding RNA LOC100507661 promotes tumor aggressiveness in thyroid cancer," Molecular and Cellular Endocrinology, vol. 431, pp. 36-45, DOI: 10.1016/j.mce.2016.05.002, 2016.

[52] D. D. Esposti, H. Hernandez-Vargas, C. Voegele, N. Fernandez-Jimenez, N. Forey, B. Bancel, F. L. Calvez-Kelm, J. McKay, P. Merle, Z. Herceg, "Identification of novel long non-coding RNAs deregulated in hepatocellular carcinoma using RNA-sequencing," Oncotarget, vol. 7 no. 22, pp. 31862-31877, DOI: 10.18632/oncotarget.7364, 2016.

[53] A. V. Villarino, Y. Kanno, J. J. O'Shea, "Mechanisms and consequences of Jak–STAT signaling in the immune system," Nature Immunology, vol. 18 no. 4, pp. 374-384, DOI: 10.1038/ni.3691, 2017.

[54] M. Fukata, A. S. Vamadevan, M. T. Abreu, "Toll-like receptors (TLRs) and Nod-like receptors (NLRs) in inflammatory disorders," Seminars in Immunology, vol. 21 no. 4, pp. 242-253, DOI: 10.1016/j.smim.2009.06.005, 2009.

[55] B. J. Chen, F. L. Byrne, K. Takenaka, S. C. Modesitt, E. M. Olzomer, J. D. Mills, R. Farrell, K. L. Hoehn, M. Janitz, "Transcriptome landscape of long intergenic non-coding RNAs in endometrial cancer," Gynecologic Oncology, vol. 147 no. 3, pp. 654-662, DOI: 10.1016/j.ygyno.2017.10.006, 2017.

[56] C. Li, H. Zheng, W. Hou, H. Bao, J. Xiong, W. Che, Y. Gu, H. Sun, P. Liang, "Long non-coding RNA linc00645 promotes TGF- β -induced epithelial–mesenchymal transition by regulating miR-205-3p-ZEB1 axis in glioma," Cell Death & Disease, vol. 10 no. 10,DOI: 10.1038/s41419-019-1948-8, 2019.

[57] B. J. Nankivell, S. I. Alexander, "Rejection of the kidney allograft," New England Journal of Medicine, vol. 363 no. 15, pp. 1451-1462, DOI: 10.1056/NEJMra0902927, 2010.

[58] H. C. Probst, K. McCoy, T. Okazaki, T. Honjo, M. van den Broek, "Resting dendritic cells induce peripheral CD8+ T cell tolerance through PD-1 and CTLA-4," Nature Immunology, vol. 6 no. 3, pp. 280-286, DOI: 10.1038/ni1165, 2005.

[59] H. C. Probst, J. Lagnel, G. Kollias, M. van den Broek, "Inducible transgenic mice reveal resting dendritic cells as potent inducers of CD8+ T cell tolerance," Immunity, vol. 18 no. 5, pp. 713-720, DOI: 10.1016/S1074-7613(03)00120-1, 2003.