Immunopathogenesis of Progressive Scarring Trachoma: Results of a 4-Year Longitudinal Study in Tanzanian Children

Trachoma is initiated during childhood following repeated conjunctival infection with Chlamydia trachomatis, which causes a chronic inflammatory response in some individuals that leads to scarring and in-turning of the eyelids in later life. There is currently no treatment to halt the progression of scarring trachoma due to an incomplete understanding of disease pathogenesis. A cohort study was performed in northern Tanzania in 616 children aged 6 to 10 years at enrollment. Every 3 months for 4 years, children were examined for clinical signs of trachoma, and conjunctival swabs were collected for C. trachomatis detection and to analyze the expression of 46 immunofibrogenic genes.

T rachoma is a neglected tropical disease and the leading infectious cause of blindness worldwide. Disease is initiated by repeated infection of the conjunctival epithelium by the intracellular bacterium Chlamydia trachomatis, which occurs during childhood in areas where trachoma is endemic. Ocular C. trachomatis infection stimulates a chronic pathological inflammatory response in a proportion of exposed individuals, which leads to scarring of the conjunctiva. Conjunctival fibrosis tightens the eyelid, drawing inwards the lid margin (entropion) and eyelashes (trichiasis) such that they cause mechanical damage to the cornea. Without intervention (epilation of the lashes or surgery to correct lid orientation) this damage results in pain, secondary infection, and ultimately blindness. In 2017, 165 million people lived in districts requiring public health interventions for trachoma and 231,447 people worldwide were managed for trichiasis (1).
Mass antibiotic distribution with azithromycin is administered for population control of C. trachomatis in districts requiring intervention, alongside facial cleanliness promotion and environmental improvements to reduce transmission. However, scarring disease progresses in previously exposed individuals in the absence of evidence of ongoing C. trachomatis infection, suggesting that service provision for trichiasis management will be required for many years in districts where it was formerly endemic (2,3). There is currently no treatment to halt the progression of established scarring, partly due to an incomplete understanding of the immunopathophysiological process.
Following pathogen recognition and initiation of inflammation by the conjunctival epithelium, an adaptive immune cell-mediated response involving Th1 cells, classically activated macrophages, and gamma interferon (IFN-␥) is believed to be essential to clear ocular C. trachomatis infection (4)(5)(6). NK cells are an additional source of IFN-␥, and there is evidence that Th17 cells and associated cytokines are involved in the antichlamydial immune response (4,7,8). Following pathogen clearance, there is an extended period of inflammation and wound healing, which is thought to be characterized by the presence and activity of neutrophils and growth and matrix factors and by a reduction in the expression of mucin genes (4,8). These data have been gathered almost entirely from cross-sectional studies, and, as a result, the factors driving healthy wound healing versus pathological inflammation and fibrosis have not been differentiated.
We recently reported the results of a cohort study that investigated the association between conjunctival C. trachomatis infection and clinically visible inflammatory episodes with scarring progression in children over a 4-year period (9). The study took place from 2012 to 2016 in a region of northern Tanzania where trachoma is endemic. Scarring progression was observed in 103/448 (23%) individuals and was strongly associated with an increasing proportion of episodes with papillary inflammation (TP) (equivalent to P2 or P3 of the Detailed WHO Trachoma Grading System [2]). There were also marginal associations between "trachomatous inflammation-follicular" (TF) and C. trachomatis with scarring progression, which were shown to be mediated through TP. This suggests that other factors causing individual differences in TP contributed to scarring progression. However, the immune mediators driving TP and their cellular origins are unknown.
The aim of this study was to determine which components of the inflammatory response were most strongly associated with 4-year scarring progression and pathological TP. We further examined whether individuals with scarring progression responded differently at the gene expression level to C. trachomatis infection relative to nonprogressors.

