Gene Expression–Based Predictive Models of Graft Versus Host Disease–Associated Dry Eye

Abstract

Purpose: To develop a predictive model based on inflammatory gene mRNA expression in conjunctival cells of graft versus host disease (GvHD)–associated dry eye (DE) patients, as well as to find meaningful correlations between gene signals and clinical signs.

Methods: Twenty GvHD-DE patients and 14 healthy controls were recruited. Patients discontinued medications for 1 week before examination. Dry eye–related symptoms and signs were recorded, and conjunctival epithelial cells were collected by impression cytology after spending 20 minutes under standard conditions within a Controlled Environmental Research Laboratory. Gene expression of inflammatory molecules was determined by polymerase chain reaction, and the results were correlated with clinical signs. Shrinkage discriminant analysis, support vector machine, and k-nearest neighbor classifier methods were used to develop predictive models that were validated considering accuracy, calibration, and discriminant capability.

Results: Out of the 84 genes analyzed, 34 showed significant differences in expression. IL-6, IL-9, CCL24, CCL18, IL-10, IFN-γ, and CCL2 were highly increased (>6-fold); 26 genes were moderately upregulated (2- to 6-fold), whereas EGFR was downregulated (2.63 fold) in GvHD-DE samples. A panel based on EGFR, IL-6, IL-9, and NAMPT had an area under the receiver operating characteristic curve of 0.994, a sensitivity of 100%, and a specificity of 92.9%. EGFR expression correlated negatively with ocular surface damage markers, while IL-6, IL-9, and NAMPT correlated positively with these tests.

Conclusions: EGFR, IL-6, IL-9, and NAMPT have the greatest potential as diagnostic biomarkers, with excellent sensitivity, specificity, and clinical relevance to the ocular surface status of GvHD.

Hematopoietic stem cell transplantation (HSCT) has rapidly evolved over the past decades, becoming a cornerstone in the management of malignant and nonmalignant hematologic disorders.¹ However, patients subjected to allogenic HSCT can develop graft versus host disease (GvHD), which constitutes the main cause of morbidity and mortality in these patients.² Graft versus host disease, which has both acute and chronic forms,³ is an immune-mediated destruction of host tissues by immunocompetent cells from the donor.⁴ 

Ophthalmic findings are more frequent in chronic GvHD, affecting up to 60% to 90% of patients.^5–7 Although the posterior segment can be affected, usually in the form of infectious complications from the concomitant immunosuppressive state, anterior segment involvement is much more prevalent. For the anterior segment GvHD, dry eye (DE) disease (keratoconjunctivitis sicca) is the most common complication. 

Graft versus host disease–associated DE (GvHD-DE) can be mild or moderate, but it is often the most severe type of DE, similar to and often worse than Sjögren-related DE. It can lead to serious abnormalities of the ocular surface, affecting the patient's quality of life and eventually leading to permanent visual loss. Graft versus host disease–DE is mainly due to aqueous tear deficiency, and histological analysis shows Sjögren syndrome-like inflammatory destruction of the conjunctiva and lacrimal gland with fibrosis, resulting in tear production deficiency.⁸ In 2005, the Schirmer test was proposed as the paramount diagnostic test for GvHD-DE patients³ despite its well-known lack of sensitivity and specificity.⁹ The general lack of adequately validated objective tests to diagnose DE was recognized in the 2007 International Dry Eye Workshop,¹⁰ and controversy remains as to whether or not GvHD-DE diagnosis should rely on tear dynamic testing, clinical exams, or new technologies.¹¹ As DE has been defined as an inflammatory disease of the lacrimal functional unit (LFU),¹² inflammatory mediators are being researched as potential biomarkers in DE in general, and in GvHD-DE in particular. In fact, different study groups have already described changes in cytokines, chemokines, growth factors, and their receptor levels in tear and/or conjunctival epithelial cells in DE patients.^13–20 Thus, there is a growing interest in the identification of the molecules that could be specifically involved in the pathogenesis of GvHD-DE. 

A major barrier to GvHD-DE research and treatment is that diagnosis and prognosis rely almost entirely on the presence of clinical symptoms and unreliable tests, like the Schirmer test. No currently used laboratory tests can predict the risk of developing the disease, responsiveness to treatment, or patient survival. Thus, the importance of biomarkers is crucial. A good diagnostic test for GvHD must be able to distinguish patients with GvHD-DE from those without GvHD-DE. Further, any biomarker should specifically and sensitively reflect the disease state and should be noninvasive, inexpensive, and standardized.²¹ Several clinical trials related to GvHD-DE patients are currently in progress, evaluating the efficacy of different anti-inflammatory drugs. Thus, the progress in GvHD-DE biomarker identification and validation is important, as clinical trial design will begin incorporating these biomarkers. 

