Blood Gene Expression Predicts Bronchiolitis Obliterans Syndrome

Bronchiolitis obliterans syndrome (BOS), the main manifestation of chronic lung allograft dysfunction, leads to poor long-term survival after lung transplantation. Identifying predictors of BOS is essential to prevent the progression of dysfunction before irreversible damage occurs. By using a large set of 107 samples from lung recipients, we performed microarray gene expression profiling of whole blood to identify early biomarkers of BOS, including samples from 49 patients with stable function for at least 3 years, 32 samples collected at least 6 months before BOS diagnosis (prediction group), and 26 samples at or after BOS diagnosis (diagnosis group). An independent set from 25 lung recipients was used for validation by quantitative PCR (13 stables, 11 in the prediction group, and 8 in the diagnosis group). We identified 50 transcripts differentially expressed between stable and BOS recipients. Three genes, namely POU class 2 associating factor 1 (POU2AF1), T-cell leukemia/lymphoma protein 1A (TCL1A), and B cell lymphocyte kinase, were validated as predictive biomarkers of BOS more than 6 months before diagnosis, with areas under the curve of 0.83, 0.77, and 0.78 respectively. These genes allow stratification based on BOS risk (log-rank test p < 0.01) and are not associated with time posttransplantation. This is the first published large-scale gene expression analysis of blood after lung transplantation. The three-gene blood signature could provide clinicians with new tools to improve follow-up and adapt treatment of patients likely to develop BOS.

inTrODUcTiOn Chronic lung allograft dysfunction is the main limitation of longterm survival after lung transplantation. It manifests mainly by abnormal remodeling of the small airways, resulting in progressive airflow obstruction called bronchiolitis obliterans syndrome (BOS) (1)(2)(3). The prevalence of BOS is around 35% at 5 years. Its late diagnosis, based on the irreversible decline of lung function, attests to irreversible and advanced degradation of the allograft. Prognosis is poor, with a 4-year median survival after disease diagnosis (4). Thus, there is an urgent need to identify the predictors of BOS, which would allow proactive and targeted strategies to slow disease progression before irreversible degradation of the allograft occurs.
Bronchiolitis obliterans syndrome is likely to arise from repeated injuries from both alloimmune and non-alloimmune mechanisms, generating fibrosis and airway obstruction (5). Tracking these inflammation and fibrotic processes has been used to identify early signs of the disease, and bronchoalveolar lavage neutrophilia, levels of regulatory T cells, chemokines/ cytokines, or matrix metalloproteases have been proposed as early biomarkers of BOS (6)(7)(8)(9)(10). More recently, expression profiling of lung biopsies pinpointed fibrosis-associated genes for the diagnosis or prediction of BOS (11). Yet, these invasive lung-centered approaches remain hampered by the accessibility to such samples and are therefore limited for routine monitoring of lung transplant recipients (LTRs). In blood, circulating fibrocytes and cytokine concentration have been proposed as predictors of BOS (12)(13)(14)(15)(16). However, these studies used a limited number of patients and have yet to be confirmed by follow-up studies. Consequently, none of these attempts have demonstrated sufficient feasibility and robustness to achieve clinical acceptance and potential routine use in the future.
Large-scale gene expression profiling of peripheral blood represents a promising tool to identify transcriptomic markers associated with the natural history of an allograft (17). The technique has been successfully tested for monitoring kidney (18), heart (19), and liver (20) transplant recipients. For the first time, we used this non-invasive approach to identify blood biomarkers of BOS in a large set of samples from LTRs. By using an independent set of patients, we were able to validate a transcriptomic signature that predicts the occurrence of BOS more than 6 months before the clinical manifestations.