RESULTS
A detailed description of the longitudinal study design and primary outcomes has been published elsewhere (4,9,10). Briefly, of the 666 eligible children, 616 were enrolled in the cohort study. Of these, 448 remained in the primary outcome analysis of factors associated with 4-year scarring progression. Scarring progression was observed in 23% (103/448) of individuals and was associated with increasing proportion of TP episodes. Gene expression was analyzed at study time points 1 to 5,7,9,11,13,15, and 17. The number of participants seen at each time point and the number in which C. trachomatis and TP were detected are shown in Table 1. In addition to the endogenous control genes HPRT1 and GAPDH, 46 genes of interest were quantified in all available samples at each of these time points. Three genes (FGF2, SERPINB4-SERPINB3, and IL22) and 61 observations were excluded from all analyses due to Ͼ10% missing data.
The association between gene expression and scarring progression was analyzed using random effects logistic regression models of the longitudinal data, clustering on participant identification number and adjusting for C. trachomatis infection, age, and sex. Of the 10 genes most strongly associated with 4-year scarring progression, only SPARCL1 (fold change [FC] ϭ 0.54; P ϭ 1.36 ϫ 10 Ϫ5 ) was downregulated in scarring progressors, whereas CXCL5, SOCS3, IL23A, IL19, CCL20, IDO1, MMP12, S100A7, and IL1B were upregulated ( Table 2). All upregulated genes had marginal fold changes of less than 1.4. With the exception of NCAM1, there was strong evidence that expression of all genes was altered with C. trachomatis infection ( Table 2). The three most strongly upregulated genes (FC Ͼ 5) in response to infection were IFNG, IL21, and CCL18, and the three most strongly downregulated genes (FC Ͻ 0.32) were MUC7, MUC5AC, and SPARCL1. There was evidence of an association between 20 genes and sex; 16 of these genes were upregulated in females. The top three genes most strongly associated with female sex were IL21, IL17A, and IFNG. Evidence was found for an association between age and all genes except ALOX5, CD247, IL12B, VIM, PDGFB, and TGFB1. Out of the genes associated with age, the majority (29/37) had negative fold changes, indicating that expression was higher in younger participants.
To further examine the genes most strongly associated with scarring progression, random effects lasso regression was performed in the longitudinal data set, clustering on participant identifier (ID). After filtering out incomplete cases, 3,622 observations remained in the analysis. In addition to genes associated with sex and age, the final iteration of the model retained 11 genes deemed to be most strongly associated with scarring progression, namely, SPARCL1, CXCL13, CCL18, CCL20, IL10, MMP12, IDO1, IL23A, S100A7, CXCL5, and IL19. All of these targets were also identified as strongly associated with scarring progression in the linear regression models (Table 2), whereas CXCL13, CCL18, and IL10 were additionally identified by the lasso regression only.
Prior analysis of the cohort data set indicated that the marginal association between C. trachomatis infection and 4-year scarring progression was mediated through TP, suggesting that an individual's response to infection-rather than simply the presence or absence of infection alone-determines their risk of scarring sequelae (9). In order to investigate this further, random effects regression models were performed for each gene in turn; these included an interaction term between scarring progression status and C. trachomatis infection in order to determine whether scarring progressors responded differently to infection relative to nonprogressors. Two genes, PDGFB and IL23A, showed some evidence of being upregulated in response to C. trachomatis infection by a greater amount in progressors than in nonprogressors (Table 3). In scarring progressors PDGFB was upregulated 1.58-fold in response to infection, whereas it was upregulated 1.31-fold in nonprogressors. Similarly, IL23A was upregulated 2.33fold in response to infection in scarring progressors and 1.77-fold in nonprogressors.
TP was identified as the major risk factor for scarring progression in our previous analysis (9); therefore, the analyses described above were repeated using TP as the primary outcome. In random effects linear regression models, all genes, with the exceptions of TGFB1 and MMP7, had evidence of an association with TP (Table S1); the three most strongly upregulated genes were CXCL13, CCL18, and S100A7, and the three most downregulated genes were SPARCL1, MUC7, and MUC5AC. Lasso regression was performed to identify the genes most strongly associated with TP (Table S2). The final iteration of the model only excluded three genes not deemed to be associated with TP, namely, MMP7, MUC1, and NCAM1. Age and C. trachomatis infection status but not sex