Moreover, there is a lack of studies focused on LFU changes in these GvHD-DE patients, and to our knowledge, only one study has been published describing gene expression analysis in the conjunctival cells of these patients.²² In this context, the use of predictive models could help in understanding not only the importance of single gene expression, but also the importance of the coordinated expression of groups of genes involved in GvHD-DE. Such predictive models rely on the identification of gene signatures that consist of a group of genes expressed in association with a biological state of interest. A properly constituted predictive model could increase the diagnostic and prognostic power of gene expression data in clinical applications. These kinds of models seek to predict a classification variable, for example, susceptibility to GvHD-DE, based on known gene expression responses. Additionally, model validation, which evaluates how close the statistical model is to reality, is necessary and gives important information about the reliability of the model and the means by which it is derived. 

The first aim of the present study was to analyze the expression of a panel of inflammatory gene mediators in conjunctival epithelial cells collected from GvHD-DE patients. The second aim was to develop and validate different prediction models based on the expression of these genes. 

Materials and Methods

Patients and Healthy Controls

This study was approved by the Institute of Applied Ophthalmobiology (IOBA) institutional review board and the University of Valladolid Clinical Hospital Ethics Committee, Valladolid, Spain, and followed the tenets of the Declaration of Helsinki. All enrolled patients and subjects were informed of the nature of this study, and written consent was obtained from all of them. 

Graft versus host disease–DE patients were selected from those referred to IOBA by the Hematology Department of the University of Salamanca Clinical Hospital, Salamanca, Spain. To ensure as much as possible that the results would not be influenced by the effects of topical medications and/or any other condition, inclusion and exclusion criteria were established to ensure that patients were stable. Consequently, the most severely affected patients were excluded. For the same reason, patients included in the study were asked to discontinue their topical therapies for 7 days before the samples were taken. Only artificial tears and lubricants were allowed. Thus, inclusion criteria for patients were as follows: (1) abnormal results of at least three DE diagnostic tests that included the ocular surface disease index (OSDI) score > 12 points, fluorescein tear breakup time (T-BUT) < 7 seconds, corneal fluorescein staining and conjunctival Lissamine green staining > 1 (Oxford scale), and Schirmer test without topical anesthesia ≤ 5 mm in 5 minutes; (2) feasibility of discontinuing topical anti-inflammatory medications for 1 week, as judged by the attending ophthalmologist (MC); and (3) patients had to be systemically stable with GvHD under control with no relapse of the patient's primary malignancies and no secondary infections as judged by the referral physician (DC). Exclusion criteria included any ocular active disease other than DE, contact lens wear, any ophthalmic surgery in the past 6 months, any eye drop other than artificial tears and lubricants, and any systemic medication that was not continuous or constant in treatment and dosage for at least 3 months prior to this study. 

Healthy volunteers, similar in age and sex to the study cases, were selected as the control group and were examined to make certain their ocular status was within normal limits. Inclusion criteria for this group included (1) absence of ocular surface-related symptoms (OSDI score ≤ 12 points) and (2) normal limits in at least two of the following four ocular surface tests: fluorescein T-BUT ≥ 7 seconds, corneal fluorescein staining and conjunctival Lissamine green staining (Oxford Scale) ≤ grade 1, and Schirmer test without anesthesia > 5 mm in 5 minutes. Exclusion criteria for this control group included any present or previous history of ophthalmic or systemic disease, any ophthalmic surgery, under any medication, pregnancy, or current contact lens wear. 

Clinical Evaluation and Sample Collection

All individuals were evaluated after 20-minute exposure in a controlled environmental chamber (VISIÓN I+D, SL; Valladolid, Spain) located in the Controlled Environmental Research Laboratory (CERLab) at IOBA, University of Valladolid, Valladolid, Spain.^23,24 The temperature was maintained at 23°C and the relative humidity (RH) at 45%, which corresponds to a comfortable indoor temperature and the average RH in Valladolid. 

