This study was reviewed and approved by the Tanzanian National Institute for Medical Research Ethics Committee, the Kilimanjaro Christian Medical Centre Ethics Committee, and the London School of Hygiene and Tropical Medicine Ethics Committee. It adhered to the tenets of the Declaration of Helsinki.

The study was conducted in three adjacent trachoma endemic villages in Kilimanjaro and Arusha regions, Northern Tanzania. The villages are relatively remote, geographically neighbors and have similar patterns of life and traditions. This area is predominantly inhabited by people of the Maasai tribe. Pastoral activities are the main occupation. The area is dry for much of the year, except for the rainy season (February to May). Water supply is therefore limited, and largely depends on a long-distance water pipe scheme from Mount Kilimanjaro. Family units are organized in Bomas, with living huts arranged in a circle around a central animal enclosure, which is often characterized by a high density of flies.

In January 2012, we recruited a cohort of children aged 6–10 years from these villages. The cohort has subsequently been followed-up every 3 months for 4 years to investigate the pathogenesis of conjunctival scarring. All children, aged 6–10 years, who were normally resident in one of the three villages, were eligible for inclusion. We chose this restricted age group as we considered that they were more likely to show evidence of incident or progressive conjunctival scarring during the 4 years of the study. At the outset community meetings were held to introduce the study. Each household was then visited to meet the parents or legal guardians of children eligible for enrolment. A field worker explained the nature of the study in detail in either Swahili or Maasai language. There was an opportunity to discuss and ask questions. Finally, if the parent or guardian agreed to allow the child to be enrolled into the study this was documented on a consent form in Kiswahili, and witnessed by a third person.

The left eye of each child was examined by an ophthalmic nurse experienced in grading trachoma. The examinations were all performed under standardized conditions, using x2.5 loupes and a bright touch. The conjunctiva was anesthetized with a drop of preservative-free proxymetacaine hydrochloride 0.5%w/v (Minims®, Chauvin Pharmaceuticals Ltd, Surrey, UK). The eyelid was everted, examined and photographed (Nikon D90 with 105 mm Macro lens). Two conjunctival swab samples (Dacron polyester, Puritan Medical Products Company, Maine, USA) were collected for C. trachomatis detection and gene expression analysis. The swabs were passed across the upper tarsal conjunctiva four times, with a quarter turn between each pass. The first swab was placed directly into a tube containing RNAlater (Thermo Fisher Scientific, Massachusetts, USA) and the second into a dry tube. The samples were placed into a cool box. Later the same day the dry swab samples were stored directly at −80°C and the RNAlater samples kept at 4–8°C overnight and then stored at −80°C. Air control swabs were collected after every 50 samples by passing a swab 10 cm from a participant's everted eye, these were labeled and processed identically to participant samples.

Clinical signs were graded using the 1981 WHO Detailed Trachoma Grading System (FPC) (Dawson et al., 1981). This sub-divides the features into several four-point severity scales: follicles (F), papillary inflammation (P), and conjunctival scarring (C). This system corresponds to the WHO Simplified Trachoma Grading System in the following way: Trachomatous inflammation-Follicular (TF) is equivalent to F2/F3 and Trachomatous inflammation-Intense (TI) is equivalent to P3 (Thylefors et al., 1987). For the purpose of this study, we consider that both P2 and P3 represent clinically significant papillary inflammation, and refer to this as “TP” (Burton et al., 2015). We followed the widely used definition of “Clinically Active Trachoma”: TF and/or TI. The digital photographs were graded by an ophthalmologist experienced in trachoma assessment, using the FPC system, supplemented by a previously described system for fine grading of conjunctival scarring, that quantifies the extent of the conjunctiva involved (Hu et al., 2011). The ophthalmologist was masked to the infection status.

Extensive public health education about trachoma was provided to the community through village level meetings and during the house-to-house visits, including the importance of face washing for children and environmental improvements. Trichiasis surgery was provided free of charge within the community. Subsequently, all residents of the three villages have been offered three rounds of annual mass antibiotic treatment with oral azithromycin (and topical tetracycline ointment for infants under 6 months and pregnant women), during the course of the longitudinal study.

To detect C. trachomatis, we extracted DNA from swab samples stored in dry tubes using the PowerSoil DNA Isolation Kit (Mo Bio Laboratories, California, USA), according to the manufacturer's instructions. The cells attached to the swabs were initially disrupted by bead beating to release their contents. C. trachomatis DNA was detected using a previously described droplet digital PCR assay (Roberts et al., 2013). All samples were tested for the C. trachomatis plasmid and the human gene RPP30. Samples in which C. trachomatis plasmid was detected were retested for C. trachomatis omcB (a single-copy gene on the chlamydial chromosome). Five microliters of template DNA was added per reaction. PCR reaction conditions were as follows: 95°C for 10 min, then 40 cycles of 95°C for 10 s and 60°C for 30 s and finally 98°C for 12 min. Droplets were then examined for fluorescence on a QX100 Droplet Reader (Bio-Rad, UK), providing a quantitative result.

