Innate immunity in ocular Chlamydia trachomatis infection: contribution of IL8 and CSF2 gene variants to risk of trachomatous scarring in Gambians

Background Trachoma, a chronic keratoconjunctivitis caused by Chlamydia trachomatis, is the world's commonest infectious cause of blindness. Blindness is due to progressive scarring of the conjunctiva (trachomatous scarring) leading to in-turning of eyelashes (trichiasis) and corneal opacification. We evaluated the contribution of genetic variation across the chemokine and cytokine clusters in chromosomes 4q and 5q31 respectively to risk of scarring trachoma and trichiasis in a large case-control association study in a Gambian population. Methods Linkage disequilibrium (LD) mapping was used to investigate risk effects across the 4q and 5q31 cytokine clusters in relation to the risk of scarring sequelae of ocular Ct infection. Disease association and epistatic effects were assessed in a population based study of 651 case-control pairs by conditional logistic regression (CLR) analyses. Results LD mapping suggested that genetic effects on risk within these regions mapped to the pro-inflammatory innate immune genes interleukin 8 (IL8) and granulocyte-macrophage colony stimulatory factor (CSF2) loci. The IL8-251 rare allele (IL8-251 TT) was associated with protection from scarring trachoma (OR = 0.29 p = 0.027). The intronic CSF2_27348 A allele in chromosome 5q31 was associated with dose dependent protection from trichiasis, with each copy of the allele reducing risk by 37% (p = 0.005). There was evidence of epistasis, with effects at IL8 and CSF2 loci interacting with those previously reported at the MMP9 locus, a gene acting downstream to IL8 and CSF2 in the inflammatory cascade. Conclusion innate immune response SNP-haplotypes are linked to ocular Ct sequelae. This work illustrates the first example of epistatic effects of two genes on trachoma.


Background
Trachoma is the leading infectious cause of blindness. Eighty-five million children have active (inflammatory) trachoma, and about 7 million people, mainly adults, are blind from late scarring sequelae [1]. For demographic reasons, the number of people blind due to trachoma may be increasing [2]. In trachoma, ocular Chlamydia trachomatis (Ct) infection causes inflammatory changes in the conjunctiva, and repeated infections sometimes lead to fibrosis and scarring of the sub tarsal conjunctiva. This may cause the upper eyelid margin to turn inwards so that the lashes rub against the eyeball (trichiasis), which damages the cornea and leads ultimately to blindness [3].
In trachoma endemic areas most individuals suffer from ocular Ct infection in childhood. The majority resolve the infection without permanent sequelae. However, some individuals develop severe and persistent clinical disease in response to infection [4,5], and are more likely to develop conjunctival scarring and trichiasis in later life [6]. The reasons for this heterogeneity in susceptibility to chlamydial infection and disease progression following a rather uniform bacterial exposure remain incompletely understood; however the host genetic background may account for part of the risk. Some susceptibility loci have been mapped to cytokine genes [7][8][9][10][11][12][13], but most of the heritable risk effects in scarring trachoma and trichiasis are uncharacterised. Identification of genetic loci which contribute to the burden of trachoma sequelae in the population may help to dissect the immunopathology of ocular Ct infection in humans.
It has been postulated that the early secretion of proinflammatory cytokines and chemokines by epithelial cells following Ct infection may initiate and sustain a chronic inflammatory process associated with pathology [14][15][16]. The chromosomal regions 4q and 5q31, which are rich in chemokine and cytokine loci respectively, have been implicated in susceptibility to several common diseases [17][18][19]. They contain the genes interleukin 8 (IL8) and granulocyte-macrophage colony stimulatory factor (CSF2) which have strong candidacy to affect risk of the scarring complications of trachoma. However, studies of diversity in these chromosomal segments in the Gambian population revealed complex genetic features, including low haplotypic diversity of IL8 [20,21] and a long-range linkage disequilibrium (LD) pattern spanning CSF2 [22], both of which extend for several hundred kilobases. Therefore, observed disease association at the IL8 and CSF2 loci may potentially map to one or more of the neighbouring genes. We report here an analysis of the risk effects at IL8 and CSF2, in a Gambian case-control association study of scarring trachoma and trichiasis, which took the LD structure at these loci into account.