Clinical evaluation was always performed by the same two clinicians (LC and VM). Each one always evaluated the same tests in the following sequence: (1) OSDI questionnaire,²⁵ (2) conjunctival hyperemia, (3) tear stability, (4) corneal integrity, (5) tear production, and (6) conjunctival epithelial cell collection. 

For the OSDI questionnaire, the eye with the higher score was selected as the most symptomatic eye. A random data table was used in those cases in which both eyes were equally symptomatic and for control subjects. The most symptomatic eye was used in the clinical evaluation and for conjunctival cell sampling. Conjunctival hyperemia was evaluated under a slit lamp and based on the Nathan-Efron scale (0–4).²⁶ Tear stability was evaluated by measuring T-BUT. Five microliters of 2% sodium fluorescein was gently applied into the outer third of the inferior fornix with a micropipette. The time between the last of three blinks and the appearance of the first dry spot was measured three times, and the mean value was recorded. 

Corneal integrity was evaluated using a slit lamp mounted with a cobalt blue filter (Topcon Corp., Tokyo, Japan) and a yellow Wratten no. 12 filter (Eastman Kodak, Rochester, NY, USA). The evaluation was done 2 minutes after instillation of 5 μL 2% sodium fluorescein. For the evaluation of conjunctival integrity, Lissamine green strips (GreenGlo; HUB Pharmaceuticals, LLC, Rancho Cucamonga, CA, USA) were wet with 25 μL sodium chloride and then gently applied into the inferior fornix. Corneal and conjunctival staining were scored using the Oxford Scale (score 0–5).²⁷ Tear production was assessed by Schirmer test without topical anesthesia. A Schirmer sterile strip (Tearflo; HUB Pharmaceuticals, LLC) was placed in the lateral canthus of the inferior lid margin. The length of wetting was measured in millimeters after 5 minutes with eyes closed. 

Conjunctival epithelial cell collection was performed by conjunctival impression cytology (CIC) under topical anesthesia. Two halves of polyethersulfone filter (Supor 200, pore size: 0.20 μm, diameter: 13 mm; Gelman Laboratory, Ann Arbor, MI, USA) were gently applied on the upper and upper-temporal bulbar conjunctiva for 10 seconds under moderate pressure. Then the filters were suspended in 1 mL lysis buffer (RLT Buffer; Qiagen, Hilden, Germany) containing 1% 2-mercaptoethanol (Merck, KGaA, Darmstadt, Germany) and stored at −80°C for subsequent isolation of total RNA. 

Gene Expression Analysis in Conjunctival Epithelial Cells by Quantitative Real-Time PCR (qPCR)

RNA Isolation and cDNA Synthesis.

For RNA isolation, total RNA was extracted from each CIC sample by a commercial kit (RNeasy Micro Kit, Qiagen) and treated with RNase-free DNase following the manufacturer's instructions. The first strand of cDNA was synthesized from the total extracted RNA using the commercial kit HT RT2 First Strand (Qiagen). One hundred nanograms of cDNA from each sample was synthesized. 

Real-Time PCR Array Analysis.

Gene analysis using a real-time PCR (RT-PCR) array was performed according to the manufacturer's recommendations using a customized array Profiler (SuperArray SABioscience, Qiagen, Izasa, S.A., Barcelona, Spain) that used SYBR Green I dye detection. We determined the expression of 84 human genes (Table 1), along with the expression of five housekeeping genes: actin beta, beta-2-microglobulin, glyceraldehyde-3-phosphate dehydrogenase, hypoxanthine phosphoribosyltransferase 1, and large ribosomal protein P0. Additionally, the array included one control to monitor genomic DNA contamination, three reverse transcription controls for first strand synthesis, and three positive PCR controls for real-time PCR efficiency. Genes included in the customized array were selected based in part on a previous study (unpublished results) in which differential gene expression was analyzed in a pool of patient samples versus a pool of control samples. Other genes were included because they were known to be important in either GvHD or DE.^19,21 

Table 1

 

Panel of Inflammatory Genes Included in Customized SuperArray Profiler

 

The assay was done following the manufacturer instructions. Briefly, 102 μL (100 ng) of the diluted first strand cDNA synthesis reaction was mixed with 1275 μL 2× SuperArray PCR Master Mix (SABioscience, Qiagen, Izasa, S.A.) and 1173 μL double-distilled H₂O. This mixture and the template cocktail (25 μL each) were added to each well of the PCR array. Real-time PCR (7500 Real Time PCR System; Applied Biosystems, Madrid, Spain) was then performed as follows: 10 minutes at 95°C, 40 cycles of 15 seconds at 95°C, and 1 minute at 60°C. A melting curve program generated the dissociation curve for each well in the entire plate to verify the identity of each gene amplification product. 