Patients
Lung transplant recipients were recruited from September 2009 to October 2013 within the multicentre Cohort of Lung Transplantation (COLT) cohort (NCT00980967). The local ethical committee (Comité de Protection des Personnes Ouest 1-Tours, 2009-A00036-51) approved the study, and all participants provided written informed consent. Inclusion criterion was at least 3 years of follow-up unless the diagnosis of BOS was made before ( Figure S1 in Supplementary Material).
The eligible patients (n = 688) were phenotyped by a blind adjudication committee as described previously (21,22), based on pulmonary function tests and chest imaging according to ISHLT/ ERS/ATS guidelines (3,23). Assessed pulmonary tests were performed before the transplantation, the day of the transplantation, 1 and 6 months after the transplantation and every 6 months thereafter up to 3 years posttransplantation. We excluded 265 patients because of other phenotypes or confounding factors, and 338 stable patients and 85 BOS patients were identified. Stable patients displayed no signs of chronic dysfunction for at least 3 years after lung transplantation. Stable and BOS patients were then further selected to constitute homogenous groups based on sample availability, absence of concurrent infection or acute rejection within 1 month before or after blood collection, and the quality of RNA (RNA integrity number ≥6.5). Eighty-nine patients (49 STA and 40 BOS) were included in the identification set and 25 in the validation set (13 STA and 12 BOS).

rna isolation
Peripheral blood samples were collected in PAXgene tubes (PreAnalytix, Qiagen, Hilden, Germany) and stored at −80°C. Total RNA was extracted using the PAXgene blood RNA system kit (Qiagen) with an on-column DNase digestion protocol according to the manufacturer's instructions. Quantity and quality of total RNA were determined using a 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA), and RNA samples with an RNA integrity number above 6.5 were selected for further analyses.

gene expression Microarray analysis
Total RNA (100 ng) were labeled using the Two-Color Agilent Low Input Quick Amp Labeling Kit and hybridized on SurePrint G3 Human Gene Expression v3 8 × 60K Microarrays following the manufacturer's instructions (Agilent Technologies). Data extraction of median feature intensity was performed with Feature Extraction software v10.7 (Agilent Technologies). To remove signal intensity bias between each array, median feature intensities were normalized with the lowess (locally weighted scatter plot smoothing) method, then spots for which half the samples exhibited a signal less than the mean of all median signals were removed. Correction between two microarray hybridization batches was performed on the 28,867 remaining spots with the Combat algorithm (24) available through the R package sva (25). These normalized microarray data were deposited in the Gene Expression Ominbus database (accession number GSE94557). Mean expression levels for the spots targeting the same genes were assessed, resulting in 16,128 unique genes. In the BOS PRED group, the closest time point from transplantation was selected in cases of several samples per patient, so that no patient duplicate was included. Genes with low variation (i.e., variance <0.2), and thus considered as invariants, were excluded, resulting in 6,581 analyzed genes. For the identification of differential genes, linear modeling with empirical Bayes statistical procedure was performed, comparing the STA group and each group of interest, using the limma package in R. Genes with p < 5% and fold change (FC) >1.5 were considered as differentially expressed. The biological significance of selected genes was assessed using GOminer software. Only gene ontology (GO) categories enriched with a false discovery rate (FDR) <10% and with at least five represented genes were selected. The cell type source of differential genes was evaluated using the web tool Enrichr (26). Relative estimation of cell abundances using gene expression data was performed through the CIBERSORT software (27), and gene set enrichment analysis (GSEA) was performed to highlight relevant biological processes using curated gene sets from the Molecular Signatures Database (1,000 permutations and FDR threshold of 25%) (28).

gene expression Data sets
Gene expression values for the three genes of interest (BLK, POU2AF1, and TCL1A) from two public microarray data sets were collected from the Gene Expression Ominbus database: GSE38267 (29) and GSE28042 (30), corresponding to two studies analyzing blood gene expression in non-transplanted patients with terminal respiratory failure.

statistical analysis
Non-parametric Kruskal-Wallis tests with Dunn's ad hoc pairwise comparisons, Mann-Whitney tests, ROC curves, log-ranked survival analyses, and Fisher's exact test for categorical variables were performed using GraphPad Prism v. 6 (GraphPad Software, La Jolla, CA, USA). Time-to-event analysis was performed using Cox proportional analysis between gene expression and time to BOS with the R survival package (Coxph function, R software v. 3.3.2).