Methods
Study population 651 subjects with scarring trachoma, identified by clinical examination using WHO criteria and unrelated controls with normal eyelids pair-matched by age, sex, ethnicity and village of residence, were recruited for the study. These were genotyped in two phases; in the first phase, 373 case-control pairs were typed using the full marker sets. 269 of the cases had trichiasis in addition to scarring. Mean age was 39 years (range 5-90), and 273 pairs were female (73%). 32% were of Mandinka ethnicity, 30% Wolof and 25% Jola. In the second phase, 278 new casecontrol samples were further typed with the reduced marker set, 38 cases had trichiasis in addition to scarring. Mean age was 34 years, and 193 pairs were female (68%). 5% were Mandinka, 56% Wollof, 24% Fula and 9% Jola. The sample populations of the two case-control phase studies (phase 1+ phase 2) differed slightly in age group and ethnic composition. In the pooled dataset, 651 Gambian subjects with scarring trachoma (307 of them also had trichiasis), and pair-matched controls with normal eyelids, 268 individuals were Mandinka (25.8%), 542 were Wolof (33.3%), 236 were Jola (22.6%), 164 were Fula(10.8%) and 90 (7.5%) were of other ethnic origin. 69% of the recruited subjects were female and the mean age of the population under study was 37 years.
Subjects gave written consent. The study and its procedures were approved by the Gambia Government/MRC Ethics Committee (SCC 729/857), the Ethics committees of the London School of Hygiene and Tropical Medicine and of Oxford University, and are in accordance with the Declaration of Helsinki. Subjects diagnosed with trichiasis were offered free corrective surgery.

DNA extraction and SNP genotyping
Genomic DNA (gDNA) was isolated from either venous blood in EDTA or buccal brush samples. Genotypes were determined by the Sequenom system using matrix assisted laser desorption/ionization time-of-flight (MALDI-TOF) mass spectrometry as previously described [23]. Primer sequences are available on request.

SNP selection
Initially 9 haplotype tagging SNPs (htSNPs) spanning a region 330 kb long across the IL8 locus in the Gambian population [21] (Table 1) and 16 informative htSNP markers capturing the haplotype block structure within a 656 kb segment of the 5q31 region in Gambians [22] ( Table 2) were genotyped on a subset of 373 case-control pairs. From these marker sets, two reduced sets of 4 htSNPs covering the IL8 and CSF2 loci respectively were typed on the remaining 278 case-control pairs.

Allele and genotype frequencies
Estimates of allele and genotype frequencies were determined by dividing the total number of each observed allele (or genotype) by the total number of chromosomes in the population sample. Genotype distributions at all loci were in Hardy-Weinberg equilibrium (HWE) in both cases and controls.

LD estimates
The program HaploXT http://www.sph.umich.edu/csg/ abecasis/gold/docs/haploxt.html was used to calculate pair-wise LD estimates for the markers, expressed as the abs D' and R2 parameters (the latter being independent of the allele frequencies in the population). The input file for HaploXT was the haplotypes for each control subject generated by PHASE v2 [24]. The HaploXT output file was used as an input for MARKER v1.0, a graphical output program written by D. Kwiatkowski at the Wellcome Trust Centre for Human Genetics (WTCHG, Oxford) http:// gmap.net/. These estimates were used to construct haplotype blocks of high LD (and reduced within block diversity) within the regions studied.

Haplotype construction
The program PHASE v2 was used to infer haplotypes from population genotype data [25] and to estimate the frequency of each inferred haplotype for cases and controls.
To check the accuracy of haplotype reconstruction PHASE inferred haplotypes were compared with those from the SNPHAP program for haplotype construction from population data http://www-gene.cimr.cam.ac.uk/clayton/soft ware/snphap.txt.