Samples were individually analyzed. However, in some GvHD-DE patients, the amount of RNA obtained from CIC was not enough for individual analysis (<100 ng); therefore, three of these samples consisted of a pool of patients. The first pooled sample was from patients 2, 3, and 4, and the second was from patients 5 and 6. The third pooled sample was from patients 17, 18, 19, and 20. 

For data analysis, the Ct method, where the target fold = 2^−ΔΔCt (User Bulletin, No. 2, P/N 4303859; Applied Biosystems), was used. The results were reported as a fold upregulation when the fold change for patients was greater than 1 compared to the healthy controls. If the fold change was less than 1, the negative inverse of the result was reported as a fold downregulation. 

Statistical Analysis

Statistical analysis was performed by a licensed statistician (IF) using the R software (R Foundation for Statistical Computing, Vienna, Austria) and packages from the Bioconductor project.²⁸ Data were determined as means ± standard deviations, unless otherwise specified. To compare samples of independent groups, the nonparametric Mann-Whitney U test was used. To assess the association between qualitative variables, Fisher's exact test was used. A P value ≤ 0.05 was considered statistically significant. 

The expression stability of two candidate reference genes, actin beta and glyceraldehyde-3-phosphate dehydrogenase, was estimated using the geNorm algorithm²⁹ with the ReadqPCR and NormqPCR packages.³⁰ The analysis of selected genes was conducted using moderated t-test,³¹ implemented in the HTqPCR package.³² To quantify correlations between gene expression levels and clinical parameters, the Spearman ranked correlation coefficient (ρ) was used, along with a P value that contrasts the null hypothesis ρ = 0; that is, there is no correlation between tear levels and clinical parameters. Moreover, a bootstrap confidence interval for this coefficient was built using 5000 bootstrap samples. 

Construction of the Gene Expression-Based Predictive Models

We conducted gene signature identification using three different classification models. The selection of relevant genes was based on the minimum redundancy-maximum relevance (mRMR) algorithm³³ through repetitions of the leave-one-out cross-validation (LOOCV) process (Fig. 1). We used incremental feature selection to determine the optimal size of the signature. The LOOCV process was repeated 84 times, the total number of candidate genes. The LOOCV-based error rate for each m-gene signature was estimated, where an m-gene signature is a group of “m” genes within which the combined expression could be used to identify GvHD-DE patients. The optimal size of the signature was defined as the size of the set with the best prediction accuracy and smallest gene number. 

Figure 1

 

Gene signature identification process. Gene signature (m) identification was conducted using three different classification models. The selection of relevant genes was based on the minimum redundancy-maximum relevance (mRMR) algorithm, through repetitions of the leave-one-out cross-validation (LOOCV) process.

Two-class discriminative models were built using three techniques: (1) shrinkage discriminant analysis (SDA), (2) support vector machine (SVM), and (3) k-nearest neighbor classifier (kNN). The classical discriminant analysis involves the determination of a linear combination of predictors of gene expression, termed the “discriminant function,” that maximizes the differences between classes. In general, this method has suboptimal performance when the sample size is small and the number of predictors is large because the commonly used estimators for the variances can become unstable and reduce the classification accuracy. One solution is to make use of shrinkage-based variance estimators. In this study, we used SDA based on James-Stein-type shrinkage estimators.³⁴ 

Support vector machine aims at searching for a hyperplane that separates the two classes with the largest margin or distance between the hyperplane and the point closest to it.³⁵ A key factor in SVM is to use kernels to construct the decision boundaries. We used linear kernels with the e1071 package³⁶ to fit this model. 

The kNN is a distance-based approach for classification.³⁷ To classify a new observation with the kNN algorithm, the setting k, as a predefined parameter, finds the k-closest observations and classifies the new observation by majority vote. Thus, it chooses the class that is most common among the k neighbors. In this work, we used the Euclidean distance, and the optimal k was decided based on the classification error estimated by a LOOCV procedure. The class package³⁸ was used to fit this model. 

Using the Brier score as a global measure, we determined the accuracy, calibration, and discrimination capability of the models.³⁹ The Brier score is based upon individual differences between predicted risks in terms of likelihood and observed final outcomes. The Brier score ranges from 0 for a perfect degree of agreement to 1 for the worst possible degree of agreement. 