Total RNA was extracted from the swab stored in RNALater using a DNA/RNA Purification Kit (Norgen Biotek Corp, Canada), following the manufacturer's instructions. RNA was reverse transcribed using the SuperScript® VILO® cDNA Synthesis Kit (Life Technologies). Quantitative real-time PCR was performed to measure the relative abundance of a panel of human transcripts, using the TaqMan® Microfluidic 384-well Low-Density Array (TLDA) (Life Technologies) on a ViiA7 real-time PCR instrument (Thermo Fisher Scientific, Massachusetts, USA). These cards have separate channels for eight samples, each containing 48 wells, pre-printed with the assay primers and probes. We used two different TLDA designs and measured the expression of 91 different genes of interest for each sample. The selected genes are listed in Table 5. We measured expression of three different reference genes: HPRT1, GAPDH, and RPLP0. HPRT1 was selected as the most suitable for normalization as it was expressed at a relatively similar level to the majority of other transcripts of interest, whereas GAPDH and RPLP0 were very highly expressed. The choice of transcripts tested was informed by previous gene expression studies including transcriptome analysis experiments on samples from The Gambia, Ethiopia and Tanzania and by a genome-wide association study from The Gambia (Holland et al., 2010; Natividad et al., 2010; Burton et al., 2011b,c, 2015; Roberts et al., 2015).

Data were managed in Access and transferred to STATA v14 for analysis. The ophthalmologist's grading of the digital photographs was used for analysis. The ddPCR results were initially exported to R for analysis; the critical cut-off value for designating a positive result for each target was 0.2 plasmid copies/μl of eluted DNA, 0.2 omcB copies/μl and >0.3 RPP30 copies/μl. Samples were designated C. trachomatis positive if plasmid and human RPP30 were detected.

The associations between each of the clinical signs (follicular inflammation, papillary inflammation and scarring) and C. trachomatis infection were assessed using univariable logistic regression with infection as the outcome variable. The associations between scarring and each of the other clinical signs, infection, sex and age were investigated with logistic regression, using presence of scarring as the outcome variable. Initially, univariable logistic regression was performed using each exposure in turn, and then a multivariable analysis was performed in order to provide an unbiased estimate of the association between exposure and outcome adjusting for all clinical features.

The gene expression data were normalized relative to the expression of HPRT1 in the same sample, to adjust for variable concentrations of cDNA. This was done by the ΔC_T method. Distributions of ΔC_T values were plotted to assess them for normality. Fold change differences in gene expression between different phenotypic groups (Infection, follicular inflammation, papillary inflammation and scarring), vs. individuals without those clinical signs or infection, were calculated using the ΔΔC_T method. Linear regression was used to compare the two groups in each comparison, with ΔC_T values as the dependent variable and each phenotype group in turn as the independent variable, adjusting for age and sex. To take account of multiple comparisons we used a false discovery rate (FDR) of 5% (Benjamini and Hochberg, 1995). The linear regression analyses were repeated, but rather than taking the four phenotypic groups individually, they were all included as exposures in a multivariable model with age and sex, providing estimates of the fold-changes associated with each phenotype.

A heatmap of ΔC_T values for each gene of each participant was produced, to indicate visually whether the expression of individual genes was consistently of higher or lower expression in groups of individuals with clinically active Trachoma (TF and/or TI) and/or infection. A principal component analysis (PCA) was performed on the gene expression data and clinical disease and infection status were overlaid on plots in order to visually indicate any associations between the first two components and infection, active trachoma and scarring. These associations were tested formally using a logistic regression analysis, with the principal components as the exposure variables. For the heatmap and PCA, only individuals with complete expression data could be included. Therefore, so as not to lose too many individuals, any target with >5% of observations missing was excluded from these analyses, as well as the lasso regression and co-expression analyses described below. Any individuals with missing observations on the remaining targets were excluded from these analyses also.

