Reduced Human Leukocyte Antigen (HLA) Protection in Gulf War Illness (GWI)

Background Gulf War Illness (GWI) is a disease of unknown etiology with symptoms suggesting the involvement of an immune process. Here we tested the hypothesis that Human Leukocyte Antigen (HLA) composition might differ between veterans with and without GWI. Methods We identified 144 unique alleles of Class I and II HLA genes in 82 veterans (66 with and 16 without GWI). We tested the hypothesis that a subset of HLA alleles may classify veterans in their respective group using a stepwise linear discriminant analysis. In addition, each participant rated symptom severity in 6 domains according to established GWI criteria, and an overall symptom severity was calculated. Findings We found 6 Class II alleles that classified participants 84.1% correctly (13/16 control and 56/66 GWI). The number of copies of the 6 alleles was significantly higher in the control group, suggesting a protective role. This was supported by a significant negative dependence of overall symptom severity on the number of allele copies, such that symptom severity was lower in participants with larger numbers of allele copies. Interpretation These results indicate a reduced HLA protection (i.e. genetic susceptibility) in veterans with GWI. Funding University of Minnesota and U.S. Department of Veterans Affairs.


Introduction
Shortly after the Gulf War (1990-91), veterans started to report a variety of health problems that began during, or soon after returning from, deployment, prompting investigation into the epidemiology and etiology of the complaints. Those investigations revealed that diffuse symptoms such as fatigue, musculoskeletal pain, mood and neurocognitive complaints, gastrointestinal problems, and rashes were most commonly reported. The constellation of symptoms, now commonly referred to as Gulf War Illness (GWI), has affected a substantial number of Gulf War veterans (Fukuda et al., 1998;Institute of Medicine, 2006, 2010Kang et al., 2009;Steele, 2000). Several population-based studies have demonstrated that these symptoms occur at significantly higher rates in deployed Gulf War veterans relative to their nondeployed peers and other veterans (Fukuda et al., 1998;Kelsall et al., 2009;Unwin et al., 1999), raising the issue about possible in-theater exposures and stress as contributing factors. However, these symptoms are also present in non-deployed military personnel (Steele, 2000), leading some to suspect other causes, including reactions to vaccine adjuvants (Israeli, 2012;Toubi, 2012). In summary, GWI is now a recognized constellation of symptoms of unclear etiology, also co-occurring with psychiatric disorders.
To date, the pathophysiology of GWI remains poorly understood. The overlap of GWI symptoms with symptoms of autoimmune disorders (e.g., chronic fatigue), together with evidence of abnormal immune activation following exercise challenge (Broderick et al., 2012(Broderick et al., , 2013, raise the possibility that GWI reflects an abnormal immune process (Broderick et al., 2012(Broderick et al., , 2013Hotopf et al., 2000;Israeli, 2012;Moss, 2013;Toubi, 2012). This probable involvement of altered immune mechanisms prompted our investigation of Human Leukocyte Antigen (HLA) in GWI.
HLA genes are located in the Major Histocompatibility Complex (MHC) of chromosome 6 and play a central role in immune recognition (Meuer et al., 1982). Most investigations of association of HLA to various diseases have focused on evaluating HLA allele frequencies in diseases of interest, as compared to the general, healthy population. Such studies have demonstrated HLA involvement with cancer, autoimmune, and infectious diseases (Rioux et al., 2009;Trowsdale and Knight, 2013). HLA Class I proteins (HLA-A, B, C) are expressed on all nucleated cells and present peptides from endogenous proteins to cytotoxic T lymphocytes engaged in immune surveillance. HLA Class II proteins (HLA-DRB1, DRB3/4/5, DQB1, DPB1) are expressed on antigen-presenting cells and present peptides derived from exogenous proteins to CD4+ helper T cells (Reche and Reinherz, 2003). A previous study of Gulf War syndrome in 27 veterans found that HLA DRB1*15 was more prevalent in cases than controls with an odds ratio of 1.66, although this association was not statistically significant (O'Bryan et al., 2003).
Here we tested the hypothesis that HLA may be a contributing factor in developing GWI by performing a stepwise linear discriminant analysis which assessed how well the number of copies of HLA alleles could classify participants in their respective group.

Study Participants
Veterans Affairs (VA) medical records were reviewed to identify potential participants. Participants were recruited if they had completed a Gulf War physical or if they had served during the Persian Gulf War. Exclusionary criteria included central nervous system disorders (e.g. Parkinson's disease, dementia, cerebral vascular accidents, a history of traumatic brain injury, etc.), lifetime psychotic or bipolar diagnoses, and current drug dependence. Veterans who might have had difficulty understanding the protocol were also excluded from recruitment. Participants provided written informed consent prior to participation and were compensated for their time. The study protocol was approved by the Institutional Review Board at the Minneapolis VA Health Care System.
We studied a total of 82 participants, 16 controls (15 men, 1 woman; age range 43-71 years; 54.9 ± 10.2 years, mean ± SD) and 66 GWI (64 men, 2 women; age range 39-76 years; 50.6 ± 7.9 years). The mean age did not differ significantly between the 2 groups (P = 0.13, t-test). Participants completed a symptom presence/severity questionnaire developed for use in Kansas Gulf War veterans (Steele, 2000) that evaluates a range of symptoms associated with GWI and permits determination of case status according to either the Centers for Disease Control and Prevention (CDC) criteria (Fukuda et al., 1998) or the Kansas GWI case definition (Steele, 2000). Participants meeting either set of Gulf War case definition criteria were included in the study. All participants had deployed during the Gulf War and were free of any autoimmune disease.
The questionnaire asks participants to indicate if they have had a persistent problem over the last 6 months with various symptoms from the following six domains: fatigue, pain, neurological-cognitivemood, skin, respiratory, and gastrointestinal. For each symptom rated as present, participants are asked to rate the severity of the symptom as mild, moderate, or severe, and to indicate whether the symptom first became problematic before, during or after deployment to the Gulf. Only symptoms that began during or after Gulf War service are counted toward diagnosis. There were six symptom domains: fatigue, pain, neurological-cognitive-mood, skin, gastrointestinal, and respiratory (Steele, 2000). Individual symptom severity was reported in a scale from 0 to 3. For each participant, an average score per domain was calculated, and a grand average across domains was used in the regression analysis below.

Data Analysis
Standard statistical tests were used to analyze the data using the IBM-SPSS statistical package (version 23) and ad hoc Fortran computer programs employing the International Mathematics and Statistics Library (IMSL; Rogue Wave Software, Louisville, CO, USA) statistical and mathematical libraries. We performed a stepwise linear discriminant analysis to identify those alleles that would correctly classify participants to their respective groups (control, GWI) using the SPSS statistical package above. In that analysis, group membership was the classification variable and the numbers of copies for each allele for each participant were the predictor variables (N = 144 alleles). This analysis identified 6 alleles that, as a group, classified correctly 84.1% of the participants (see below). We validated the robustness of this allele set by calculating the classification rate using the leave-one-out method in the original sample and by performing an extensive analysis based on the bootstrap (Efron and Tibshirani, 1993), as follows. Each participant was classified 10,000 times using bootstrap samples of participants and associated data from the 6 alleles above as predictors. Specifically, each bootstrap sample was generated by random sampling with replacement from each one of the 2 groups (control and GWI) to a total of N = 200 participant values per group (in separate additional analysis, bootstrap samples of N = 100 participants were used). We then performed a linear discriminant classification analysis to classify the given participant to one of the 2 groups, based on the 2 bootstrapped participant samples. Finally, we calculated the correct classification rate for each group. The quality and strength of classification was assessed using standard methods, including chi-square analysis, calculating confidence intervals on the sensitivity and specificity using the binomial theorem and Wilson's conservative test (Wilson, 1927), and calculating the receiver operating characteristic (ROC) curve and its associated statistics.
A different issue concerns the possible effect of unequal group sizes on the results of the discriminant analysis (Sanchez, 1974). Ideally, the group sizes should be equal, or as closely equal as possible, because, then, the null hypothesis for chance classification is the same or very similar across groups. To remedy this sample inequality, Sanchez (1974) proposed the following procedure aimed to perform the discriminant analysis using the same (or very similar) sample sizes. Let G1 and G2 be two groups to be discriminated with sample sizes m and km, respectively. The standard discriminant analysis would be somewhat hampered by the unequal sample sizes, km N m. To remedy this, Sanchez (1974) proposed that, instead of a single discriminant analysis, k such analyses be performed, where one group is always the same N(G1) = m, whereas the other group is a random subsample (without replacement) from G2, of size NðG2 the same to that of G1, with equal chance probabilities of 0.5. The final classification rate is the average of the k classification outcomes. We applied this approach to our data as follows. The sample size of the control group was 16, and that of the GWI was 66. Since 66 cannot be divided exactly by 16, we performed the following steps.
Step 1: The original control group (N = 16) was used always as one of the two groups to be classified.
Step 3: Four consecutive subsamples of n = 16 each were identified and four discriminant analyses performed between the original control group and each subsample above. This yielded four classification rates the average of which provided an overall classification rate against a chance probability of 0.5, since sample sizes were equal.
Step 4: This procedure left out two participants from the GWI group. To remedy this, we performed Steps 2 and 3 10,000 times, each with a different permutation. In this way, all GWI participants contributed to the analysis.
Step 5: The final classification rate was derived as the average of the 10,000 rates above and its value assessed against the expected chance rate of 0.5 (50%).
We assessed the direction of the effect of those alleles (protective vs. predisposing) by carrying out four additional analyses. First, we compared the frequencies of those alleles between the two groups using a paired t-test. Second, we calculated the odds ratio (ω) for each allele, took its natural logarithm, ln(ω), and tested the null hypothesis that its average does not differ significantly from zero, using a one-sample t-test. (If any cell in the 2 × 2 table had zero counts, the odds ratio was calculated after the constant 0.5 was added to all counts in the table, to avoid taking the logarithm of zero or dividing by zero.) Third, we performed a linear regression analysis where the overall symptom severity was the dependent variable and the number of copies of the 6 alleles was the independent variable; this analysis was also performed separately for each symptom domain. Significantly lower average allele frequencies in GWI, a significant negative ln(ω) and a significant negative slope (and correlation) in the regression analysis would argue for a collective protective effect of the 6 alleles. Finally, we compared the frequencies of those 6 alleles in our two groups (GWI, controls) to the frequencies that have been reported in the literature. Specifically  Rossman et al., 2003). Analyses and comparisons with our data were carried out separately for each one of the 3 databases above.

Results
A total of 144 unique HLA alleles (74 Class I and 70 Class II, Table 1) were identified from a total of 82 participants. The stepwise linear discriminant analysis identified 6 alleles that yielded 84.1% correct classification of the 82 subjects to their respective groups. These results were robust. The leave-one-out classification rate for the original sample was 79.1% and for the 2 bootstrap analyses N 80% (80.1% for bootstrap samples of N = 100 and 81.2% for N = 200). We also performed the same stepwise discriminant analysis in a gender-homogeneous sample for men only (N = 79). The correct classification rate was 86.1% and the leave-one-out classification rate was 81.0%. The overall classification rate obtained by applying the procedure proposed by Sanchez (1974) (see Methods) was 82.3% (range: 78.1%-85.9%, N = 10,000 permutations), which is way substantially and consistently above the expected chance rate of 50%.
The alleles and associated statistics are given in Table 2 and details about the goodness of classification based on the results using the full GWI sample are presented in Table 3. It can be seen that the discriminant classification was highly statistically significant and effective, and that the classification was correct at a high level (N80%), was very similar with respect to sensitivity (0.848; 84.8%) and specificity (0.812; 81.2%), exceeded chance even by the most conservative Wilson's test (Table 3E), and yielded a ROC curve (Fig. 1) that was highly significant and considerably above chance (Table 3F).
The analysis of frequencies of the 6 alleles (Table 4, A, B) pointed to a systematic protective effect in the control group and, by extension, lack of protection in the GWI group. This is evidenced by the fact that all allele frequencies were lower in the GWI group (Table 4A), as compared to the control group (t = − 5.789, DF = 5, P = 0.002, paired t-test) and that all the odds ratios were less than one, i.e. negative ln(ω): ln(ω) = − 1.792 ± 0.383 (mean ± SEM), t = − 4.671, DF = 5, P = 0.005, one-sample t-test against the null hypothesis that mean ln(ω) = 0). In addition, the percentage of participants with a given allele was systematically lower for all 6 alleles (Table 4B).
This collective protective effect of the 6 alleles above was further corroborated by the results of the linear regression of overall symptom severity against the number of allele copies (Fig. 2) which revealed a strong and highly significant negative relation (t = −4.148, DF = 80, P = 0.000083, R 2 = 0.177): GWI symptom severity ¼ 4:55−2:044 number of 6 allele copies We also performed a multiple linear regression analysis separately for each symptom domain. We found a highly significant effect of the number of copies of individual alleles on symptom severity of Pain (P = 0.01, R 2 = 0.199), Fatigue (P = 0.006, R 2 = 0.210), and Neurological-Cognitive-Mood (P = 0.004, R 2 = 0.225) symptoms; remarkably, all slopes (i.e. partial regression coefficients) of individual allele copies vs. symptom severity were negative in each regression model, indicating a consistent and robust effect. Finally, the regression analysis was not significant for skin (P = 0.911), gastrointestinal (P = 0.576), and respiratory (P = 0.598) symptoms.
Next, we analyzed our data with respect to the three databases (SF, MN, CDC) from the literature. Fig. 3 shows the mean frequencies across the 6 alleles as found in the 3 databases, our controls and the GWI group. It can be seen that the values for the 3 databases are similar, whereas the values for the GWI and control groups are lower and higher, respectively, than those of the databases. Next, we carried out a statistical analysis of these data by calculating the ln(ω) for each allele between each database and control, and each database and GWI populations. This yielded 6 alleles × 3 databases × 2 groups = 36 values. We assessed the main effects of Database and Group, and the Database × Group interactions by performing a univariate ANOVA. We found that the Group main effect was highly significant ( Fig. 4; P = 0.00019), whereas neither the Database main effect (P = 0.911; data not shown) nor the Database × Group interaction ( Fig. 5; P = 0.934) were statistically significant (F-test in the ANOVA). These results document the significant difference between the control (increased) and GWI (decreased) mean allele frequencies, relative to the 3 databases, while positing the similarity among the 3 databases and the similarity of the group effect across these databases. Finally, we searched the "Allele*Frequencies in Worldwide Populations" (http://www.allelefrequencies.net/) website for known haplotypes related to the 6 discriminating alleles, within the 3 databases that we used. Haplotype information was available only for the SF database, as shown in Table 5. It can be seen that no identified haplotype contains more than one of the 6 discriminating alleles, indicating their unique contribution.

Discussion
These results document for the first time significant differences in the HLA makeup between veterans with GWI and Gulf War era veterans without it. All the evidence from this study points to an enhanced vulnerability (or lack of protection) of the GWI veterans and, conversely, a reduced vulnerability (or additional protection) of the healthy veterans who served in the Gulf War but did not suffer from GWI. Collectively, our findings are in keeping with other evidence for an immune    dysfunction in GWI (Broderick et al., 2012(Broderick et al., , 2013Craddock et al., 2015;Hotopf et al., 2000;Israeli, 2012;Moss, 2013;Toubi, 2012;Whistler et al., 2009), in addition to other factors, including inflammatory components (Broderick et al., 2012(Broderick et al., , 2013Johnson et al., 2013), mitochondria dysfunction (Koslik et al., 2014) and genetic variants regarding butyrylcholinesterase enzyme activity (Steele et al., 2015). In fact, our findings provide a genetic susceptibility framework within which environmental triggers (e.g. vaccines, exposure to chemicals, stress, etc.) can be discussed, investigated and evaluated. All six HLA alleles that were singled out by the discriminant analysis to classify successfully (Fig. 1) control and GWI participants belonged in Class II and came from all three major genes (DQB1, DPB1, DRB1). Two alleles (DQB1*02:02 and DRB1*08:11) were absent from the GWI population, and the frequencies of the remaining four were all lower in the GWI group than in the control group (Table 4A). In addition, the percentage of participants with any of the six alleles was lower than in the controls (Table 4B). These results, together, document the lower frequency of occurrence of these alleles in the GWI group, as compared to the control group. The results of the regression analysis further documented the protective association of those alleles, given the significant negative slope in Fig. 2. In addition, this analysis further differentiated this effect among specific symptom domains, since it showed highly significant overall protective effects for Pain, Fatigue, Neurological-Cognitive-Mood domains. Remarkably, the effect of each allele on symptom severity in each one of these domains was protective, as evidenced by the universally negative partial regression coefficients in those three separate regression analyses. In contrast, the effect was not significant for the Skin, Gastrointestinal and Respiratory domains.
The comparison of the observed allele frequencies to those reported in published databases (Ovsyannikova et al., 2005;Rossman et al., 2003;Skibola et al., 2012) (Figs. 3-5) documented an additional finding, namely that both GWI and controls differed systematically and in opposite directions (higher for controls, lower for GWI,Figs. 4 and 5) with respect to the published allele frequencies. This suggests the following hypothesis regarding the role of the six discriminating alleles we identified. First, we assume that the population of GW-era veterans came from a larger population whose allele frequencies would be very similar to those reported in the three databases we used in this study. This is a reasonable assumption. Second, it is known that GW veterans were administered various vaccines, possibly together and/or multiply (Institute of Medicine, 2000;Petersdorf et al., 1996;Schwartz et al., 1997;Steele, 2000) and they were exposed to various chemical agents (Institute of Medicine, 2000;Steele et al., 2015). Third, a proportion of GW veterans developed GWI consisting of chronic symptoms affecting multiple organ systems (Fukuda et al., 1998;Institute of Medicine, 2006, 2010Kang et al., 2009;Kelsall et al., 2009;Steele, 2000), mostly affecting veterans who were deployed but also present in nondeployed or minimally exposed veterans (Steele, 2000). Based on those facts, we hypothesize that vaccinations and/or chemical exposures of GW era veterans served as environmental triggers ("hits") that contributed to developing GWI in genetically (HLA) "vulnerable"  veterans. Specifically, we propose that the presence of certain HLA alleles in higher frequencies conferred protection, whereas their relative paucity conferred vulnerability. This, in turn, resulted in the two distinct subpopulations of our study, namely control, deployed GW veterans with higher allele frequencies and absence of GWI on the one hand, and deployed GW veterans with lower frequencies and presence of GWI (Fig. 3). It should be noted that both vaccinations (Hotopf et al., 2000;Israeli, 2012;Toubi, 2012;Rook and Zumla, 1997) and chemical exposures (Moss, 2013) have been implicated previously as contributing factors for the development of GWI. The results of our study simply add a genetic susceptibility framework within which the effects of the factors above could be interpreted and investigated in future work. The nature of this postulated protection is likely to relate to autoimmune as well as inflammatory processes, since HLA has been implicated in both (Trowsdale and Knight, 2013;Johnson et al., 2013;Blackwell et al., 2009). In addition, HLA genetic underpinnings in immune responses to vaccines have been well established (Ovsyannikova et al., 2006;Ovsyannikova and Poland, 2011;Poland et al., 2007Poland et al., , 2008aPoland et al., , 2008b. Specific HLA makeup could thus manifest as autoimmune reactions, aberrant immune response to vaccinations, and/or increased susceptibility to infections (Nicolson, 1998); all three of them have been implicated in GWI.
A major limitation of our study is the small number of participants studied which, nevertheless, yielded systematic, significant and robust differences between the control and GWI groups but also between those two groups and published HLA allele databases. In this context, it is especially noteworthy that high correct classification rates were obtained in our extensive bootstrap analyses and permutation analyses (Sanchez, 1974) which documented the robustness of our results, i.e. their extension to larger sample populations. However, extensive bootstraps cannot substitute for validation in new, independent samples, which remains to be done. Finally, it is reasonable to assume that the six alleles identified in the present study are only part of a protective HLA makeup. We anticipate that a similar but large scale study will identify more such alleles and will provide a firm background to investigate the molecular mechanisms by which such protection/vulnerability may be mediated (Wahren-Herlenius and Dorner, 2013).

Contributions
LMJ and BEE contributed to data collection and clinical evaluation; APG and JJ contributed to data analysis; MYM and AG contributed to DNA extraction and arranging for HLA genotyping; APG, LMJ, AG wrote the paper. All authors contributed to extensive editing of the paper.

Declaration of interests
The authors do not report any financial disclosures or conflicts of interest.