We used two methods to evaluate the calibration of the model, measured as the degree of agreement between the predicted and observed values. With the calibration in the large (CL) method, a perfectly calibrated model will have a CL value of 0 and a calibration slope of 1. For this method, the scores are not bounded, so the model will be badly calibrated if the values depart from the optimal values. We also used the Hosmer-Lemeshow test, which is significant for badly calibrated models. 

We used the receiver operating characteristic (ROC) curve analysis to assess the discriminant ability of the fitted models. They were evaluated and compared according to the area under the ROC curve (AUC). In addition, the sensitivity and specificity for each ROC curve were obtained by setting an optimal threshold using the pROC package.⁴⁰ 

Results

Clinical Evaluation

A total of 20 GvHD-DE patients (13 males, 7 females; 56.8 ± 11.6 years of age; range, 34–72 years) and 14 healthy volunteers (8 males, 6 females; 51.9 ± 14.7 years of age; range, 30–75 years) were recruited. There were no significant differences in sex (P = 0.916) or age (P = 0.3533) between the groups. Clinical history, including previous diagnosis for HSCT, time to GvHD-DE diagnosis, and systemic and topical therapies, is shown in Table 2. 

Table 2

 

Clinical Data of Graft Versus Host Disease Dry Eye (GvHD-DE) Patients

 

The GvHD-DE patients had significantly more frequent and intense ocular surface symptoms than controls, as reflected by a higher OSDI questionnaire score (Table 3). Ocular surface examination, including conjunctival hyperemia, T-BUT, corneal fluorescein staining, Lissamine green conjunctival staining, and Schirmer test, revealed that GvHD-DE patients always had significantly altered clinical tests compared with controls. However, the values were not extremely altered, reflecting the fact that patients were selected from among those who could stop topical medications for 1 week and who therefore had only a moderate level of inflammation. 

Table 3

 

Ocular Examination Parameters

 

Conjunctival Epithelial Cell Gene Expression Analysis

Out of the 84 genes studied (Table 1), 33 were significantly upregulated in GvHD-DE patients compared to controls (Table 4). Of those, the genes most upregulated were IL-6, IL-9, CCL24, CCL18, IL-10, IFN-γ, and CCL2, with a fold increase in value between 17.50 and 6.26. The expression of IL-17E, IFNA2, IL-12A, TNF, BMP2, CCL3, CXCL12, CSF2, IL-18, IL-17B, CCL11, IL-4, CSF3, IL-9R, CXCL3, TGFB3, IL-12B, IL-1B, IL-5, TGFB2, TNFSF11, IL-22R, CXCL13, VEGFA, IL-15, and NAMPT genes was upregulated 2- to 6-fold. Additionally, the expression of IL-5RA, IL-2, CCL4, CXCL1, and TNFSF10 was upregulated with borderline statistical significance (0.05 < P value ≤ 0.08). Finally, EGFR gene was the only gene significantly downregulated (−2.59 fold) in GvHD-DE patients compared to healthy controls. 

Table 4

 

Fold Up- or Downregulation of Gene Expression in Conjunctival Epithelial Cells of Graft Versus Host Disease–Dry Eye Patients Compared to Healthy Controls

 

Predictive Model Results

Three different classification methods, SDA, SVM, and kNN, were used to analyze a panel of molecules expressed in conjunctival epithelial cells that could help diagnose and monitor GvHD-DE in a better way than any single-gene evaluation (Fig. 2). Based on the output of the mRMR algorithm, we tested the accuracy of each of the three models with m genes as predictors. The m-parameter was determined by incremental feature selection, starting with m = 1 and adding one gene at time until m = 84. The LOOCV-based performance for each of the m-gene models was estimated and used to determine the optimal m in each case. The accuracy reached the maximum value when four genes, EGFR, IL-6, IL-9, and NAMPT, were included in the SVM-based model (Fig. 2). For the kNN- and SDA-based models, a panel of five and eight genes, respectively, was optimal. For the kNN model, IFN-γ was included along with the four genes used in the SVM model. For the SDA model, IFN-γ, TGFB2, CCL3, and VEGFA were included with the four genes from the SVM model. 

Figure 2

 