For each of infection, active trachoma and scarring, a multivariable logistic regression was performed using this subset of gene expression levels as exposure variables, with the aim of retaining those expressions most strongly associated with the outcome of interest, adjusting for all other expression levels. Due to the large number of exposure variables, a penalized logistic regression was performed using the lasso technique (Tibshirani, 1997). The list of genes most strongly associated with each outcome was entered into ConsensusPathDB (http://cpdb.molgen.mpg.de/) gene set over-representation analysis, with all genes tested entered as background. Enriched pathway-based sets with a minimum overlap between pathway and input gene list of two and a p < 0.05 were identified.

A network graph based on the specimen-to-specimen Pearson correlation was generated using miru (https://kajeka.com/miru). The overall expression correlation matrix and graph were constructed from the raw cycle threshold expression values. The filtered dataset including only individuals and genes with complete expression data was used. Pearson correlation coefficients (r) ≥0.85 were retained and used as cut-offs in network construction. Nodes in the graph are individual mRNA transcripts linked by an edge if r was ≥0.85. The graph was then clustered using a Markov Clustering algorithm with an inflation value of 2.2. The partitioned clusters of expression contain sets of transcripts that exhibit a very high degree of co-expression across the sample set. The co-expression clusters or modules were then investigated for enrichment at the pathway level using ConsensusPathDB as described above, using all genes tested as background. A single value for each cluster of each individual was defined using the first principal component and this value was used to investigate differential expression associated with disease and infection phenotypes.

A total of 506 children between the ages of 6 and 10 years were recruited and examined. Their demographic and clinical characteristics are described in Table 1. There were almost equal numbers of male and female children, with a mean age of 7.5 years. Most (97.4%) of the children were of the Maasai ethnic group.

Trachomatous inflammation-Follicular (TF, F2/F3) was present in 170 (33.6%) children (Table 1). However, a further 135 (26.7%) had evidence of mild follicular conjunctivitis (F1). Significant conjunctival papillary inflammation (TP, P2/P3) was observed in 160 (31.7%) children, of which 64 (12.7%) had intense papillary inflammation (TI, P3).

C. trachomatis plasmid was detected by ddPCR in 78 (15.4%) individuals. There was a strong association between the presence of TF and the detection of C. trachomatis plasmid; 62/78 (79.5%) of C. trachomatis positive individuals had TF (Table 2). Of the 16 infected children who did not have TF, 11 (68.8%) had signs of mild follicular conjunctivitis (F1). However, only 62/170 (36.5%) individuals with TF had detectable infection. There was a similar strong association between the presence of TP and chlamydial infection (Table 2).

Infection load and the concentration of RPP30 endogenous control DNA were quantified. The concentration of plasmid DNA ranged from 0.22 to 3023 copies/μl, omcB ranged from 0.2 to 742 copies/μl and RPP30 ranged from 0.32 to 1,081 copies/μl of eluted DNA. All samples had detectable RPP30. There was a marked positive trend in both the proportion infected and the load of infection with increasing F-Score and increasing P-Score (Table 2). There was a consistent ratio of plasmid to omcB copies/μl: geometric mean ratio of 4.57 to 1 (95% CI 3.91–5.34) and a median of 4.79 to 1 (95%CI 4.01–5.21). The correlation between plasmid and omcB copies/μl is illustrated in Figure 1. The ratio did not vary with either the follicular or papillary clinical severity scores (data not shown).

Conjunctival scarring was relatively frequent, being found in 144 (28.5%) individuals; this was mostly mild scarring (Table 1). There was a significant association between the detection of C. trachomatis and scarring [p = 0.002, OR = 2.23 (95%CI 1.4–3.7)]. The load of infection also increased with scarring severity (Table 2). There were univariate associations between conjunctival scarring and TF, TP (P2/P3), C. trachomatis detection and female sex (Table 3). However, in a multivariable model only TP, female sex and age were independently associated with scarring.

There was strong evidence (p = 0.0047) that even at this young age females are proportionately more likely to have signs of established scarring than males (Table 4). Similarly, there was some evidence (p = 0.048) that females had slightly increased odds of TP. No evidence was found of an association between sex and the odds of TF or chlamydial infection (p = 0.12 and 0.25, respectively) (Table 4).

The expression of 91 target genes was measured relative to that of HPRT1. Out of the 506 participants, samples from 12 individuals failed the qPCR for all target genes, leaving 494 individuals for gene expression analysis. The number of participant samples with detectable expression for each target is shown in Supplementary Table 1 (N). Multivariable linear regression analysis was performed using data for all 91 genes from 494 individuals.

Four sets of comparisons were performed: (1) infected v non-infected, (2) TF v no TF, (3) TP v no TP, and (4) TS v no TS. The relative Fold Change (FC) between these paired groups, calculated by the ΔΔC_T method, are presented in Supplementary Table 1 along with the p-values for the linear regression (adjusted for age and sex) for each comparison.

As clinical signs and infection may be highly correlated in their association with the expression of specific targets we tested for independent associations with each target's expression using multivariable linear regression models, Table 5. More genes were significantly differentially expressed in individuals with C. trachomatis infection (55/91) relative to individuals with TP (33/91), TF (11/91), and TS (17/91) in these multivariable models (P < 0.05). C. trachomatis Infection was associated with increased expression of multiple cytokines and chemokines (CCL2, CXCL13, CCL18, CSF2, FOXP3, IFNG, IDO1, IL1B, IL6, IL8, IL10, IL12B, IL17A, IL19, IL21, IL22, IL23A), cell cycle components (CDC25C, TTK, TYMS), matrix modifiers (MMP7, MMP9, MMP12, TGFB1), NK cell markers (CD247, NCR1, NCAM1) and intracellular signaling molecules/regulators (CD274, IKZF1, RHOH, SAMSN1, SERPINB3, SERPINB4, STAT1, STAT4, SOCS1, TBX21), Table 5. Infection was associated with particularly marked increases in expression of IFNG (FC: 7.12, p-value: 1.32E-42) and IL22 (FC: 6.49, p-value: 1.6E-23). There was reduced expression of several factors including mucins (MUC1, MUC4, MUC5AC, MUC7) and SPARCL1.

Clinical signs of active trachoma were associated with increased expression of multiple cytokines, chemokines, antimicrobial peptides, matrix modifiers, following a similar pattern to infection (Table 5). TP was associated with more substantial increases in a wider range of factors than TF (Table 5): antimicrobial peptides (DEFB4A, S100A7), cytokine/chemokines (CXCL5, CXCL13, CCL18, CCL2, IL1B, IL8, IL10, IL17A, IL19, IL23A, PTGS2), matrix modifiers (MMP7, MMP9, MMP12) and intracellular signaling molecules/regulators (CD274, SAMSN1, SERPINB3, SERPINB4, SOCS3). Several cytokines and chemokines (CCL18, CXCL13, IL19, and S100A7) had a FC >2 in children with TP. Conjunctival scarring was associated with modestly increased expression (>1.2 FC) and borderline significance of several chemokines and cytokines (CCL2, CCL20, CXCL5, CSF3, IL1B, IL8, IL10, IL21, IL33, FGF2) in the multivariable linear regression models (Table 5). SPARCL1 was significantly downregulated in TF and TP but not TS.

Only 108 individuals had detectable gene expression for all 91 targets, so in order to retain a reasonable number of individuals in the PCA and lasso analyses, which require complete records for all individuals, targets which were missing >5% of observations were excluded. This resulted in the exclusion of eight transcripts (USP6, IL13, CTNND2, FGF2, SERPINB4, IL22, PREX2v1, PREX2v2) for which there was no detectable expression in >5% individuals. In individuals with detectable expression of these eight targets, raw C_T values were generally very high (C_T > 33). These eight targets were excluded from the heatmap, PCA, lasso regression and co-expression network analyses, retaining a complete expression dataset for 83 genes. Among the 494 individuals, 36 were excluded as they had missing data among the 83 included targets, resulting in a final dataset of 458 individuals and 83 genes. Excluding further genes offered minimal gains in sample size. A heatmap of ΔC_T values separated by clinical and infection phenotype is presented in Figure 2. Infected individuals with or without active trachoma had a visibly distinct pattern of gene expression, whereas the ΔC_T values of normal healthy controls and individuals with active trachoma in the absence of infection appeared more similar overall.

PCA analysis was performed for 458 individuals that had complete gene expression data for 83 targets and active trachoma/infection and scarring status were overlaid on plots, Supplementary Figures 1A,B. Principal component 1 (PC1) was strongly associated with infection and active trachoma; for a one unit increase in PC1 there was an estimated 35% decrease in the odds of infection (95% CI 27–41%, p < 0.001) and an estimated 23% decrease in the odds of active trachoma (95% CI 9–28%, p < 0.001). PC1 loadings showed that as genes in the chemokines/cytokines and regulator/signaling pathways groups were up-regulated, PC1 became smaller, suggesting that the odds of infection/active trachoma was greater among individuals where these genes are up-regulated. The expression of mucins, CDH1, SPARCL1, and S100A4 showed the reverse relationship, as when these genes were down-regulated the PC1 was smaller. Principal component 2 (PC2) was associated with increased odds of active trachoma alone, with a unit increase in PC2 being associated with an estimated 15% increase in the odds (95% CI 7–24%, p < 0.001), however there was no evidence of an association between PC2 and infection. There was evidence of a weaker association between PC1 and scarring, with an estimated 10% decrease in the odds of scarring (95% CI 6–14%, p < 0.001) per unit increase in PC1, Supplementary Figure 1B. There was no evidence of an association between scarring and PC2.

Three lasso logistic regressions were performed, first using infection (with or without active disease) as the outcome and the 83 targets, age and sex of 458 individuals as exposures. Then, excluding all infected individuals, active trachoma (uninfected only) was used as the outcome against uninfected individuals without active trachoma and the same 83 targets, age and sex as the exposures of interest. Finally, the regression was performed using scarring as the outcome and the same set of exposures in all 458 individuals. The transcripts most strongly associated with each of the three outcomes, and biological pathways with members that were over-represented in each relative to the background (83 targets) are listed in Table 6. Genes most strongly associated with infection were associated with Th1 cell development and NK cell pathways, whilst genes most associated with active trachoma in the absence of C. trachomatis infection were enriched for members of leukocyte transendothelial migration, TNF receptor, MMP activation, collagen formation and collagen fibrils assembly pathways. Genes most strongly associated with the TS phenotype were enriched for members of immunoregulatory interactions between lymphoid and non-lymphoid cells.

Networks of co-expression, independent of differential expression, were explored in the complete expression data of 83 transcripts and 458 individuals using Miru. The undirected graph contained 51 nodes connected by 121 edges. Markov clustering partitioned the network into 5 separate clusters of co-expressed genes that accounted for 38/51 (74.5%) of the transcripts in the original network. Cluster 1 contained 11 co-expressed transcripts, cluster 2 contained 10 transcripts, there were 6 transcripts in each of cluster 3 and 4, and 5 transcripts in cluster 5 (Supplementary Table 2A), the remaining nodes were not connected in the graph. The most connected gene transcript or hub in each module was MyD88/REL (cluster 1), IKZF1/VIM (cluster 2), STAT1 (cluster 3), SRGN (cluster 4), and IL10/IL17A/IL21 (cluster 5).

Pathway enrichment analysis was performed for each cluster relative to the background of 83 transcripts (Supplementary Table 2B). ΔC_T values for the transcripts in each cluster were collapsed into a single value for each of the 458 individuals, represented by the first principal component. Multivariable linear regression was then performed for each cluster, using the first principal component as the outcome variable and TF, TP, TS, infection, age, and sex as independent variables (Supplementary Table 2C). Cluster 1 was enriched for cell cycle and apoptosis pathways, however the combined expression value was not differentially expressed in any condition. Cluster 2 was significantly associated with infection and TF and was marginally enriched for the retinoblastoma pathway. Cluster 3 was highly enriched for Th1/2 cell differentiation, IFNγ, IL-12, CD8, and NK cell signaling pathways and was strongly associated with infection and age and marginally associated with sex, but interestingly not with TP or TF. Conversely, cluster 4 was enriched for TNF, inflammasome, TLR and NFkB signaling pathways and was associated with TP, TS, and age but not infection. Cluster 5 was enriched for the allograft rejection pathway and was associated with infection, TP, TS and age.

In this study, conducted in communities prior to azithromycin MDA, there was a relatively high TF (34%) and moderate C. trachomatis (15.4%) prevalence in 6–10 year olds. There was a strong relationship between infection and clinical signs. The large majority (79.5%) of individuals with C. trachomatis detected by ddPCR had TF (F2/F3). Moreover, most individuals (68.8%) with detectable chlamydial infection but without TF had a mild follicular conjunctivitis (F1). There were few cases of infection in the absence of clinical signs. This is consistent with earlier studies which found a stronger correlation as the underlying prevalence of infection increased (Ramadhani et al., 2016b). However, in common with many other studies reporting the relationship between disease and infection, a minority (36%) of individuals with TF had detectable infection, probably due to the shorter duration of infection episodes relative to disease (Grassly et al., 2008; Burton et al., 2011a,c; Lee et al., 2014). With increasing severity of both TF and TP scores, there was an increasing proportion with infection and increasing loads of infection. Consistent with earlier studies, increasing load appeared to be more closely related to increasing TP rather than TF (Burton et al., 2003; Solomon et al., 2003; Michel et al., 2011; Derrick et al., 2016a).

In the present cross-sectional study conjunctival scarring was associated with TP, female sex and increasing age but not C. trachomatis infection, which is consistent with several earlier cohort studies (Dawson et al., 1990; West et al., 2001; Wolle et al., 2009b; Burton et al., 2015; Hu et al., 2016). However, longitudinal data supporting an association between TF or C. trachomatis infection and progressive scarring are limited (Ramadhani et al., 2016a). We found females had more TP and TS than males, whereas C. trachomatis infection and TF prevalence were not significantly different between the sexes (Table 4). This might suggest that there are additional determinants of TP and TS development beyond the initial C. trachomatis infection. It is generally believed that the difference in the proportion of males and females developing scarring sequelae of trachoma is attributable to a greater life-time exposure to repeated C. trachomatis infection among females (Taylor et al., 2014). While that may well be the case, the data from this study also suggest that even at a relatively young age, females appear more susceptible to developing TP and TS than males, despite comparable levels of C. trachomatis infection. This finding is consistent with the observation that in general, females generate stronger immune responses than males, making them more susceptible to diseases resulting from immune-mediated pathology (Klein and Flanagan, 2016).

Among infected individuals, the C. trachomatis plasmid:omcB ratio (4.57:1) was very similar to the ratio reported from a range of ocular and genital C. trachomatis isolates (mean 4.0:1) and from a population-based trachoma study in Guinea Bissau (5.3:1) (Pickett et al., 2005; Last et al., 2014). In common with the study from Guinea Bissau, we did not find evidence of variation in the plasmid copy number with increasing disease severity scores.

IFNG was strongly upregulated in association with C. trachomatis infection, consistent with previous findings in trachoma endemic populations (Burton et al., 2004; Holland et al., 2006; Natividad et al., 2010). Several NK cell markers were also upregulated in association with C. trachomatis infection: NCR1, CD56 (NCAM1), and CD247. IFNγ, NCAM1, and NCR1 were among the genes most strongly associated with increased odds of C. trachomatis infection, and Th1 and NK cell pathways were enriched in this gene set (Table 6). A strong IFNG response and transcriptional signatures associated with NK cell activation and cytotoxicity were closely associated with C. trachomatis infection and active trachoma in a West African population (Natividad et al., 2010), and infiltrates suggestive of NK cells were identified in scarred conjunctival biopsy tissue from individuals with trachomatous trichiasis from Tanzania (Hu et al., 2016). In peripheral blood from adults with trichiasis and controls that were stimulated ex-vivo with C. trachomatis antigens, NK cells and to a lesser extent CD8+ T cells were the major sources of IFNγ with <50% of the IFNγ produced by CD4+ T cells (Gall et al., 2011). In a non-human primate live-attenuated ocular chlamydial vaccine study, depletion of CD8+ cells abrogated protective immunity, however the contribution of CD4+ and NK cells was not investigated (Olivares-Zavaleta et al., 2014). IFNγ is a crucial component of the Th1 response to restrict C. trachomatis survival via tryptophan depletion, macrophage activation and promoting leukocyte adhesion, migration and activation (Beatty et al., 1994; Mosser and Edwards, 2008; Redgrove and McLaughlin, 2014). Interestingly, whilst IFNG was strongly upregulated in individuals with C. trachomatis infection, there was much less upregulation in individuals with disease but not infection, suggesting that the IFNG response is quickly down-regulated once C. trachomatis has been cleared. Supporting this observation, cluster 3, which was enriched for Th1, NK, IL-12, and IFNγ pathways, was strongly associated with infection but not with TF or TP. Cluster 3 expression was also marginally associated with female sex, and IFNG, IDO1, IL12B, and IL19 were each upregulated in females relative to males (Table 5). Both the number of activated CD4+ cells and their relative expression of pro-inflammatory genes including IFNG have previously been shown to be higher in females than males (Hewagama et al., 2009; Sankaran-Walters et al., 2013). Our data suggest that increased expression of IFNG and related pathways in females may contribute to the observed differences in tissue responsiveness (TP and TS) between the sexes in the context of equal infection exposure. In support of this hypothesis, a longitudinal study conducted in The Gambia that visited participants every 2 weeks found that the likelihood of developing infection or active disease and the duration of infection were independent of sex, however males had a reduced duration of active disease relative to females (Grassly et al., 2008). Whilst Th1 responses and IFNG are essential for controlling chlamydial infection, excessive responses are thought to cause collateral tissue damage and contribute to chlamydial pathology (Van Voorhis et al., 1997; Rank et al., 2000). It is therefore possible that a female predisposition to overproduce Th1 responses and IFNG may tip the balance from protection toward pathological tissue damage. However, further research is required to clarify this hypothesis and a direct causal relationship cannot be determined from this cross-sectional analysis.

The cytokine genes IL1β, IL6, IL17, IL21, IL22, IL23, and TGFβ1 were strongly expressed in individuals with C. trachomatis infection. Th17 cells are generated in the presence of TGFβ1, IL-1β, IL-6, IL-21, and IL-23 and classically produce IL-17, IL-22, and IL-21 (Manel et al., 2008; Zielinski et al., 2012). IL21 expression was also upregulated in individuals with TF and TS. Individuals with TP significantly upregulated IL17 and IL23 expression and IL21 at marginal significance. The combined expression of cluster 5, containing IL10, IL17A, IL21, CXCL13, and MMP12, was significantly associated with C. trachomatis infection, TF, TP, and TS. Th17 cells have important roles in mucosal immunity against bacteria and fungi, however they are also associated with the inflammatory pathology of diseases such as Crohn's disease, psoriasis and rheumatoid arthritis (Tesmer et al., 2008). Epithelial cells stimulated by IL-17 and IL-22 upregulate chemokines, cytokines, antimicrobial peptides and growth factors such as IL-6, IL-8, S100A7, and GM-CSF (CSF2) (Wolk et al., 2006; Eyerich et al., 2010); these transcripts were upregulated in individuals with C. trachomatis infection (S100A7 at marginal significance). IL8 and S100A7 were also upregulated in individuals with TP. IL-8 and S100A7, amongst other factors induced by Th17 cells, are chemotactic for neutrophils, an influx of which is associated with C. trachomatis infection and contributes to pathology in animal models of urogenital and ocular infection (Frazer et al., 2011; Lacy et al., 2011). IL-22, which was highly upregulated with infection, is also important in the maintenance of epithelial health and barrier function (Radaeva et al., 2004; Rutz et al., 2013). The relative contributions of IL-17 and associated responses to protection from C. trachomatis infection and maintenance of barrier function vs. pathological tissue damage remain unclear. In murine models of genital C. trachomatis infection higher levels of IL-17/Th17 cells have been associated with greater pathology (Lu et al., 2012; Vicetti Miguel et al., 2016), whereas in IL-23 knock-out mice which had greatly reduced levels of IL-17 and IL-22 there was no difference in infection burden or oviduct pathology compared to wild type mice (Frazer et al., 2013). Longitudinal investigation of these factors is expected to further our understanding of the role of Th17 cell associated cytokines in human chlamydial pathogenesis.

STAT1 and STAT4, which promote Th1 cell development (Nishikomori et al., 2002; Mikhak et al., 2006; Zhu et al., 2010), were upregulated in response to C. trachomatis infection. STAT1 and STAT4 have previously been shown to be associated with active trachoma and ocular C. trachomatis infection, and STAT1 was upregulated in cervical epithelial cells infected with C. trachomatis in vitro which was found to inhibit chlamydial growth (Lad et al., 2005; Natividad et al., 2010). STAT3, which is required for Th17 development (Zhu et al., 2010; Yang et al., 2011), was not differentially regulated, contrasting the upregulation of Th17 cytokines. SOCS1, IL10, and FOXP3 were upregulated in individuals with C. trachomatis infection and SOCS3 and IL10 were upregulated in individuals with TS and TP, indicative of inflammatory regulation. Upregulation of IL10 expression has previously been reported in individuals with active trachoma and a genetic polymorphism in IL10 was associated with trachomatous scarring (Mozzato-Chamay et al., 2000; Burton et al., 2004; Natividad et al., 2005; Faal et al., 2006; Skwor et al., 2008). However, the role of Tregs and IL-10 in trachoma remains unclear; dampening inflammation might result in less tissue damage, but conversely it might also prolong survival of C. trachomatis (Zhang et al., 2009).

MMP9, MMP12, TGFβ1, and PDGF were significantly upregulated in individuals with C. trachomatis infection. MMP7, MMP9, and MMP12 were upregulated in individuals with TP, MMP9 was upregulated in TF and FGF2 and PDGF were upregulated in individuals with TS. MMP7 was downregulated in infected individuals. MMP7 has consistently been found to be upregulated in scarring trachoma and trichiasis alongside MMP9 and MMP12 but there is less evidence for a role of MMP7 in active trachoma (Holland et al., 2010; Burton et al., 2011b, 2015; Hu et al., 2012). Degradation of the extracellular matrix and basement membrane by MMP9 and MMP12, in part derived from macrophage and neutrophil degranulation (Shapiro et al., 1993; Kjeldsen et al., 1994), facilitates connective tissue remodeling and leukocyte migration into the epithelium (Ozdemir et al., 1999; Zeng et al., 1999; Misko et al., 2002). Whilst MMPs are an essential part of the wound healing process (Nwomeh et al., 1998; Wolf et al., 2017), MMP9 overexpression has been linked to aggressive contractile scarring in proliferative vitroretinopathy and failure of rat corneal cells to re-epithelialize following corneal injury (Fini et al., 1996; Kon et al., 1998).

SPARCL1 was down-regulated in individuals with C. trachomatis infection, TF and TP, consistent with previous studies (Hu et al., 2012; Burton et al., 2015). SPARCL1 codes for a secreted glycoprotein which regulates extracellular matrix (ECM) and cell-matrix adhesion (Girard and Springer, 1996) SPARCL1 inhibits cell proliferation and migration and reduced SPARCL1 expression has been linked to cancer metastases (Nelson et al., 1998; Sullivan and Sage, 2004). In a murine corneal injury model, knock-out of SPARCL1 led to excessive accumulation of collagen type IV from myofibroblasts and irregular fibrotic ECM formation (Chaurasia et al., 2013). Addition of exogenous SPARCL1 rescued the wild type phenotype with a cessation of collagen IV production, replacement with collagen I from keratinocytes and restoration of normal ECM architecture. Consistent with this, immunohistochemistry and histology studies using tissue from individuals with scarring trachoma have demonstrated a reduction in collagens I and III, an increase in collagens IV and V and progressive disorganization of collagen fibers in scarring trachoma relative to healthy conjunctival tissue (Abu El-Asrar et al., 1998; Hu et al., 2013b; Derrick et al., 2016b). Genes that were most strongly associated with active trachoma in the absence of C. trachomatis infection were enriched for pathways of leukocyte transendothelial migration, TNF receptors, MMP activation, collagen formation and assembly (Table 6). Active trachoma signs often persist after the clearance of C. trachomatis infection (Grassly et al., 2008; Burton et al., 2011a,c; Lee et al., 2014) and inflammation is linked to progressive scarring (Ramadhani et al., 2016a). We found that TP but not infection was associated with TS in the multivariable model (Table 3). These results suggest that the inflammatory phenotype following the clearance of C. trachomatis infection is characterized by collagen formation and assembly, MMP-mediated matrix remodeling and leukocyte influx, and that these pathways might be causative factors in the scarring process.

We found significantly reduced expression of MUC1, MUC4, MUC5AC, and MUC7 in individuals with C. trachomatis infection, and MUC5AC expression was reduced in those with TP. MUC1 and MUC4 are primarily expressed by apical cells of the stratified epithelia, MUC7 by lacrimal gland epithelia and MUC5AC by goblet cells (Gipson and Argueso, 2003). MUC7 expression was also significantly reduced in females relative to males. MUC5AC and MUC7 down-regulation has previously been reported in individuals with trachomatous trichiasis, however expression of MUC1 and MUC4 was increased in this later stage of disease (Burton et al., 2011b). IL-22 stimulates the upregulation of MUC4 and MUC1, therefore, one might expect a corresponding increase in mucin expression in infected individuals (Turner et al., 2013). This inconsistency and the broad reduction in mucin expression during infection may reflect a loss of epithelial homeostasis and a potential reduction in barrier function. This in turn could result in increased access of bacteria to the epithelial cell surface, promoting inflammation. In patients with ulcerative colitis, reduction of goblet cells, mucin production and microbial diversity were identified as causative factors in disease etiology (Alipour et al., 2016). Microbial dysbiosis can itself be caused by inflammation, creating a positive feedback loop (Winter et al., 2013). Decreased diversity has previously been found in the conjunctival microbiome of individuals with scarring trachoma relative to healthy controls (Zhou et al., 2014). Furthermore, penetration of the mucin barrier and contact between epithelial cells and bacteria was demonstrated in ulcerative colitis patients with acute inflammation, whereas the mucin barrier of healthy controls was impenetrable (Johansson et al., 2014). It is tempting to speculate that reduced barrier function resulting from inflammation and loss of mucins causes microbial dysbiosis and contact between bacteria and epithelial cells, leading to persistent conjunctival inflammation, which in turn contributes to trachomatous scarring.

Consistent with our previous studies in various clinical stages of trachoma (Burton et al., 2011b,c, 2012, 2015; Hu et al., 2012), the antimicrobial peptide S100A7 was upregulated in individuals with TP and was marginally associated with C. trachomatis infection. S100A7 is secreted by epithelial cells in response to microbial challenge and expression can also be stimulated by IL-22 (Wolk et al., 2006). It is chemotactic for neutrophils and interferes with pathogen membrane permeability and enzyme function (Ganz, 2003; Guilhelmelli et al., 2013). The upregulation of S100A7 in various clinical stages of trachoma could support the loss of epithelial barrier function hypothesis described above and could contribute to chronic inflammation.

There was differential expression of several cell cycle genes (CD53, CDK13, CDC25C, PHYH, TYMS, TTK, TP53) in individuals with C. trachomatis infection and/or TP (Table 5). These expression changes could reflect regulation of the cell cycle in response to infection or inflammation, such as the proliferation or contraction of infiltrating inflammatory cells and cells effecting the wound healing response. Studies have reported that progression of the cell cycle is slow in cells infected with C. trachomatis (Kun et al., 2013; Elwell et al., 2016) and C. trachomatis is known to interfere with host cell apoptosis as a survival mechanism (Fischer et al., 2004; Ying et al., 2007; Siegl et al., 2014). Cluster 1 was enriched for genes in apoptosis pathways, however the combined expression of this cluster was not associated with clinical phenotypes or infection.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.