0.014
Psoriasin-1 (S100A7)  were included in the final model. The genes whose expression was most strongly associated with TP were MZB1, MMP9, CXCL5, PDGFB, TGFB1, IL17A, S100A4, IL8, S100A7, SPARCL1, and IDO1. In order to investigate whether individuals prone to TP (TP was detected at any time point) responded differently to C. trachomatis infection than individuals not prone to TP (in whom TP was never detected throughout the study duration), random effects linear regression models were repeated and included an interaction term between infection and whether any TP was detected. In individuals prone to TP, evidence was found that the expression of PDGFB, S100A7, IL23A, IL8, MMP9, CCL2, IL19, CCL18, and IL6 was upregulated, and the expression of MUC7, IDO1, S100A4, and MMP7 was downregulated in response to C. trachomatis infection, relative to individuals in which TP was never detected (Table S3A). Of these 13 genes, only S100A7 had a fold change greater than Ϯ1.5. Stratum-specific fold changes are shown in Table S3B.

DISCUSSION
In this study, we measured the expression of 46 immune response genes at 11 time points over a 4-year period, analyzing gene expression in relation to 4-year scarring progression status and pathological TP. The genes found to be associated with scarring progression included those encoding proinflammatory chemokines (CXCL5, CCL20, CXCL13, and CCL18), cytokines (IL23A, IL19, and IL1B), matrix modifiers (MMP12 and SPARCL1), immune regulators (IDO1, SOCS3, and IL10), and a proinflammatory antimicrobial peptide (S100A7). All genes except SPARCL1 were upregulated. A summary of the putative functions and interactions of these immune mediators is shown in Fig. 1. SPARCL1 had the greatest fold change and was the gene most strongly associated with scarring progression in this study. It was also strongly associated with TP and C. trachomatis infection. The downregulation of SPARCL1 has previously been reported to be associated with trachomatous scarring and inflammation (2,4,11). SPARCL1 is a nonstructural secreted protein that regulates the interaction between cells and the extracellular matrix (ECM). Its downregulation has been associated with cell proliferation and metastases in several cancers (12,13). Knockout of SPARCL1 in a murine corneal injury model led to the accumulation of inflammatory infiltrates, neovascularization and irregular ECM deposition at the site of injury (14). In contrast to wild-type mice, in which matrix metallopeptidase (MMP) activity was increased shortly after injury and then stabilized, MMP activity (detected through collagenase assay) was significantly greater in knockout mice, and activity increased throughout the duration of the experiment. This increased MMP activity was attributed to the excessive production of irregular collagen (14). MMP12, which is produced by monocytes, ocular epithelial cells, and fibroblasts following injury, was also significantly upregulated in scarring progressors in this study, and MMP9 was strongly associated with TP. In another murine corneal injury model, MMP12 was found to enhance early wound repair through increasing neutrophil infiltration and epithelial cell migration, suggesting that overexpression in trachoma might contribute to leukocyte infiltration (15). Histological analysis of conjunctival tissue from trichiasis patients has revealed inflammatory cell infiltrates and disrupted collagen structure consistent with the SPARCL1 knockout corneal injury model (14,16,17).
The chemokines CXCL5, CCL20, CXCL13, and CCL18, derived from epithelial and antigen-presenting cells on exposure to microbial stimuli and proinflammatory cytokines, recruit lymphocytes and neutrophils. IL-1␤ and S100A7 are innate proinflammatory mediators that are upregulated in the epithelium upon microbial exposure; S100A7 has antimicrobial and chemotactic properties (18). CXCL5, CCL18, IL1B, and S100A7 have consistently been found to be associated with trachomatous inflammation and scarring (2,8,11,(19)(20)(21)(22). The upregulation of these chemokines and proinflammatory mediators is indicative of ongoing inflammation in the conjunctival epithelium and the recruitment and activation of neutrophils and lymphocytes. The overexpression of these immune mediators in scarring progressors after adjusting for C. trachomatis infection could reflect differences in exposure to other microbes or to irritants such as dust or smoke (23,24). The strong association of S100A7, IL8, and CXCL5 with TP emphasizes the importance of innate epithelial responses in driving pathological inflammation.
The association of scarring progression with IDO1, SOCS3, and IL10, all of which are regulators of inflammation, is likely a result of the host's attempts to limit ongoing and pathological inflammation. While this may reflect that these individuals experienced greater inflammation than nonprogressors, it could also suggest that scarring progressors produce an excess of anti-inflammatory factors, leading to a poorer ability to clear infection. Genetic polymorphisms in interleukin 10 (IL-10) that lead to increased cytokine production and diminished cell-mediated immune responses to C. trachomatis have previously been reported (25,26). However, in this scenario, one might expect the expression of these immune regulators to be increased in response to C. trachomatis infection in scarring progressors relative to that in nonprogressors, which was not the case.
IL23A and IL19 were upregulated in individuals with scarring progression. IL23A was also upregulated in response to C. trachomatis infection in individuals with scarring progression and those prone to TP. IL-23A is a proinflammatory cytokine that is largely produced by dendritic cells and is essential for the survival and expansion of IL-17producing Th17 cells (27). In this study, IL17A was strongly associated with TP but not with scarring progression. Furthermore, IL21, which was strongly upregulated in response to C. trachomatis infection (and to a lesser extent to TP), is an autocrine factor that sustains Th17 cells (28). Psoriasis is an autoimmune disease characterized by excessive cytokine production that is triggered by environmental stimuli on a back-ground of genetic susceptibility (29). In psoriasis, keratinocytes recruit dendritic cells via CCL20, and the production of IL-23A by keratinocytes and dendritic cells recruits and activates IL-17A-producing Th17 cells, CD8 ϩ T cells, innate lymphoid cells (ILC), and ␥␦ T cells (30). IL-19 is produced by monocytes and epithelial cells and is also strongly upregulated in the keratinocytes of psoriasis patients (31,32). IL-19 has been proposed as a member of the inflammatory IL-23A/IL-17A cascade in psoriasis; its upregulation in keratinocytes is driven by IL-17A, and it acts in an autocrine fashion on keratinocytes to amplify the effects of IL-17A (31). One such effect is the induction of S100A7 production, which contributes to sustaining the cycle of inflammation (33). In the intestine, IL-23A is thought to act synergistically with IL-1␤ to promote pathogenic ILC and Th17 responses following infection with Helicobacter hepaticus (34). IL-23A-responsive ILCs were also found to mediate intestinal pathology in a murine colitis model (35), and IL-1␤-responsive IL-17A-producing ILCs have been reported in the conjunctiva (36). Together, these data suggest that the overexpression in scarring trachoma of CCL20, IL23A, IL19, IL1B, and S100A7 could reflect mechanisms of pathogenesis similar to those observed in psoriasis and intestinal inflammation, whereby the upregulation of IL23A and IL1B promotes the recruitment and activity of pathogenic IL-23A-responsive cells (Th17, ILC, or ␥␦ T cells, for example), which drive proinflammatory responses in the conjunctival epithelium, leading to chronic inflammation and fibrosis. The upregulation of IL23A in scarring progressors and TP-prone individuals in response to C. trachomatis infection could also suggest that greater polarization toward pathological IL-23Aresponsive cell types could be involved in initiating and sustaining pathological proinflammatory cycles in the epithelium. Evidence from murine models suggests that chlamydial antigens can be maintained in distal tissues such as the gut for long periods of time (37), which could offer an additional explanation for prolonged inflammation through the continued stimulation of circulating cells.
PDGFB was upregulated in response to C. trachomatis infection in scarring progressors relative to nonprogressors; however, it was not independently associated with scarring progression. PDGFB is a key mediator of the wound-healing process; it can promote the recruitment and activity of neutrophils, fibroblasts, and macrophages, and it induces the production of matrix molecules from fibroblasts and induces fibroblast and myofibroblast contraction (38,39). PDGFB could therefore directly contribute to scarring progression. PDGFB was upregulated in response to C. trachomatis infection and TP and was one of the genes most strongly associated with TP. The lack of direct association between scarring progression and PDGFB in a model adjusting for infection could indicate that its upregulation is tightly correlated with the presence of infection. The reason for its upregulation alongside IL23A in scarring progressors in response to C. trachomatis relative to expression in nonprogressors is unclear. A genome-wide association study of scarring trachoma suggested that host cell cycle, cell surface receptor signaling, and immune response pathways were associated with scarring, although no specific cytokine risk loci were identified (40).
The strengths of this study lie in the large sample size and the use of 11 longitudinal time points at which high-quality gene expression, infection, and clinical data are available. Of the genes associated with scarring progression, most had marginal fold changes of Ͻ1.4; however, small differences in multiple genes could act synergistically to favor pathological pathways (34). Although mucins MUC7 and MUC5AC were strongly downregulated in response to infection, consistent with our results at baseline (4), we found little evidence in this study for dysregulation of mucins (hypothesized to lead to loss of epithelial barrier function) being associated with scarring progression. Supporting previous results from ourselves and others (4)(5)(6), IFNG was strongly upregulated in response to infection but was not associated with scarring progression (and was not strongly associated with TP), suggesting that a Th1-cell or NK-cell IFN-␥ response is beneficial in the clearance of C. trachomatis infection. Upregulation of IL21 and IL17A was strongly associated with infection and female sex, suggesting the involvement of Th17 cells in the antichlamydial immune response and that Th17 cell responses might be heightened in females relative to those in males. Sex differences in immune responses have previously been reported, including an increase in Th17 responses in the intestinal mucosa of females (41,42). Given the potential role of the IL-23A/IL-17 axis in trachoma pathogenesis, this could offer a potential explanation as to why females are at greater risk of scarring trachoma than males (3). The majority of the genes tested were upregulated in younger participants, probably reflecting the decline in infection and disease incidence with increasing age. A caveat of this study is that gene expression data do not necessarily translate to changes in effector responses, and further research should aim to identify functional pathways of trachoma pathogenesis, including the phenotype and function of the cells responding to IL23A. We previously reported that mass azithromycin distribution (MDA) had an anti-inflammatory effect on conjunctival gene expression independent of the clearance of C. trachomatis infection, which was detectable by 3 (but not by 6) months posttreatment (10). However, analysis of the impact of MDA on gene expression in relation to scarring progression and whether it has a protective effect was outside the scope of this study.
Conclusions. Collectively, these data suggest that innate proinflammatory signals from the epithelium that drive leukocyte infiltration, IL-23A-responsive cells, and SPARCL1-, MMP12-, and PDGFB-mediated matrix reorganization and contraction are key pathways driving trachomatous scarring sequelae. The factors driving innate epithelial inflammation and causing scarring progressors to produce more IL23A and PDGFB in response to C. trachomatis infection remain unclear; they could have underlying genetic or epigenetic bases and/or be due to the presence of other organisms or irritants influencing the local immune response (43,44). However, despite several studies, thus far there is limited evidence for any major infectious or genetic risk factors (40,45). One could speculate that a complex combination of genetic or infectious risk factors increases an individual's susceptibility, such that upon the trigger of C. trachomatis infection, individuals possessing these risk factors develop a bias toward pathological IL-23A-responsive cells, which lead to sustained proinflammatory responses in the conjunctival epithelium that are exacerbated by other stimuli. The concept of infectious triggers causing sustained inflammation and fibrosis has been illustrated in the intestinal epithelium, with a number of distinct molecular pathways (46). Further research should seek to verify these pathways of trachoma immunopathogenesis at the functional level. Nevertheless, IL23A, SPARCL1, and PDGFB may be key mediators driving pathological inflammation and fibrosis in trachoma, and molecules that inhibit their action could hold therapeutic potential in preventing scarring progression.

MATERIALS AND METHODS
Ethical approval. This study was reviewed and approved by the Tanzanian National Institute for Medical Research, the Kilimanjaro Christian Medical Centre, and the London School of Hygiene & Tropical Medicine Ethics Committees and it adhered to the tenets of the Declaration of Helsinki. Written informed consent from a parent or legal guardian was requested from all study participants after detailed explanation in Swahili or Maa in the presence of a third person. A witnessed thumbprint was acceptable for consent if the individual was unable to read or write.
Study design and trachoma control. Study participants were recruited from three rural and predominantly Maasai villages in northern Tanzania. The study design and population have been described in detail in several earlier reports (4,9,10). In brief, a cohort of 616 children, aged 6 to 10 years at the beginning of the study in February 2012, were enrolled in the longitudinal cohort study and were visited every 3 months for 4 years, for a total of 17 time points.
The SAFE strategy (surgery for trichiasis [in-turned eyelashes], antibiotics, facial cleanliness, and environmental improvement) was implemented in the study villages by the field team in collaboration with district eye coordinators, following approval from the Tanzanian Ministry of Health (MoH). Education about environmental improvements and facial hygiene was provided by the field team, free trichiasis surgery was offered, and all members of the three villages (including study participants) were offered azithromycin for trachoma control during the August of the years 2012, 2013, and 2014.
Clinical examination and sample collection. At each time point, all available and consenting cohort participants were examined. Eyes were initially anaesthetized using preservative-free proxymetacaine hydrochloride 0.5% eyedrops. Each participant's left eyelid was everted, and the tarsal conjunctiva was examined by an ophthalmic nurse experienced in trachoma grading using 2.5ϫ loupes and a torch. Clinical signs were graded using the WHO detailed FPC grading system (47). Using the FPC grading system, F2/F3 corresponds to "trachomatous inflammation-follicular" (TF) and P3 corresponds to "trachomatous inflammation-intense" (TI) of the WHO simplified grading system (48). We consider P2 to also represent significant clinically apparent inflammation, and we therefore refer to P2/P3 as "TP" and use this designation in all analyses (2,9,10). At each time point, high-resolution conjunctival photographs were taken with a Nikon D90 camera with a 105-mm macro lens.
At baseline (time point 1, or time point 2 if not seen at time point 1) and final time points (time point  17, or time point 16 if not seen at time point 17), conjunctival photographs were independently graded by an ophthalmologist using a detailed scarring grading system (49). Baseline and final photographs for each individual were subsequently compared side by side in order to determine whether there was no progression (no scarring, or no progression of existing scarring) or incidence/progression of trachomatous scarring. Participant entry into the longitudinal study was permitted at time points 1 and 2, after which no new participants were enrolled.
At each time point swab samples were collected from the left tarsal conjunctiva using sterile polyester-tipped swabs (Puritan), as described previously (4). The first swab was collected into 250 l RNAlater (Invitrogen), and the second was stored dry. Swabs were stored on ice in the field. RNAlater swabs were stored at 2 to 8°C overnight before transfer to Ϫ80°C for long-term storage. Dry swabs were stored immediately at Ϫ80°C.
Chlamydia trachomatis detection. At time point 1, DNA was extracted from dry-stored swabs using a PowerSoil DNA isolation kit (Mo Bio Laboratories), and this DNA was used for C. trachomatis detection by droplet digital PCR (ddPCR), as described previously (4). At time points 2 through 17, RNA and DNA were extracted from RNAlater-stored swabs using RNA/DNA purification kits (Norgen Biotek) following the manufacturer's instructions, and this DNA was used for C. trachomatis detection by quantitative PCR (qPCR) (9, 10). Triplex qPCR was performed targeting chlamydial chromosomal (omcB) and plasmid (pORF2) targets and a human endogenous control gene (RPP30) (50). Samples were tested in duplicate and were determined positive if RPP30 in combination with pORF2 and/or omcB targets amplified for Ͻ40 cycles in one or both replicates. Norgen-extracted DNA from time point 2 was tested by qPCR and ddPCR for the C. trachomatis plasmid target for comparison, and the kappa score for agreement was 0.84 (10).
Gene expression. Norgen-extracted RNA from time points 1 to 5,7,9,11,13,15, and 17 was reverse transcribed using SuperScript VILO cDNA synthesis kits (Thermo Fisher Scientific) following the manufacturer's instructions. The relative abundances of 48 genes of interest, including those of GAPDH and HPRT1 endogenous control genes, were quantified in each sample by qPCR using TaqMan microfluidic 384-well array cards (Thermo Fisher Scientific). qPCR was performed on a ViiA7 thermal cycler with TaqMan Universal mastermix, as described previously (4). The 46 genes of interest were shortlisted from a total of 91 genes tested at time point 1 (4), based on those most strongly associated with C. trachomatis infection and/or clinical signs. The original 91 genes were selected based on the results of previous cross-sectional studies and were centered around key biological processes hypothesized to underlie the immunopathogenesis of trachoma, including antimicrobial peptides, cell cycle regulators, cytokines and chemokines, biomarkers of epithelial-mesenchymal transition, matrix modifiers, the response to microbiota, mucins, NK cell markers, pattern recognition receptors, and signaling pathway regulators.
Analysis. Data were stored in Microsoft Access and were analyzed in STATA v15 and R (www.R -project.org). Gene expression data were normalized using the cycle threshold (ΔC T ) method (51), normalizing the expression of each gene to the expression of HPRT1 in the same sample. For quality control purposes, both genes and observations (defined as a sample from a participant at one time point) with Ͼ10% missing data were excluded. This resulted in the exclusion of three genes (FGF2, SERPINB3-SERPINB4, and IL22) and 61/4,853 observations, leaving 4,792 observations and 43 genes of interest (excluding HPRT1 and GAPDH) in the analysis.
In order to assess the association between longitudinal gene expression and scarring progression, random effects linear regression models were performed for each gene in turn, using gene expression as the dependent variable and adjusting for infection, age, and sex, clustering on participant ID. These analyses were subsequently repeated for TP, adjusting for infection, age, and sex.
Random effects lasso regression models using the glmmLasso package in R (52) were used to select the genes most strongly associated with scarring progression and TP. This analysis does not permit any missing data, and therefore only complete cases were retained from the filtered data set described above. This resulted in a data set of 3,622 observations and 43 genes for the scarring progression analysis and a data set of 4,510 observations and 43 genes of interest for the TP analysis. Age, sex, and infection were adjusted for in each model.
In order to determine whether scarring progressors responded differently from nonprogressors to C. trachomatis infection at the gene expression level, the random effects linear regression analyses described above were repeated, including an interaction term between scarring progression status and C. trachomatis infection. These analyses were further repeated to investigate whether individuals prone to TP (defined as those who had TP at any time point) responded differently to infection relative to individuals in whom TP was never detected throughout the study duration. To provide an objective threshold for reporting the genes most strongly associated with the outcome, we highlight those genes with a P value below the level that controls the false discovery rate (FDR) to be less than 5%, using the Benjamini-Hochberg method (53).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, XLSX file, 0.03 MB.