Leave-one-out cross-validation method for estimating accuracy of the SDA, SVM, and kNN predictive models for GvHD-DE. The number of genes used for classification is represented on the x-axis, and the y-axis corresponds to the estimated accuracies evaluated by LOOCV of each of the three fitted models. The LOOCV-based error rate for each m-gene signature was estimated. The optimal size of the signature was defined as the size of the set with best prediction accuracy and smallest number of genes. LOOCV, leave-one-out cross-validation; SDA, linear discriminant analysis, SVM, support vector machine; kNN, k-nearest neighbor classifier.

Internal Validation

We used the LOOCV procedure to validate our results (Table 5). Three aspects were evaluated in each case: accuracy, calibration, and discrimination ability. All models were accurate, with Brier scores close to 0, and none of them presented serious problems of calibration. With respect to the discriminative performance, the SVM-based model was somewhat better than the others, with an AUC of 0.995, a sensitivity of 100%, and a specificity of 92.9%. 

Table 5

 

Internal Validation Output

 

Taking into account all of these results, we determined that the model obtained by the SVM method, based on EGFR, IL-6, IL-9, and NAMPT gene expression levels, produced the best model to predict the occurrence of GvHD-DE. Thus, this panel consisting of four biomarkers effectively discriminated between patients with and without GvHD-DE. 

Correlations Between Selected Gene Expression and Clinical Data

There were significant correlations between the expression of EGFR, IL-6, IL-9, and NAMPT and the clinical tests (Fig. 3). Considering the entire sample of patients and controls, the expression of EGFR was negatively correlated with the OSDI questionnaire score, conjunctival hyperemia, fluorescein corneal staining, and Lissamine green conjunctival staining. In contrast, EGFR expression was positively correlated with T-BUT and Schirmer test results. IL-6 and IL-9 were positively correlated with hyperemia and fluorescein staining, and negatively correlated with T-BUT and Schirmer test results. Additionally, IL-6 was positively correlated with OSDI questionnaire score. Finally, NAMPT gene expression was positively correlated with OSDI, hyperemia, and corneal staining while inversely correlated with T-BUT and Lissamine green staining. 

Figure 3

 

Correlations of EGFR, IL-6, IL-9, and NAMPT gene expression with clinical parameters. Correlations were determined for the entire study population, both patients and controls together. ▴, GvHD-DE patients; ○, control group subjects; values of hyperemia, corneal fluorescein staining, and conjunctival Lissamine green staining were jittered to reduce overplotting. Dashed lines represent the 95% confidence interval; rho, Spearman ranked correlation coefficient.

Discussion

In the present study, we analyzed the expression of a panel of inflammatory mediator genes in conjunctival epithelial cells of patients with GvHD-DE. Inflammation-related genes were selected as both pathologies, DE and GvHD, have been demonstrated to have an inflammatory etiology, and proteins encoded by them have shown an important role in mobilization and recruitment of immune cells in affected tissues.^12,21 We then developed and validated different prediction models based on the expression of these genes. 

Clinical evaluation of the enrolled patients revealed that all of them suffered a moderate-to-severe DE as shown by the clinical test scores. The GvHD-DE patients had significantly more DE symptoms, lower Schirmer test scores, and considerably decreased T-BUT values compared to healthy controls. Ocular surface integrity, as evaluated by fluorescein and Lissamine green vital staining, was significantly altered compared to controls; however, it was not so severely damaged that the patients could not stop their anti-inflammatory topical medications for a week as required by the inclusion criteria. 

Of the 84 genes studied in conjunctival epithelial cells, 34 had a statistically significant different expression in the GvHD-DE patients compared to the controls. Those that demonstrated the most significant differences between groups were IL-6, IL-9, CCL24, CCL18, IL-10, IFN-γ, CCL2, and EGFR. Some of the proteins encoded by these genes have already been related to GvHD or to DE.^18,41,42 When we look at the changes in mRNA expression, these results give us information about which inflammatory pathways may be active in the disease. 