Association analysis
An adjusted univariate analysis was first carried out by using Mantel-Haenszel (M-H) chi-square statistics (or Fisher's exact test if appropriate) to test for differences in minor allele, genotype and haplotype frequencies between cases of trachomatous scarring (TS) and of trichiasis (TT) with their pair-matched unrelated controls. Subsequently conditional logistic regression (CLR) analysis for disease association, taking into account matching of case-control pairs, was performed. In the presentation of results, reference genotypes were generally selected to be those that were present at the highest frequency in our study population. Genotypes conferring susceptibility or protection are accordingly represented by odds ratios (OR) of greater than or less than one respectively. A test for trend in the odds ratios was carried out to check for a dose response effect relationship between genotype/haplotype and disease allowing for potential confounders. All analysis was performed using STATA (v8.0) software.

Risk localisation analysis
To localise observed risk effects within the inferred population risk haplotype, PHASE assignments from a total of 764 population chromosomes from unrelated Gambian individuals were made for differing DNA segments by progressively reducing the number of nearby SNPs used to infer the population haplotypes. Association analyses were repeated each time until the risk effects (as measured by ORs) associated with each DNA segment delimited by those SNPs started to decrease.

Epistasis
The simplest system of two interacting loci using conditional logistic regression was used to assess epistasis between MMP9Q279R genotypes, which we previously reported as conferring risk to scarring trachoma and trichiasis, and genotypes at the unlinked IL8 and CSF2 loci.
Because the combination analysis leads to small cell sizes, we also tested for allelic epistasis, whereby joint risk effects of two alleles were tested for deviation from additivity. To assess evidence of epistasis between alleles at two loci (SNP1 and SNP2) a likelihood ratio test (LRT) comparing nested models with and without the interaction was performed. The LRT is based on the likelihood ratio statistic (LRS) with a χ 2 distribution with 1 degree of freedom: LRS = 2(L1-L0) where L1 is the maximum log likelihood including an interaction term between SNP1 and SNP2 and L0 is the log likelihood under the null hypothesis of no effect modification of SNP2 on the risk effects of SNP1, L0 = α+β1 (SNP1)+β2 (SNP2). Genotypic interactions were analysed similarly but allowing for the different types of main effect for each genotype (table 1) Values of p = 0.05 were considered statistically significant. P-values have not been adjusted for multiple testing. Analysis was performed using the STATA v8 package.

(i) Estimation of risk effects across a 330-kb region of chromosome 4q
DNA from 373 Gambian case-control pairs was genotyped at 9 SNP sites across 330 kb in chromosome 4q (table 1, Figure 1). Full-length population haplotypes inferred using PHASE v2 and population haplotype frequencies determined among cases and controls are shown in table 2. Pair-wise comparisons of R 2 between the 9 typed markers are similar to those identified in the Yoruba ( Figure 1) and Caucasian populations [20,21] with two broad clusters of association of high-to-moderate pairwise LD (R 2 > 0.1); cluster 1 including AFP and AFM genes and cluster 2 including IL8 and its downstream flanking region (data not shown).
Three out of 11 inferred haplotypes were present at population frequency ≥10% and accounted for 60% of all the haplotypes present in the population (table 2), suggesting a low haplotype diversity in this genomic region in Gambians. The haplotype AFP+8865G/AFM+1666A,/ AFM+4530A/AFM+15790T/IL8-251A/IL8+396T/ IL8+37674T/IL8+39739A/IL8+40050C, or GAATATTAC, was associated with increased risk of scarring trachoma, at borderline of statistical significance, and the haplotype AGGTATTAC was associated with significant protection from disease (table 2). In CLR analysis the risk of scarring increased with copy number of the GAATATTAC haplotype and it decreased with increasing copy number of the AGGTATTAC haplotype (table 3). The risk effects detected by the high risk haplotype were not affected by those observed for the low risk haplotype (and vice versa) (data not shown), suggesting that at least two independent risk effects may operate in the region.
Pair-wise measures of linkage disequilibrium and regression modelling (figure 1, table 3) mapped the risk associations to the four contiguous SNPs bounded by AFM+1666 and IL8-251 (AFM+1666A/AFM+4530A/ AFM+15790T/IL8-251A OR(95%CI) = 2.00 (1.30, 2,94) and 1.610 (0.99, 2.61) for scarring and trichiasis respectively). CLR analyses suggested that the AFM+1666 and AFM+4530 SNP sites explain mainly the high risk effects detected by the risk AATA (H2) haplotype whereas IL8-Top: Hapmap data of LD spanning a 330 kb segment of chromosome 4q in the Yoruba (Nigerian) population (LD map drawn using MARKER (see methods)) Figure 1 Top: Hapmap data of LD spanning a 330 kb segment of chromosome 4q in the Yoruba (Nigerian) population (LD map drawn using MARKER (see methods)); LD between each SNP pair is colour-coded (red dots represent an absolute R 2 > 0.3, green dots indicate R 2 of 0.1 to 0.3, grey dots indicate R 2 of 0.05 to 0.1 and absence of dots represents pairs with R 2 ≤ 0.05). Middle: Genes on the segment are indicated as black boxes, the position of the 9 SNP markers as vertical lines below the boxes, with arrows indicating direction of transcription. Bottom: Haplotype structure of the genomic segment in Gambians illustrating four haplotype blocks of high within-and low between-cluster diversity and the haplotype frequencies within each block beneath. Dots and shading indicate the configuration of the protective and risk haplotypes respectively. The dark shading spanning the IL8 locus shows the configuration in blocks 3 and 4 which is shared by the risk and protective haplotypes, differing only at SNP positions 1, 2 and 3 across the AFP and AFM loci but sharing allelic configuration at sites within, and telomeric to, IL8. 251 SNP tags some of the risk effects captured by the GGTA (L2) haplotype (data not shown).
When the risk haplotype tagging SNPs AFM+1666, AFM+4530, AFM+15790 and IL8-251 were retested for disease association in the complete case-control sample set of 651 case-control pairs, the AATA haplotype, present in 23% of the population, was found to be associated with an increased risk of trichiasis (OR = 1.39, 95%CI = 1.03, 1.88) with risk effects increasing 40% for each copy of the haplotype (OR = 1.40, 95%CI = 1.04, 1.87). On the other hand the rare IL8-251TT genotype (2.4% population frequency) was associated with reduced risk of scarring trachoma (OR for scarring 0.29, 95% CI = 0.09, 0.87, p = 0.027 and for trichiasis 0.50, 95%CI = 0.04, 5.51, p = 0.571 for IL8-251 TT versus IL8-251 AA+AT). The relative risk of developing scarring among IL8-251 TT homozygotes was estimated to be 2 times lower in the younger group (<35 years old) than in the older group (≥ 35 years), although this difference in risks was not statistically significant possibly due to small sample sizes (data not shown). The risk effects of the AATA haplotype were independent of the IL8-251 TT protective effects suggesting that the risk tagged by the haplotype may be independent of that detected by the IL8-251 TT genotype (data not shown).

(ii) Epistasis between IL8-25and MMP9 Q279R SNPs
There was some evidence for interaction between the MMP9 Q279R and IL8-251 alleles affecting risk of scarring trachoma (LRT χ 2 = 4.14, p = 0.042). The protective effect of the IL8 T allele was significantly increased, by more than 70%, in subjects carrying the MMP9 G protective allele (OR, 95%CI = 0.352 (0.14, 0.87) but suppressed in the presence of the MMP9 Q279R A risk allele (OR, (95%CI) = 1.90, (0.67, 1.21). No significant interaction effects were noted at the genotype level, although the number of IL8-251TT homozygotes is small. There was no evidence for gene-gene interaction involving MMP9 Q279R and either AFM+1666 or AFM+4500 (data not shown).

2(i) Estimation of risk effects across a 656 kb region of chromosome 5q
In brief, due to the irregular LD between the IL3/CSF2 and IL13/IL4 regions (Figure 2), genotyping data were used to construct 8 marker haplotypes across each of these two regions in 746 unrelated chromosome pairs (Table 4). LD statistics obtained from the haplotype data were in agreement with published data [22]. A significant association was observed between two 8 SNP-haplotypes delimiting a genomic block of high LD across IL3 and CSF2, and risk of scarring (Figure 2; tables 5 and 6). No significant risk effects were detected across the IL13 and IL4 loci (data not shown).
Following a similar approach to that of the IL8 region, PHASE assignments on progressively reducing SNP sets fine-mapped SNP-haplotype risk effects to 3 contiguous haplotype tagging SNPs bounded by CSF2_27348 and CSF2_27438 spanning the IL3 and CSF2 loci, (figure 2, tables 5 and 6). The high and low risk extended haplotypes shared the SNP alleles in IL3 (tables 5 and 6).
The CSF2_27348A/CSF2_2069614T/CSF2_27438A haplotype (ATA), of 17% population frequency, was found to be associated with lower risk of scarring and trichiasis (OR, (95%CI) of 0.82, (0.65, 1.04) and 0.61, (0.42, 1.88) for scarring trachoma and trichiasis respectively), decreasing with copy number (OR, 95%CI for test for trend of 0.82(0.66, 1.01) and OR of 0.65 (0.47, 0.89) for scarring trachoma and trichiasis respectively). Among common population haplotypes, ATA is unique in that it exclusively carries the CSF2_27348 A low risk allele. The risk effects detected by haplotype and single-marker analyses are very close (tables 5, 6, 7 &8) adding evidence to the idea that ATA haplotype and the CSF2_27348 A allele may capture the same risk effects within the region. This was further confirmed by logistic regression analyses showing that the haplotype risk effects are not significant when adjusted for the risk effects of CSF2_27348 (and vice versa) (data not shown). The risk associated with the CSF2_27438 SNP is best described by a recessive model (GG vs AG + AA genotypes). A protective effect at CSF2_27348 was also most marked in a recessive model for the minor genotype AA vs TT+AT. There was some evidence for a dose-response effect towards decreasing risk of disease with increasing number of copies of the CSF2_27348 A minor allele (tables 7 and 8).

Discussion
By means of combining LD mapping, case-control and epistatic analyses we have identified a genetic association between severe outcome of ocular Ct infection and the pro-inflammatory genes IL8 and CSF2 within the chemokine and cytokine gene clusters in chromosomes 4q and 5q31 respectively. These two pro-inflammatory mediators characterise an early host response to Ct infection thought to act to augment the innate immunity essential for host defence or, if unchecked, cause cytotoxic damage to the epithelium [14,15,26].
Homozygotes for the IL8-251T promoter region allele showed a greatly reduced risk of scarring trachoma. However the prevalence of the IL8-251TT genotype was low (2.4%) in our population indicating that a larger casecontrol series will be needed to get accurate estimates of risk effects. The degree of protection was greater in younger compared to older age groups suggesting that effects linked to IL8-251 may be particularly relevant to the development of scarring trachoma and trichiasis at an SNPs given in bold were selected for typing in the complete sample set of 651 case-control pairs a SNPs given in bold were selected for typing in the complete sample set of 651 case-control pairs.   Figure 2 Top: LD map across a 656 kb segment of the 5q31 region in the Gambian population [22]. Middle: Distribution of the 16 most informative markers ( Table 4) that capture the haplotype block structure extending from IL3 to IL4. Bottom; Haplotype structure defined by 8 SNP markers across IL3 and CSF2 illustrating two haplotype blocks and haplotype frequencies within each block beneath. Dots and shading indicate the configuration of the protective and risk haplotypes respectively; the dark shading spanning the IL3 locus shows the haplotype configuration shared by the risk and protective haplotypes.
early age. Cases and controls were matched by age and village of residence, so that this age difference is not accounted for by different treatment history (treatment is programmed per village) or environmental influences. The IL-251A common allele has been correlated with high IL8 expression in vitro [27] and also with an increased susceptibility to respiratory syncytial virus (RSV) bronchiolitis [21,27], a chronic inflammatory disease characterised by high levels of IL8 and neutrophil infiltration of the airway epithelium. It is therefore plausible that, in IL8-251T homozygotes, genetically controlled production of low IL-8 levels by Ct infected cells in the conjunctival epithelium leads to decreased activation and migration of neutrophils to the site of infection. The subsequent reduction in inflammation and cytotoxic epithelial damage, may ultimately protect the conjunctiva from the development of severe trachoma sequelae. The genetic risk effects for scarring and trichiasis observed in this genomic region are similar, which suggests that the above molecular mecha-nism may be common to the development of both phenotypes, rather than implicated in the progression of scarring to trichiasis.
Fine mapping of susceptibility to RSV bronchiolitis at this locus [21] suggested that additional risk variants in the region mapped to a region excluding the AFM gene (which codes for an albumin-like protein expressed in the liver) but including the novel gene RASSF6, a member of the Ras superfamily of small GTPases that can mediate bacterial entry by controlling changes in the actin cytoskeleton [28][29][30]. Here we show that the risk effect on trachoma maps to the region bounded by the AFM+1666 and IL-251 markers which contains RASSF6 ( Figure 1). Interestingly, by means of gene expression microarray analyses, we have found that the RASSF6 expression pattern maps to a cluster of genes involved in reorganization of the actin cytoskeleton whose expression is up-regulated  The functionality of the 2 intronic risk variants detected here across CSF2 is not obvious from genome annotation and it is possible that these SNPs act only as markers for a yet unidentified functional element. The CSF2_27348 marker has been found to be linked to the rs721121 SNP, located in an enhancer site for the inducible transcription factor AP1, an intergenic enhancer required for the correctly regulated activation of both CSF2 and IL3 gene expression in T cells [31]. Among other genes neighbouring CSF2 that need to be considered as candidates for the origin of the risk effects is SLC22A4 (Figure 2). This gene codes for an organic cation transported in lymphoid organs [32] which has been associated with risk of rheumatoid arthritis [32]. We have recently found that the transcriptional levels of SLC22A4 are increased in the conjunctiva of subjects with inflammatory trachoma and ocular Ct infection, mapping to a gene cluster characteristic of an innate immune response pattern of expression (unpublished data). Similarly to RASSF6, the role of SLC22A4 in the inflamed conjunctiva remains to be elucidated.
The statistical interaction (epistasis) observed between the risk variants within IL8 and CSF2 and the Q279R exonic SNP in MMP9, a collagenase thought to be involved in the inflammatory and scarring processes of Ct infection [11,33,34], may further suggest how genetic variation detected in IL8 and CSF2 could affect risk of trachoma: IL8, CSF2 and MMP9 are co-expressed in the Ct infected conjunctiva [34], unpublished data). These gene products could interact at the site of infection to augment and sustain inflammatory processes. MMP9 has been found to greatly enhance the activity of IL8 by amino-terminal processing [35], whereas activation of neutrophils by IL-8 may trigger the release of pro-MMP9 [36,37], creating a potential for a positive feedback loop. Similarly, CSF2 release by Ct infected epithelial cells [15,16] may mediate the influx and activation of inflammatory cells at the site of infection. Secreted CSF2 may trigger MMP9 production by monocytes [38], which could act to enhance and sustain the pro-inflammatory cascade initiated by CSF2. Thus the epistatic genetic risk effects could reflect biological interactions between their products, allowing a mechanism by which genetic variation affecting IL8 and CSF2 expression/activity change the risk effects associated with MMP9 variants and vice versa. Epistasis has not been reported before in studies of susceptibility to trachoma but may be an important feature of highly regulated gene networks.

Conclusion
This study provides the first evidence that genetic susceptibility determinants for severe sequelae of ocular Ct infec-tion may exist among innate response genes. This work adds to recent observations [39,40] suggesting that innate immune responses in early inflammatory events during acute infection with Ct may have both favourable and detrimental effects in the development of irreversible sequelae of Ct infection. Although we narrow down the limits of disease association to a few candidates within these genomic regions, which are complex and rich in immune genes, the results should be interpreted with caution. Functionally important SNPs directly influencing risk of disease remain to be uncovered; the complex and locusspecific architecture of the 4q and 5q31 genomic regions requires 'deeper' genotyping with a higher density of markers and increased sample sizes.