resUlTs lung Transplant recipients
Lung transplant recipients included in this study were recruited from the multicentre COLT cohort, allowing for longitudinal follow-up and 6-month interval biocollection from transplantation. On the basis of this longitudinal follow-up, we retrospectively defined two classes of BOS samples depending on the time between blood collection and BOS diagnosis (defined as the time point with a decline of ≥20% in FEV1 from baseline; Figure 1A). Blood samples collected at least 6 months before BOS diagnosis were included in the prediction class (PRED), and blood samples collected at the time or after BOS diagnosis (up to 13 months after initial diagnosis) were included in the diagnosis class (DIAG). For the group of patients with stable graft function (STA), blood was collected 6 and 12 months after transplantation, and a comparison of these two time points was performed to exclude irrelevant genes possibly altered during this interval after transplantation ( Figure 1B). No patient duplicate was included within any classes. Clinical parameters of patients in the identification and validation sets are presented in Table 1. LTR groups were homogeneous regarding age, sex, BMI, type of transplantation, induction treatment, smoking status, and infection and rejection events. A significant difference in azithromycin exposure was observed in the identification set between the DIAG and STA and PRED (76.9 vs. 22.4 and 37.5%, p < 0.0001). Predicted FEV1 was significantly lower in the DIAG in both cohorts (p < 0.05). Time of blood collection differed between the DIAG and PRED (p < 0.001).
identification of gene signatures associated With BOs  BOS groups, excluding differential transcripts between STA at the 2 time points (Figure 2; Table S1 in Supplementary Material). We identified 33 transcripts differentially expressed between STA and DIAG. While no GO term was significantly enriched, GSEA analysis highlighted 14 enriched gene sets (FDR < 25% and nominal enrichment <−2; Table S2 in Supplementary Material) among which 3 were related to B cell biology (31,32) and one to upregulated genes in lung tissue of smokers with chronic obstructive pulmonary disease (COPD) (33). Twenty-nine transcripts were differentially expressed between STA and PRED. According to the GO analysis, five genes participated to the enrichment of the immune system process Based on unsupervised hierarchical clustering of expressed genes, these genes clustered with known B cell-related genes such as MS4A1 (membrane-spanning 4-domains, subfamily A, member 14 also called CD20 molecule), BANK1 (B cell scaffold protein with ankyrin repeats 1), and CD40, reinforcing the potential association of B cell-related genes with prediction of BOS. Estimation of memory and plasma B cells relative abundances with gene expression using CIBERSORT did not highlighted significant alteration, while naive CD4 + T cells abundance was decreased in PRED compared to STA ( Figure S2 in Supplementary Material). It is noteworthy that 14 transcripts, including TCL1A, immunoglobulin lambda-like polypeptide 5 gene (IGLL5), and various immunoglobulin lambda and kappa light chain variable regions, were associated with both the DIAG and the PRED.
While these results highlight genes differentially expressed between STA and each of the two BOS groups, BOS is a timedependent phenomenon after lung transplantation. Thus, we reanalyzed differential genes from the STA and PRED comparison using a time-dependent analysis, the Cox proportional hazards test. According to the univariate analysis, all 29 differentially expressed genes displayed a significant association with BOS occurrence with time (p value of Wald test <0.05; Table S1 in Supplementary Material).

Validation of POU2AF1, TCL1A, and BLK as Predictive Biomarkers of BOs
Among the differentially expressed genes between STA and PRED, five were selected on the basis of their p values, FC magnitude, and expression level in microarrays. These five genes were measured by qPCR in an independent set of 25 patients (13 STA and 12 BOS), with 11 and 8 samples in the PRED and DIAG, respectively ( Table 2). Downregulation of 3 genes, POU2AF1 (p = 0.0065), TCL1A (p = 0.0257), and BLK (p = 0.0221) in the PRED was validated by qPCR ( Figure 3A). By contrast, the downregulation of CD19 and IGLL5 were not confirmed. Expression of POU2AF1, TCL1A, and BLK was constant in the STA between 6 and 12 months posttransplantation ( Figure S3 in Supplementary Material) but was not significantly downregulated at the time of transplantation between STA and PRED (n = 6 and 9, respectively; Figure S4 in Supplementary Material). For diagnostic purposes, we confirmed the downregulation of TCL1A (p = 0.0265; Figure 3B).
Because POU2AF1, TCL1A, and BLK were differentially expressed in the PRED, that is, more than 6 months before the clinical diagnosis of BOS, we evaluated the predictive performance of these three markers using ROC curve analysis. POU2AF1 (AUC = 0.832, 95% CI = 0.638-1.026), TCL1A (AUC = 0.773, 95% CI = 0.553-0.993), and BLK (AUC = 0.780, 95% CI = 0.569-0.991) discriminated well between STA and BOS patients (Figure 4A). Global performances of the prediction presented in Figure 4 show an accuracy higher than 80% for the three markers. Due to the high level of correlation between the three genes, performance of the prediction was not improved by their combination (Figure S5 in Supplementary Material).  In the survival analysis, these three genes were highly associated with BOS occurrence with time in the discovery set (BLK, p = 0.0013; POU2AF1, p = 0.0028; and TCL1A, p = 0.0031; Table S1 in Supplementary Material). The association between BOS-free survival and POU2AF1, TCL1A, and BLK was assessed through Kaplan-Meier analyses. As shown in Figure 4C, expression levels of POU2AF1, TCL1A, and BLK under 0.45, 0.34, and 0.505, respectively (corresponding to best expression thresholds according to ROC curves) were associated with significant reduction of BOS-free survival after lung transplantation (p < 0.01). In this validation cohort, 3 of 8 recipients developing BOS had donor-specific antibody (DSA) at blood collection vs. none of the 13 STA (Fisher p = 0.031). Using univariate logistic analysis, DSA was not significantly associated with BOS.