When developing a disease predictive model, not only should differential gene expression be taken in account, but also the biological significance of the results must be assessed. Multiple biomarker-based predictive models may increase sensitivity, specificity, and predictive or diagnostic performance compared to the use of single-molecule models. Multiple biomarker models are a frequently applied technique in the study of systemic GvHD. Paczesny et al.⁴³ determined that a panel of four biomarkers in serum (IL-2Rα, TNFR1, IL-8, and HGF) could discriminate patients with acute GvHD from controls. Most recently, a three-biomarker panel in serum, consisting of IL-2Ra, TNFR1, and elafin, was tested prospectively at different times, and it showed reasonably good specificity (75%) and fair sensitivity (57%) for the development of acute GvHD.⁴⁴ Furthermore, other research groups have established different biomarker-based models in GvHD.^45,46 We intended to apply this technology to the accurate diagnosis and prediction of ocular complications of GvHD, and for that we used minimally invasive techniques like CIC for ocular sample collection and biomarker determination. To find the best model, we selected genes for which the activity or expression was not redundant or a reflection of the increase or activation of another gene. Additionally, the three predictive models that we used, SDA, SVM, and kNN, are appropriate for small sample size and large number of feature datasets, as in our case.⁴⁷ According to the results obtained by these methods, the panel of four genes, EGFR, IL-6, IL-9, and NAMPT, produced the best model to predict the occurrence of GvHD-DE. 

EGFR, which encodes the protein epidermal growth factor receptor, was the only mRNA that was downregulated in GvHD-DE patients. The EGFR protein is one of the most targeted receptors in the field of oncology, where inhibitors have demonstrated clinical success in specific cancers.⁴⁸ Moreover, EGFR is expressed on corneal epithelial cell surfaces, where it is considered to play an important role in cell proliferation during corneal wound healing.⁴⁹ Its activation appears to be effective in reducing the time of re-epithelialization and has potential uses in severe DE.⁵⁰ In fact, soluble EGFR is significantly increased in tears of DE patients compared to healthy controls.⁵¹ To our knowledge, it has not been described in ocular samples of GvHD-DE patients or elsewhere in systemic GvHD patients. 

In contrast to EGFR, IL-6 is a well-known molecule in GvHD. It is a multifunctional cytokine involved in the regulation of the immune response, hematopoiesis, and inflammation, and there is evidence suggesting that IL-6 plays a pivotal role during the transition from innate to acquired immunity.⁵² Furthermore, under proinflammatory stimuli, conjunctival epithelial cells are an important source of this cytokine.⁵³ Regulatory T cell function is inhibited by IL-6 produced by activated dendritic cells,⁵⁴ and anti-IL-6 receptor antibody administration in an acute GvHD mouse model induced T-regulatory cells⁵⁵ and decreased GvHD-related mortality.⁵⁶ Therefore, IL-6 plays a critical role in GvHD pathophysiology.⁵⁷ Many studies have shown increased levels of IL-6 in patients with severe GvHD.^58,59 The median level of tear IL-6 from ocular GvHD patients is significantly elevated compared to controls,⁶⁰ and additionally, expression is increased in CIC samples of DE patients.^15,16 Based on these data, IL-6 could be a sensitive disease biomarker, but it seems to lack specificity, as it is increased in diverse post-HSCT complications, such as infections, mucositis, and venous occlusive disease.⁶¹ 

IL-9 has been most frequently associated with allergic inflammation and immunity to extracellular parasites. Additionally, it plays a role in Th1/Th17-mediated inflammation and in T-regulatory response.⁶² IL-9 was first recognized as a Th2-related cytokine, but it is now clear that it is produced by various CD4+ T-cell subsets, including the newly defined Th9 cells.⁶³ In the present study, IL-9 seems to play an important role in the pathogenesis of GvHD-DE, but to our knowledge, serum or tear levels have not been studied yet. Tawara et al.⁶⁴ concluded that deficiency of the four classical Th2 cytokines (IL-4, IL-5, IL-9, and IL-13) in mice enhances T-cell proliferative responses and aggravates GvHD. The importance of this finding has to be confirmed, and further investigation is needed. If substantiated, this molecule could play an important role in the pathogenesis of GvHD-DE because an allergic disease was definitely excluded in all of our patients. IL-9 might be a good biomarker for this disease, as its specificity is higher than that of other molecules studied. 

The fourth gene in the predictive model, NAMPT (nicotinamide phosphoribosyltransferase, also called visfatin), produces an enzyme that catalyzes nicotinamide adenine dinucleotide (NAD) synthesis. NAMPT is also a proinflammatory cytokine that promotes B cell maturation and inhibits neutrophil apoptosis.⁶⁵ It has been implicated in the susceptibility to and pathogenesis of a number of human pathologic conditions, including malignancies, obesity, and diabetes. Interestingly, NAMPT has been considered a novel marker of chronic inflammation, and the mRNA is upregulated in mononuclear cells from patients suffering chronic inflammatory diseases like psoriasis, rheumatoid arthritis, or inflammatory bowel disease. The high level of mRNA expression can be used to discriminate these patients from healthy controls.⁶⁶ Previous studies have determined that NAMPT inhibition has catastrophic consequences in activated T cells, leading to functional impairment and cell death. This approach could be a novel strategy to selectively eliminate activated T cells and block immune reactions such as those occurring in GvHD.⁶⁷ To our knowledge, NAMPT levels have not been described in GvHD patients in general or in ocular surface-related pathology. Based on our results, this still poorly understood cytokine looks promising, and more studies should be done in this direction. 

In this study, we correlated differential gene expression of the four selected genes with the clinical tests performed in our sample of GvHD-DE patients and healthy subjects. We found significant correlations for EGFR, IL-6, IL-9, and NAMPT with multiple tests that are clinically meaningful. For instance, lower expression of EGFR was associated with lower levels of T-BUT and Schirmer test scores, denoting worse tear stability and production, respectively. Conversely, lower EGFR expression was also associated with higher OSDI questionnaire scores, hyperemia, and ocular surface staining, denoting a worsening in symptoms, redness, and ocular surface damage, respectively. Similarly, increased expression of the proinflammatory genes IL-6, IL-9, and NAMPT correlated in general with higher symptoms, more ocular surface damage, and lower tear production and stability. These meaningful clinical correlations strengthen the potential role of these four genes as biomarkers of disease. 

Gene signatures, consisting of groups of genes expressed in a pattern that is uniquely characteristic of a biological phenotype or medical condition, are widely used in oncology,⁶⁸ hematology, and also in chronic GvHD studies.⁶⁹ These signatures, in association with a very large set of functional association data, may help reveal other genes that are related to a set of input genes. The gene signature that we found for GvHD-DE is an attempt to apply this technology to ophthalmology in general and in particular to the ocular surface inflammatory diseases. The predictive model that we have developed goes beyond testing currently available in the clinical practice, mainly Schirmer test, which lacks the sensitivity and specificity required to make a correct diagnosis. If it works properly in an external validation setting, it will provide a more accurate and earlier diagnosis that will allow physicians to promptly address the ocular manifestations of these patients, avoiding sight-threatening and permanent complications. 

The present study was limited by the relatively small sample of GvHD-DE patients. However, the number of participating patients was not inordinately small considering the low (although increasing) frequency of ocular GvHD. It was also limited by the necessity to exclude patients who were either so severely ill that they were unable to attend our clinic or could not discontinue their topical anti-inflammatory medications. Additionally, it should be considered that in some GvHD-DE patients the amount of RNA obtained from CIC was not enough for individual analysis, and this amount might be lower in more severe patients of GvHD-DE, which could be a limitation of the test in clinical practice. 

In summary, we showed that the LFU of GvHD-DE patients presented alterations in the expression of several inflammatory molecule mRNAs in conjunctival epithelial cells that could help explain the pathophysiology of this ocular disorder. Furthermore, the biomarker-based models that we developed showed that the expression of four promising genes could become diagnostic biomarkers of ocular GvHD-DE. Future directions include evaluation of these biomarkers with a larger patient population in prospective studies.⁷⁰ The next (and necessary) step in our work will be performing an external validation of the selected gene-predictive model, including not only the new cohort of GvHD-DE patients, but also other groups of patients, such as GvHD patients without DE, GvHD-independent DE patients (Sjögren-related DE, for instance) and a group of patients with a different ocular surface inflammatory disease (like allergic keratoconjunctivitis, for instance). This work is warranted before using this predictive model. Such studies could facilitate the successful design of subsequent clinical trials. Furthermore, biomarkers may represent novel therapeutic targets that could be inhibited by specific drugs, thus enabling control of these processes and improving patient quality of life. 

Acknowledgments

Supported by a national grant from Ministry of Economy and Competitiveness, Madrid, Spain, SAF-2010 15631 (AES). The authors alone are responsible for the content and writing of the paper. 

Disclosure: L. Cocho, None; I. Fernández, None; M. Calonge, None; Allergan, Inc. (C); V. Martínez, None; M.J. González-García, None; D. Caballero, None; L. López-Corral, None; C. Garcia-Vazquez, None; L. Vazquez, None; M.E. Stern, Allergan, Inc. (E); A. Enríquez-de-Salamanca, None