POU2AF1, TCL1A, and BLK are Downregulated in Blood of Patients with end-stage chronic respiratory Diseases
Given the relationship between POU2AF1, TCL1A, and BLK expression and the development of BOS, we investigated the expression of these three genes in public data sets from blood of patients with other causes of respiratory failure (29,30). The three genes were downregulated in the blood of non-transplanted patients with end-stage chronic respiratory diseases (cystic fibrosis, idiopathic pulmonary fibrosis, or pulmonary hypertension) compared to healthy volunteers ( Figure 5).

DiscUssiOn
Bronchiolitis obliterans syndrome is a major limitation of long-term survival after lung transplantation. Several previous attempts to identify early predictors of BOS remain limited by the absence of validation studies (6,(8)(9)(10)(12)(13)(14). Furthermore, they generally relied on invasive procedures incompatible with a routine clinical monitoring. We opted for a non-invasive largescale molecular profiling approach to identify predictors of BOS in whole blood from 89 LTRs. A three-gene molecular signature differentiating BOS and STA was validated in an independent set of 25 patients. To the best of our knowledge, this is the first published study combining two independent cohorts for the identification and validation of predictors of BOS. We showed that 29 and 33 transcripts were differentially expressed between STA and DIAG and STA and PRED. Twelve transcripts were common between both comparisons. This suggests a progression of the molecular profile during the development of the disease. Comparison of blood collected 6 and 12 months after transplantation for the STA group ruled out an early time-dependent alteration of gene expression. POU2AF1, TCL1A, and BLK were associated with the occurrence of BOS with time in a survival analysis and were validated as predictors of BOS more than 6 months before its diagnosis. The validation cohort used in our study revealed the strength of these biomarkers for BOS prediction. They displayed similar performances to predict BOS according to their correlation of expression, with AUC values reaching 0.83, 0.77, and 0.78, respectively, and exceeding the performance of biomarkers already published (14,34,35).
Recent literature highlights the correlation between de novo DSAs with the occurrence of chronic dysfunction (36)(37)(38)(39)(40). While in our validation cohort DSA was not significantly associated with BOS, an independent study is required to decipher whether these genes are associated with de novo DSA and BOS appearance and which parameter is the most predictive. Furthermore, our investigation focused on the BOS subtype, and further work will have to be conducted to identify predictors of the restrictive syndrome. Besides, stringent patient selection criteria were applied for a clear discrimination between groups, and all selected patients were not included in the analysis because of feasibility matters. This reduction in patient number may limit the generalizability of our results, whereas enough power was reached to highlight significantly statistical differences for the three genes. Although our predictions were confirmed on an independent set of patients, they will have to be validated on an external cohort as well.
Our primary objective was to identify a molecular biomarker for clinical monitoring in whole blood, which limits the interpretation of mechanistic investigations. Nevertheless, B cell signatures and the presence of B cells with suppressive properties are evidenced in other situations of transplantation and particularly in tolerance in kidney transplantation (41)(42)(43) and/or associated with long-term graft survival (44,45). In a longitudinal study in lung transplantation, BOS appearance is associated with a lower level of CD24 hi CD38 hi IgD hi IgM hi transitional B cells with a regulatory phenotype 18 months after lung transplantation compared to STA (46). We also observed that 12 and 13 transcripts related to immunoglobulins are downexpressed in BOS (PRED and DIAG, respectively) compared to STA, while gene expression of BAFF (B cell activating factor coded by TNFSF13B) and BAFF-R (BAFF receptor coded by TNFRSF13C), required for B cell survival and maturation, were not different between the three situations (data not shown). These data suggest a dysregulation of B cells in the BOS situation.
In addition, several genes merit further exploration to better understand the underlying mechanism behind BOS. By using GO analysis, we evidenced differentially expressed genes between STA and PRED associated with immune system process, such as CD19, HLA-DQA1, and HLA-DQA2. Interestingly, the gene enrichment and unsupervised hierarchical clustering analyses reinforce this association of the B cell cluster with the development of BOS. POU2AF1 is a B cell transcriptional coactivator involved in B cell development and function (47). Yet, POU2AF1 is expressed in T cells as well, and recent evidence revealed its role in Th1-Th2 polarization (48) and in the mounting of T-celldependent B cell responses (49). BLK is a member of the Src family of tyrosine kinases and encodes a non-receptor protein tyrosine kinase involved in the regulation of B cell receptor signaling (50). We found a downregulation of TCL1A in the BOS before disease onset. TCL1A is notably expressed by B and T lymphocytes, where it promotes cell proliferation and survival (51). Interestingly, in renal transplantation, TCL1A is downregulated at the time of acute allograft rejection (52), whereas it is overexpressed in operationally tolerant patients, an ideal situation where recipients are off of immunosuppression with a functioning allograft (53). The exact contribution of these genes in the development of BOS remains to be investigated. Some polymorphisms have been described for BLK, POU2AF1, and TCL1A. BLK polymorphism is a risk factor for developing autoimmune diseases (54,55). Genetic variations leading to reduced BLK expression are associated with several autoimmune diseases via lowering the threshold for B cell activation (54). A TCL1A single-nucleotide polymorphism is associated with downstream expression of cytokines and chemokines and the nuclear factor-κB transcriptional activity resulting in the modulation of inflammation and immune response (56). POU2AF1 polymorphism is associated with susceptibility to lymphoma (57), and mice deficient in POU2AF1 exhibit reduced numbers of mature B cells and defective immune responses to antigens (58,59). The fact that none of our patients evidenced clinical symptoms or biological associated modifications (such as decrease of total B cells), in addition with the absence of differential expression of these three genes at the time of transplantation do not suggest existence of such polymorphisms in our study. Intriguingly, however, the analysis of gene expression data from two independent studies (29,30) revealed that BLK, TCL1A, and POU2AF1 are downregulated in blood from patients with end-stage chronic respiratory diseases. GSEA analysis also highlighted the enrichment of a gene set related to COPD in lung from smokers in DIAG compared to STA, including the SPIB gene (33). Altogether, these results suggest common mechanisms between BOS and other chronic respiratory diseases, as previously noted (60), and support the need to further decipher the roles of POU2AF1, BLK, and TCL1A in the development of lung pathologies.
In conclusion, by using non-invasive whole blood profiling, we identified and validated POU2AF1, TCL1A, and BLK as three predictive biomarkers of BOS, more than 6 months before diagnosis. These genes allow stratification based on the BOS risk and could be easily monitored to provide clinicians with new tools to improve follow-up and adapt treatment of patient likely to develop BOS before clinical manifestations and allograft damage arise.

acKnOWleDgMenTs
We are most grateful to the GenoBiRD Core Facility for its technical support and to the biological resource center for biobanking CHU Nantes, Hôtel Dieu, Centre de ressources biologiques (CRB), Nantes, F-44093, France (BRIF: BB-0033-00040).

FUnDing
This research program was supported by Vaincre la Mucoviscidose, Association Gregory Lemarchal, Agence de Biomédecine, INSERM, Région Pays de La Loire, Institut de Recherche en Santé Respiratoire des Pays de la Loire, Fonds de Recherche en Santé Respiratoire, Fondation du Souffle, and Fondation pour la Recherche Médicale. This work was realized in the context of the IHU-Cesti project, the DHU Oncogreffe and the LabEX IGO thanks to French government financial support managed by the National Research Agency via the "Investment into the Future" program (ANR-10-IBHU-005 and ANR-11-LABX-0016-01).