The conjunctival transcriptome in Ethiopians after trichiasis surgery: associations with the development of eyelid contour abnormalities and the effect of oral doxycycline treatment

Background: Surgery to correct trichiasis is a key component of the World Health Organisation trachoma control strategy, however unfavourable outcomes such as eyelid contour abnormalities (ECA) following surgery are relatively common. This study aimed to understand the transcriptional changes associated with the early development of ECA and the impact of doxycycline, which has anti-inflammatory and anti-fibrotic properties, upon these transcription patterns. Methods: One thousand Ethiopians undergoing trichiasis surgery were enrolled in a randomised controlled trial following informed consent. Equal groups of randomly assigned individuals were orally administered with 100mg/day of doxycycline (n=499) or placebo (n=501) for 28 days. Conjunctival swabs were collected immediately prior to surgery and at one- and six-months post-surgery. 3’ mRNA sequencing was performed on paired baseline and one-month samples from 48 individuals; 12 in each treatment/outcome group (Placebo-Good outcome, Placebo-Poor outcome, Doxycycline-Good outcome, Doxycycline-Poor outcome). qPCR validation was then performed for 46 genes of interest in 145 individuals who developed ECA at one month and 145 matched controls, using samples from baseline, one and six months. Results: All treatment/outcome groups upregulated genes associated with wound healing pathways at one month relative to baseline, however no individual differences were detected between groups. The summed expression of a highly coexpressed cluster of pro-fibrotic genes was higher in patients that developed ECA in the placebo group relative to controls. qPCR validation revealed that all genes in this cluster and a number of other pro-inflammatory genes were strongly associated with ECA, however these associations were not modulated by trial arm. Conclusions: The development of post-operative ECA is associated with overexpression of pro-inflammatory and pro-fibrotic genes including growth factors, matrix metalloproteinases, collagens and extracellular matrix proteins. There was no evidence that doxycycline modulated the association between gene expression and ECA.


Introduction
Trachoma is a neglected tropical disease and the world's leading infectious cause of blindness. Disease is initiated during childhood by repeated infection of the conjunctival epithelium with the intracellular bacterium Chlamydia trachomatis. This can lead to recurrent episodes of pathological inflammation which persist throughout the lives of some individuals and are a major risk factor for the development of scarring 1 . Scarring of the tarsal conjunctiva causes the lid margin and eyelashes to turn inwards, known as entropion and trichiasis, respectively, such that the lashes brush against the cornea causing mechanical damage, pain and ultimately blindness. In 2017 over 200,000 people worldwide were managed for trichiasis; the disease is found predominantly in sub-Saharan Africa 2 . Trachomatous inflammation and scarring can progress in the absence of detectable chlamydial infection 1,3 , therefore it is expected that trichiasis surveillance and management provisions will be required for many years in formerly endemic districts.
Trachoma is controlled at the population level through implementation of the SAFE Strategy: Surgery to correct trichiasis, Antibiotics for infection control, and Facial cleanliness and Environmental improvements to prevent infection transmission. However, reported post-operative trichiasis rates following trichiasis surgery range between ~7% at one year to 62% at three years 4, 5 and unfavourable outcomes such as the development of eyelid contour abnormalities (ECA) following surgery are relatively common, occurring in around 10% of patients 6,7 .
ECA involves notching or distortion of the eyelid margin. When ECA is severe the eyelid may no longer cover the globe of the eye when closed (lagophthalmos), increasing the risk of corneal damage and vision loss. Unfavourable aesthetic outcomes such as ECA may also dissuade patient's friends or relatives from undergoing trichiasis surgery and may be stigmatising [8][9][10] . The pathophysiology of ECA is unknown, although it is likely due to an aberrant post-operative wound healing response involving excessive fibrosis and contraction of the tarsal conjunctiva. Indeed, the immediate post-operative appearance is associated with both the return of trichiasis and the development of eyelid contour abnormalities 11 . Previous studies have shown that a number of innate epithelial pro-inflammatory mediators (IL1B, S100A7, CXCL5), matrix factors (SPARCL1, CTGF, collagens) and matrix metalloproteinases ( MMPs 7,9,12) are associated with progressive trachomatous scarring, trichiasis and post-operative trichiasis 1,12,13 . MMPs are thought to facilitate fibrosis through remodelling the extracellular matrix (ECM) and enabling inflammatory cell infiltration. Failure of glaucoma filtration surgery is associated with excessive ECM deposition and contraction of subconjunctival tissue 14 , and a number of collagen genes were found to be upregulated in the conjunctivae of patients with fibrotic outcomes following glaucoma surgery 15 . These pro-inflammatory and matrix factors may also have a role in the development of ECA.
Doxycycline is a tetracycline-based antibiotic that is used clinically to treat sexually transmitted C. trachomatis infection. Several reports have presented evidence that doxycycline reduces inflammatory gene expression [16][17][18][19] and it inhibits MMPs by reducing their expression, blocking their activation and inhibiting protease activity through directly binding zinc and calcium ions in the catalytic site 20 . Doxycycline is used clinically to reduce inflammation and prevent fibrosis in diseases including rosacea, idiopathic pulmonary fibrosis and corneal erosion 21-23 . Li and colleagues presented data showing that doxycycline inhibited matrix remodelling and contraction and reduced MMP expression in primary fibroblasts isolated from individuals undergoing trichiasis surgery 24 . These data therefore suggest that doxycycline could have a protective effect against the development of post-operative trichiasis and ECA.
In order to address this question a randomised, double-blind, placebo-controlled trial was conducted to determine the impact of doxycycline treatment on post-operative trichiasis. One thousand Ethiopians eligible for trichiasis surgery were enrolled in the trial and the findings have been published in full elsewhere 7 . Surprisingly, the cumulative proportion of individuals who developed post-operative trichiasis by 12 months after surgery was the same in both groups; 58/498 (12%) in the doxycycline group and 62/501 (12%) in the placebo group. There was also no difference in the proportion of individuals with mild, moderate or severe ECA at one, six or twelve months (number with mild and moderate/severe ECA at one month (placebo: 179/500 (35.8%), doxycycline: 181/492 (36.8%)), six months (placebo: 82/493 (16.6%), doxycycline: 88/481 (18.3%)) and twelve months (placebo: 55/494 (11.1%), doxycycline: 66/488 (13.5%)).
This study sought to understand the transcriptional changes associated with the early development of ECA following

Amendments from Version 1
The introduction has been updated to discuss surgical appearance and subsequent ECAs. We have added a limitation of the study in not recording digital photographs within one day post-operation and clarified that samples that failed QC for gene expression did not differ from those that passed QC by age, sex, or group. We included citations to important work we had previously omitted and corrected other citations and text that were inaccurate. We have added additional information of the prior experience of surgeons. Noted that full trial details are published elsewhere and stated this in the methods. We have justified the study design as embedded in a clinical trial and explained why additional post-hoc power/samplesize calculations were not carried out. We have updated the discussion to address how the results impact understanding of the pathophysiology of ECA development and the future development of targeted anti-inflammatory or anti-fibrotic therapies.
Any further responses from the reviewers can be found at the end of the article REVISED trichiasis surgery and the impact of doxycycline upon these transcription patterns. It was nested within the overall clinical trial and took place in parallel.

Ethics
This study was approved by the ethics committee of the London School of Hygiene & Tropical Medicine, Emory University Institutional Review Board, the Ethiopian National Health Research Ethics Review Committee and the Ethiopia Food, Medicine and Healthcare Administration and Controls Authority. The study was conducted in accordance with the Declaration of Helsinki and the International Conference on Harmonisation-Good Clinical Practice guidelines. Written informed consent or a witnessed thumbprint was requested from all study participants following detailed explanation of the study in Amharic.

Surgery and examination
One thousand Ethiopians eligible for trichiasis surgery were enrolled in a randomised, double blind, placebo-controlled trial (registered with the Pan African Clinical Trials Registry, trial number PACTR201512001370307, registration date 29 th November 2015, registration URL: https://pactr.samrc.ac.za/TrialDisplay.aspx?TrialID=1370). Eligible individuals had evidence of trichiasis (one or more lashes touching the eye or evidence of epilation) and trachomatous scarring of the tarsal conjunctiva, were aged ≥ 18 years and had not had previous trichiasis surgery. Trichiasis surgery was performed using the posterior lamellar tarsal rotation (PLTR) procedure by experienced nurse-surgeons who had recently undergone refresher training. Six surgeons were standardised, prior to operating in the trial. One operated on only 4 eyes (moved to another region). The remaining five operated throughout the trial. These 6 surgeons had been performing TT surgery for a median of 4 years (range 1 -9); they had performed a median of 752 surgery (range 100 -1500). All study participants received 1% tetracycline eye ointment for self-administration twice a day for two weeks post-surgery. Participants were randomly assigned (1:1) to receive either 100mg/day doxycycline or 100mg/day placebo, orally administered, for 28 consecutive days after surgery. Full details of trial design and methodology, covering inclusion/exclusion criteria, are published in the main trial paper 7 . The clinical trial was powered to detect a difference in the primary outcome of 12-month cumulative postoperative trichiasis. Sample size/power was based on a previous randomised trial in Ethiopia 25 . The a priori design of this nested study was the assessment of gene expression profiles in all people developing ECA (and a matched control group). As the numbers that develop ECA were not known beforehand a separate sample size calculation for this outcome was not performed and a post-hoc power calculation is not recommended.
Study participants were examined at baseline (immediately prior to surgery), 10 days, one month, six months and twelve months post-surgery. Examinations were performed by two standardised examiners with very strong agreement on primary outcome (k=0.92) (at baseline and twelve months (EH) and at one and 6 months), using 2.5X loupes and a torch. Clinical signs of trachoma (follicles, papillae and cicatrice; FPC) were graded using the detailed WHO grading system 26 . Clinical signs or evidence of trichiasis, epilation and corneal scarring were recorded as previously described 7 . At baseline, scarring was graded using the detailed FPC grading system, whereas after surgery scarring was graded using the post-operative scarring grading system 25 . Sutures were removed on day 10. The primary outcome of the overall trial was the cumulative proportion of individuals who developed post-operative trichiasis by 12 months. Post-operative trichiasis was defined as minor (1-5 eyelashes touching the eye or evidence of epilation in less than a third of the eyelid margin) or major (≥6 eyelashes touching the eye or evidence of epilation in a third or more of the eyelid margin). ECA at one, six or 12 months was included as a secondary trial outcome. ECA was graded as previously described 6,27 : mild (vertical deviation from the natural contour less than 1 mm in height and affecting more than one third of horizontal eyelid length); moderate (vertical deviation from the natural contour 1-2 mm in height or affecting one third to two thirds of horizontal eyelid length) and severe (vertical deviation from the natural contour more than 2 mm in height or a defect more than two thirds the horizontal eyelid length). Mild ECA is considered clinically non-significant whereas moderate and severe ECA are considered clinically significant. Drug adherence was monitored through interviews and by counting the number of capsules remaining. Full trial details and results are published elsewhere 7 .

Sample collection
Swabs from the tarsal conjunctiva were collected at baseline and at one, six-and twelve-months post-surgery. A sterile polyester swab (Puritan) was passed across the conjunctival surface four times, with a quarter-turn on each pass. Swabs were stored in RNAlater (Thermo Fisher Scientific) at 4°C overnight and then long-term at -80°C. Swabs were shipped (on dry ice) to Kilimanjaro Christian Medical Centre (KCMC), Tanzania for processing. RNA and DNA were extracted from swabs using Norgen RNA/DNA purification kit (Norgen Biotek) following manufacturer's instructions.
Chlamydia trachomatis testing DNA samples were tested for C. trachomatis using a published qPCR assay 28 . Samples were tested in duplicate using the triplex assay which detects a human endogenous control gene (RPP30), C. trachomatis plasmid (pORF2) and chromosomal (omcB) targets. Each reaction of 20μl contained 10μl TaqMan Multiplex Master Mix (Thermo Fisher Scientific), 300nM of each primer and probe (Thermo Fisher Scientific), 4μl water and 4μl sample DNA. Samples were tested on a ViiA7 thermal cycler (Thermo Fisher Scientific) with the following conditions: a 20 second hold at 95°C, followed by 40 cycles of 95°C for one second and 60°C for 20 seconds. Samples were determined C. trachomatis positive if RPP30 and both chlamydial targets (pORF2 and omcB) were positive in either or both replicates. Target concentration was calculated by extrapolating from a standard curve.
3'RNA sequencing 3'RNA sequencing was performed on 96 samples (paired baseline and one-month) from 48 individuals; 12 in each treatment/outcome group (Placebo-Good outcome, Placebo-Poor outcome, Doxycycline-Good outcome, Doxycycline-Poor outcome). Based on an average read count of differentially expressed genes (50), dispersion estimate (0.2) the total number of genes for testing (11000) data, and the number of differentially expressed genes between baseline and 1-month (1000), we expect 12-paired samples will have 88% power to detect a fold change of 2.2 at a FDR threshold of 10% or 90% power at a FDR threshold of 5%. All 48 individuals had moderate trachomatous scarring and no eyelid distortion at baseline (Detailed FPC scarring grade 2). Individuals with Poor outcomes had scarring distortion (SC4 or SC5) at one and six months, grade 2 or 3 ECA at one month, and grade 1, 2 or 3 ECA at six months. Individuals with Good outcomes had no scarring distortion or ECA at one and six months. Potential participants were selected at random (randomisation was performed in STATA v15) from each of the four treatment/outcome groups.
Library preparation was performed using QuantSeq 3' mRNA-Seq Library Prep Kit FWD for Illumina (Lexogen GmbH, Austria), following the manufacturer's instructions. 3' RNA sequencing generates one read per transcript at the 3' end of polyadenylated mRNA and maintains strand-specificity. Sample RNA concentration was measured using Qubit HS RNA kits (Thermo Fisher Scientific) and RNA quality was tested using Agilent RNA 6000 Pico chips (Agilent Technologies, USA). Sample concentrations were normalised to 10ng/μl and 5μl RNA was used for library preparation. The PCR Add-on kit (Lexogen GmbH) was used during sample preparation to generate libraries that were not under or over-cycled. Final library concentration was tested using Agilent DNA chips and libraries were normalised to 6nM prior to pooling. Libraries were shipped to Lexogen GmbH for quality control and sequencing on two lanes of an Illumina HiSeq (Illumina).
Bam files were converted to FastQ using Picard and FastQ files were uploaded to BlueBee Genomics Platform (BlueBee). FastQ files were trimmed and quality filtered using the default Integrated QuantSeq FWD workflow. In brief, Bbduk was used to trim low-quality tails, adapter contamination and polyA read-through, reads were aligned to human genome version GRCh38.77 using STAR Aligner and reads were counted using HTSeq-count with kit-specific options. Raw and processed RNAseq data are deposited within the NCBI GEO public database (accession number GSE135455).

Gene expression
Of the 1000 participants in the clinical trial, a total of 145 had developed moderate or severe ECA (grade 2 or 3) at the one-month time-point 7 . These 145 individuals and 145 age, sex and trial arm matched individuals without ECA (none or mild ECA (grade 0 or 1)) at one month were selected for qPCR validation of gene expression, using samples from baseline, one month and six months; making a total of 868 samples. For TLDA qPCR assays 46 genes and 2 controls were selected for testing. If the desired minimum fold change is 1.5, we will need to study 80 subjects in contrast at probability (power) 0.8 using exact test and an FDR of 1%.
RNA was converted to cDNA using SuperScript VILO cDNA synthesis kits (Thermo Fisher Scientific) following the manufacturer's instructions. The relative abundance of 46 genes of interest and endogenous control genes HPRT1 and GAPDH were quantified in each sample by qPCR using TaqMan Microfluidic 384-well Array Cards (Thermo Fisher Scientific). qPCR was performed using TaqMan Universal Master Mix (Thermo Fisher Scientific) on a ViiA7 thermal cycler following the manufacturer's instructions. Genes of interest were selected based on 3'RNA sequencing results, those hypothesised to be involved in the pathogenesis of scarring trachoma and genes reported to be immunomodulated by doxycycline 1,12,[16][17][18][19]24,29,30 .
After running the first 239 samples and noting poor amplification or failure of a number of samples on TaqMan Array Cards, all remaining samples were tested for HPRT1 expression by single-plex qPCR for quality screening. qPCR was performed using a HPRT1 TaqMan Assay (the same as that used on the Array Card (assay ID: Hs02800695_m1; Thermo Fisher Scientific) and TaqMan Universal Master Mix on a ViiA7 thermal cycler following the manufacturer's instructions. A cycle threshold value of 31 was determined as the screening cut-off based on data from Array Cards that had already been run. Of the 631 remaining samples that had not yet been run on Array Cards, 63 had HPRT1 qPCR CT values >31. After excluding failed samples and those that did not pass HPRT1 screening, 774/868 samples remained with gene expression data. Samples that were excluded did not differ from the remaining samples by treatment/outcome group (Pearson Chi2 P value = 0.591), age (mean age of included: 55.9; excluded 56.5, t test P value = 0.665) or sex (included: 587/774 female (75.6%); excluded: 72/94 female (76.6%), Pearson Chi2 P value = 0.162).

Data analysis 3'RNA sequencing analysis.
Read count data were downloaded from BlueBee and imported to R (www.R-project.org, R version 3.5.1) for analysis. Read count data were analysed for differential expression using DEseq2, adjusting for sex and age 31 . Genes were retained in the analysis that had a read count of 10 or more in 48 or more samples (n=11424). The DEseq2 package internally corrects for library size and a plot of the Cook's distance between samples revealed no sample outliers (data not shown) therefore no samples were excluded from the analyses. Samples were assigned to four groups; doxycycline and good outcome, doxycycline and ECA, placebo and good outcome, placebo and ECA. The differences in gene expression at one month relative to baseline in paired samples from each group was compared initially. The differences in gene expression between baseline and one month within each group were then contrasted between groups in pairwise comparisons. P values were adjusted by the Benjamini and Hochberg method within the DEseq2 package. Principal component analysis was performed within DEseq2 using read count data. The matrix of raw read counts for all 96 samples and 11424 genes was entered into Miru v1.0 (https://kajeka.com; now called Graphia) to investigate clusters of samples with highly correlated expression values. The Pearson correlation cutoff was increased to the maximum possible whilst still retaining all 96 samples (cutoff = 0.93).
Genes with an adjusted P value (Padj) <0.1 and FC>1.5 within each group at one month relative to baseline were entered into Gene Set Over-representation analysis in Consensus Path Database (CPDB, http://consensuspathdb.org), using all 11424 genes as background, minimum overlap with input list=2 and P value cutoff=0.01.
The matrix of raw read counts for all 96 samples and 11424 genes was transposed and entered again into Miru to investigate networks of co-expressed genes. Genes were clustered by Pearson correlation coefficient using a cutoff >=0.88. The genes in each cluster were also entered into Gene Set Overrepresentation analysis in CPDB using 11424 genes as background. For each individual, the read counts for all genes in each cluster were summed, creating a single summative value per cluster. These summative values were then contrasted between one-month post-surgery treatment/outcome groups using linear regression, using the placebo/good outcome group as the reference group and adjusting for baseline (pre-surgery) expression, age and sex.
qPCR. Gene expression results were imported into STATA v15 for analysis. Cycle threshold values in each sample were normalised to the endogenous control gene (HPRT1) to generate delta CT values 32 . For quality control purposes, both genes and observations (a sample from a participant at one time-point) with 10% missing data were excluded. This resulted in the exclusion of two genes, PAPPA and TRGV9, and 22/772 observations, leaving 46 genes (44 genes of interest in addition to HPRT1 and GAPDH) and 751 observations in the analysis.
A binary ECA variable was created for analysis purposes based on the clinical significance criteria at each of the one-and six-month time-points: none or mild ECA (grades 0 or 1) = 0 and moderate and severe ECA (grades 2 or 3) = 1. In order to determine whether any of the shortlisted genes were associated with ECA, mixed effects logistic regression models were constructed with binary ECA as the dependant variable and the delta CT values of each gene in turn as the independent variable, adjusting for trial arm, age and sex, and using data from baseline, one and six months. Participant identification number was included as a random effect to account for samples from the same individuals being used from multiple timepoints. C. trachomatis was not adjusted for in these analyses due to the very low number of infections detected. Following these initial models, the models were performed again including an interaction term between gene expression and treatment arm, to determine whether treatment arm modified the effect of gene expression on ECA outcome. A likelihood ratio test was then performed to determine whether the model including the interaction explained the data better than the model without the interaction. The Benjamini and Hochberg method was used to control for a false discovery rate of 5% 33 .

RNA sequencing
3'RNA sequencing was performed on samples from 48 individuals at baseline (immediately prior to surgery) and at one-month post-surgery, making a total of 96 samples. All 48 individuals had moderate trachomatous scarring and no eyelid distortion at baseline (FPC scarring grade 2). Of these 48 individuals, there were 12 in each treatment/outcome group: Placebo-Good outcome, Placebo-Poor outcome, Doxycycline-Good outcome, Doxycycline-Poor outcome. The demographic and clinical characteristics of these 48 individuals are shown in Table 1. There was no difference in sex between groups (Chi2 P-value=0.763), however there was a slight difference in age between groups (ANOVA P-value = 0.0143).
After filtering out genes that had a read count <10 in 48 or more samples, 11,424 genes were retained. RNA sequencing analysis revealed that within each of the four groups (Placebo-Good; Placebo-Poor; Doxycycline-Good; Doxycycline-Poor), a large number of genes were differentially regulated at one month relative to baseline (Table 2). In each group, roughly twice as many genes were upregulated as were downregulated. More genes were differentially regulated in the Poor groups relative to the Good groups in both placebo and doxycycline treatment arms. Whilst the number of differentially expressed genes was similar in both Poor groups, there were a higher number of differentially expressed genes in the doxycycline-Good group relative to the placebo-Good group.
Differentially expressed genes (Padj<0.1, FC>1.5) within each treatment/outcome group at one month relative to baseline were entered into pathway analysis. The ten most enriched pathways for each group are shown in Table 3. All pathways were reflective of wound healing responses, including ECM organisation, collagen formation and interactions between cells and the ECM. The two Poor outcome groups were enriched for several immune response pathways that were not present in the top ten enriched pathways of the Good outcome groups ( Table 3).
The differences in gene expression at one month relative to baseline within each treatment/outcome group were then contrasted between groups. Pairwise comparisons revealed almost no significant differences in gene expression between groups (Table 4). Four genes (PSD3, ZEB1, RAD17, RNF181) were differentially expressed in the placebo group, and one (NUAK2) in the doxycycline group, between Good and Poor outcomes (Padj<0.1, FC>1.5). There were no overall differences between placebo and doxycycline groups (irrespective of outcome). Four genes were differentially expressed between Good and Poor outcomes (irrespective of treatment group); two known genes (PTP4A1, FBXO31) and two uncharacterised genes (ENSG00000256940, ENSG00000257742).
Principal component analysis of read count data revealed considerable within-group variation in the one-month samples ( Figure 1). The raw read count matrix for all 96 samples and  11424 genes was entered into coexpression analysis using a Pearson correlation cut-off of 0.93. This generated three clusters of samples (data not shown). Nine samples were not included in any cluster and had no class. Treatment/outcome group samples at baseline and one month appeared to distribute evenly across all three clusters, with the exception of the placebo-Poor one-month samples, which were all found in cluster one (with the exception of one sample that had no class).
The raw read count matrix for all samples was inverted and analysed again for networks of co-expressed genes. Using a Pearson correlation cut-off of 0.88, 40 clusters of genes were identified. Genes in each of the ten largest clusters (containing ≥9 genes) were entered into pathway analysis using the background of 11424 genes included in differential expression testing. The largest cluster, cluster 1, was highly enriched for pathways associated with wound healing (Extended data).
With the exception of cluster 8, which was enriched for phagocytosis-related genes, all other clusters were enriched for cell cycle and pathways of general cellular processes (Extended data). In order to determine whether these tightly co-expressed gene clusters were differentially expressed between treatment/outcome groups, the read counts for all genes in each cluster were summed for each individual. The expression of each cluster-sum was then compared between groups at one month by linear regression, adjusting for cluster-sum expression at baseline (prior to surgery), age and sex.
The summative expression of the 20 genes in cluster one, which was strongly enriched for pro-fibrotic pathways, was higher in those with a Poor outcome relative to those with a Good outcome in the placebo group (P=0.017). The summative expression of cluster one genes was not different in the doxycycline group in those with Good or Poor outcomes relative to the placebo-Good group however. There was no evidence of differential expression of clusters 2-10 between treatment/ outcome groups.

Gene expression
The genes in cluster 1, which were enriched for pro-fibrotic pathways and were more highly expressed in individuals who developed a Poor outcome relative to Good outcome in the placebo arm, were tested by qPCR in a larger sample of individuals to validate their differential expression. A number  Two genes were filtered out prior to analysis due to having >10% missing data (PAPPA and TRGV9), leaving 44 genes of interest in the analysis. Twenty-six of these genes were strongly associated with ECA (at the FDR adjusted P value threshold of 0.0036), all of which were upregulated, reflected by an odds ratio less than one (delta CT values were used in the analysis) ( Table 6).
Of the 20 genes in cluster 1, 17 were tested by qPCR (three collagen alpha variants were excluded; COL1A2, COL5A2, COL6A3). Furthermore PAPPA, which was in cluster 1 and was tested by qPCR, was excluded during quality control. All 17 genes from cluster 1 that were tested by qPCR were strongly associated with ECA (P<0.0036). Of the genes that were selected for qPCR testing based on prior literature (total = 29), 10 were associated with ECA. These included proinflammatory chemokines and cytokines CCL2, CCL3, CXCL5, IL1B and IL6, matrix metalloproteinases MMP1 and MMP9, CHUK, FUT4 and PTPRC.
In order to determine whether any of the associations between gene expression and ECA were modulated by doxycycline, regression models were run again using data from baseline and one month only (since treatment was only administered for 28 days after surgery). There was no strong evidence that any models including the interaction fitted the data better than the models without the interaction term ( Table 7), suggesting that doxycycline treatment did not modulate the association between conjunctival gene expression and ECA. There was weak evidence for interactions in CD68, IL10, IFNG and MMP12.

Discussion
Here we present evidence that a number of pro-fibrotic and pro-inflammatory genes were upregulated in the conjunctivae of patients who developed ECA following trichiasis surgery. Analysis of 3'RNA sequencing data revealed that all treatment/outcome groups upregulated wound healing pathways at one-month post-surgery, however surprisingly, there were no differences in the expression of individual genes between groups. Visualisation of global gene expression by PCA revealed considerable variation amongst one-month samples relative to baseline samples. Clusters of highly co-expressed genes were identified in the raw data, the largest of which was enriched for pro-fibrotic pathways. The summed expression of this cluster of co-expressed genes was upregulated in individuals with a Poor outcome relative to those with a Good outcome in the placebo group. However, this cluster was not upregulated in the doxycycline-treated Poor outcome group. The expression of the genes in this cluster and a number of other immune response genes hypothesised to be associated with ECA and reported to be immunomodulated by doxycycline were evaluated in a larger number of participants by qPCR at baseline, one and six months. The results showed that all genes tested within this pro-fibrotic cluster and a number of the other pro-inflammatory genes were strongly associated with ECA, adjusting for trial arm, age and sex. No strong evidence was found to suggest that doxycycline treatment had an immunomodulatory effect on the association between gene expression and ECA.
The genes strongly associated with ECA included all tested members of cluster one (enriched for collagen biosynthesis/formation, integrin interactions and extracellular matrix organisation pathways), pro-inflammatory cytokines and chemokines (CCL2, CCL3, CXCL5, IL1B and IL6), matrix metalloproteinases MMP1 and MMP9, CHUK (also known as IKK-α, an inhibitory regulator of NFκB transcription factor), FUT4 (also known as CD15, a cellular marker of granulocytes), and PTPRC (also known as CD45 or common leukocyte antigen). All associated genes were upregulated in individuals with ECA.
The upregulation of pro-inflammatory mediators and cell markers (CCL2, CCL3, CXCL5, IL1B, IL6, CHUK, FUT4 and PTPRC) is indicative of higher levels of inflammation and leukocyte infiltration in the conjunctival tissue of patients who developed ECA following trichiasis surgery. Inflammation is a crucial stage of the wound healing process, necessary to attract leukocytes to the wound site to remove debris and pathogens and to trigger the subsequent proliferative stage. Excessive inflammation can lead to fibrotic pathology; however, an increase in pro-inflammatory cytokines was associated with poor prognosis following glaucoma filtration surgery 34 . The factors underlying the increased inflammation in these individuals is unclear; there was no difference in baseline clinical inflammation grade between those that did and did not develop ECA. Risk factors associated with ECA development include baseline trichiasis severity and older age and immediate post-operative appearance 11,35,36 , however there was no difference in baseline age or trichiasis severity in those with and without ECA in this study. Immediate post-operative appearance and its association with ECA development and transcriptional patterns were not assessed in this study and this potentially requires further study. Inter-surgeon variability and distance between sutures were also associated with ECA development following PLTR surgery 36 , however there was no difference between surgeon and ECA outcome in this trial.
A number of pro-fibrotic genes were upregulated in individuals with ECA, including VCAN, FN1, SPARC and FERMT2 (proteins that constitute the ECM or regulate the interactions between cells and the ECM), matrix metalloproteinases MMP1 and MMP9 and pro-fibrotic growth factors CTGF, TGFB1 and TGFB1|1. MMPs are involved in multiple stages of wound healing and are crucial for matrix remodelling during the tissue regeneration process 37 . However, excessive production of MMPs is often associated with pathology and overexpression of MMPs 1, 2, 8, 12 and 13 correlated with levels of corneal fibrosis in patients requiring penetrating and deep anterior lamellar keratoplasty 38 . TGFβ is a key regulator of the wound healing process and regulates cellular differentiation and proliferation, ECM production and wound contraction. TGFβ is also thought to orchestrate fibrotic outcomes following glaucoma filtration surgery 39 and it induces tissue contraction and fibronectin production by conjunctival fibroblasts 40 . Furthermore, increased expression of fibronectin and integrin α5β1 are thought to mediate TGFβ-induced myofibroblast differentiation 41 , which have high contractile activity. The upregulation of these pro-fibrotic and matrix factors in patients with ECA likely reflects excessive remodelling and contraction of the conjunctival tissue.
Type I, III, V, VI and VII collagen transcripts were strongly upregulated in individuals with ECA and were all found in co-expression cluster one. Collagens are fibrous proteins that make up the connective tissue and collagen type I is the major constituent of the ECM 42 . Early on during the wound healing process myofibroblasts are thought to produce collagen type III, which is later replaced by collagen type I 43 . Collagen types III, V and VI have been found to be increased in pathological scarring conditions 44-47 . In tissue from normal conjunctiva examined by immunohistochemistry, collagen types I and III were limited to the substantia propria and collagen type V was absent, whereas in inflammatory and scarring trachoma increased collagen type IV was detected in the thickened basement membrane and new collagen type V was deposited in the upper substantia propria 48,49 . A number of collagen transcripts and proteins, in particular collagen types I, VIII and XI, were associated with pathological scarring in the conjunctiva following glaucoma filtration surgery in mice and humans 15 .
Contrary to our hypothesis, we found no strong evidence for modulation of the association between gene expression and ECA development by doxycycline treatment. In support of this finding there was no evidence that doxycycline treatment reduced the risk of developing post-operative trichiasis or ECA in the overall trial 7 . These findings are somewhat surprising given the evidence supporting the clinical use of oral doxycycline to reduce inflammation and fibrosis in various diseases including rosacea, idiopathic pulmonary fibrosis, chronic blepharitis, corneal erosion and corneal melting post-Pseudomonas infection 21-23,50-52 . Furthermore, doxycycline treatment in vitro reduced MMP1, 9 and 12 expression by primary fibroblasts from trichiasis patients and reduced collagen matrix contraction 24 . Treatment was taken every day for 28 days prior to collection of the one-month sample therefore levels of active drug in the body should have remained high, and the positive impact of oral doxycycline on corneal and other skin diseases suggests that it successfully penetrates the epithelium. Oral doxycycline administered for three weeks was found to have no benefit in patients with persistent symptoms following treatment for neuroborreliosis (a neurological manifestation of tick-borne Lyme disease); there was no difference in symptoms relative to a placebo group and no change in systemic cytokine responses 53 . It is possible that in the context of trichiasis surgery, doxycycline either has no immunomodulatory effect or it does not immunomodulate the factors involved in post-operative trichiasis and ECA.
Although a large number of genes were upregulated in all treatment/outcome groups at one-month post-surgery by RNA sequencing, no differences were detected between groups.
This was probably due to the small sample size and high degree of variation in the one-month samples, evident by PCA. Reducing this variation through summed expression of co-expressed gene clusters revealed differences in wound healing pathways that were validated by qPCR. A strength of this study was the original trial design, including standardisation of surgery, patient treatment adherence and follow up. This study was designed and targets were chosen for follow-up to investigate the immuno-fibrogenic factors associated with the development of ECA, therefore it was not appropriate to reanalyse the results in relation to post-operative appearance, post-operative trichiasis (reported elsewhere 12 ) or ECA regression; in addition, these analyses would have had small group sizes.
This study has revealed that the development of post-operative ECA is strongly associated with the upregulation of pro-fibrotic genes, including growth factors (and the canonical fibrogenic growth factor TGFB1), matrix proteins and collagens. A number of MMPs and pro-inflammatory cytokines were also upregulated. Contrary to our hypothesis, there was no strong evidence to suggest that the anti-inflammatory and anti-fibrotic drug doxycycline modulated the relationship between gene expression and ECA. Our results suggest that other drugs with broad anti-inflammatory and/or anti-fibrotic properties, or targeted inflammatory/fibrotic pathway regulators might show promise in preventing ECA development, however the specific pathways of relevance remain unclear.
Although it is reassuring that ECA regressed in roughly a third of patients by 6 months, the short-term aesthetic appearance may still dissuade others from undergoing trichiasis surgery and anti-fibrotic treatments to prevent ECA development remain desirable.

Open Peer Review
Methods: The nurse-surgeons have been described as 'experienced'. Perhaps a quantification of the number of years of experience (e.g. median or range) of the nurses should be included. Also, how many nurse-surgeons were involved in this study?

2.
Methods: As this is within an RCT, the method of randomisation should be described.

3.
Methods: Were there specific exclusion criteria (e.g. patients already on specific antiinflammatory drugs)?

4.
Methods: Of the 98 samples that failed QC, which groups were the samples from? Was the number of failed samples relatively equal in the moderate/severe ECA group compared to samples taken from individuals without/mild ECA?

5.
Methods: Details of the power calculation for both analyses should be included. Were the sample sizes adequate to detect differences between the groups for the RNA seq experiments? Also, the power of the study is relevant as the study reported results showing no evidence of doxycycline treatment affecting gene expression and ECA development.

6.
Results: The first paragraph of the results have been described in the Methods section and can probably be left out.

7.
Results: For table 1, were there significant differences in the baseline characteristics between the groups? 8.
Results: For table 2, it appears that more genes were up-regulated in poor outcome groups compared to control groups. Also, there were a higher number of differentially expressed genes in the 'doxycycline-good' group compared to 'placebo-good' group. Were these statistically significant?

9.
Results: For table 5, were there significant differences in the baseline clinical and demographic characteristics between the groups? 10.
Discussion: The authors could perhaps include a short paragraph on the relevance of their transcriptome findings on understanding the pathophysiology of ECA development and how this is important in the future development of targeted anti-inflammatory or antifibrotic therapies.

Are sufficient details of methods and analysis provided to allow replication by others? Yes
If applicable, is the statistical analysis and its interpretation appropriate? number of failed samples relatively equal in the moderate/severe ECA group compared to samples taken from individuals without/mild ECA?
As described above, the samples that failed did not differ from those that passed QC by age, sex, or group. A sentence has been added to clarify this and the numbers in each group are shown to reviewer 1 Methods: Details of the power calculation for both analyses should be included. Were the sample sizes adequate to detect differences between the groups for the RNA seq experiments? Also, the power of the study is relevant as the study reported results showing no evidence of doxycycline treatment affecting gene expression and ECA development.
RNAseq experiments were conceived to be hypothesis generating for TLDA target selection. The number of samples selected was limited by RNA quality and sufficiency for RNAseq and the total budget available. The a priori design of this nested study within a clinical trial, was the assessment of gene expression profiles in all people developing ECA (and a matched control group). As the numbers that would go on to develop ECA were not known beforehand a separate sample size calculation for this outcome was not performed and a post-hoc power calculation is not a recommended practice." The first paragraph of the results have been described in the Methods section and can probably be left out.

This has been removed
For table 1, were there significant differences in the baseline characteristics between the groups? Samples were selected for RNA-seq based upon trachomatous scarring and eyelid distortion characteristics at baseline, one and six months, as described on line 183. There was no difference in sex between groups (Chi2 P value 0.763), however age distribution differed slightly between groups (ANOVA P value 0.0143).
For table 2, it appears that more genes were up-regulated in poor outcome groups compared to control groups. Also, there were a higher number of differentially expressed genes in the 'doxycycline-good' group compared to 'placebo-good' group. Were these statistically significant This observation is described in the first paragraph of the results section; As this was a descriptive testing for differences between total numbers of differentially expressed genes was not performed.
For table 5, were there significant differences in the baseline clinical and demographic characteristics between the groups? Groups with no or mild ECA at one month were matched to those with moderate/severe ECA by age, sex and trial arm. There was no difference between groups in baseline conjunctival scarring (Chi2 P value 0.192), age (ANOVA P value 0.0737) or sex (Chi2 P value 0.162). C. trachomatis was detected in only 8 samples therefore differences between groups were not tested. and possibly the regression of ECAs seen in the early post-operative period, but not at later visits. Methods: Exam section, last paragraph: Please clarify that it is the conjunctival scarring that was graded with the FPC system.
Last paragraph before data analysis: Were the samples that failed randomly distributed or were they more commonly from one or more groups -either one of the 4 groups, or did they vary by age and/or sex.
There is a bit of a mix of methods in the results and results in the methods section. Recognizing that sometimes it is important to describe results of intermediate steps in the methods when they influence the next set of methods, please review and consider not repeating specific aspects. For example, the first paragraph of the results is nearly a copy/paste from the methods.
Discussion: Bottom of page 12: Please include discussion regarding prior research that has shown a strong association between immediate post-op appearance and early outcomes. [Merbs, PLOS NTD 2012 reports that 85% of the time ECAs were predicted based on immediate post-op appearance] Minor suggestions: In the introduction the word "recurrence" is used instead of post-operative trichiasis. In keeping with WHO recommendation, it would be good to use the latter instead.
The introduction reports 12-22% "recurrence" within 12 months after surgery and then cites the study used for this analysis and one other by this group. There is a much broader range of PTT rates, and studies have been reported by other groups. This statement should be modified to reflect the range of rates reported in the literature.
Please check references in the 3 rd paragraph of the introduction. At present, the second sentence references the grading scheme for ECAs, but is discussing the corneal damage and vision loss associated with ECAs. Same paragraph, "patient's" should be "patients'" Second to last paragraph of introduction: please add numbers or %s regarding ECAs, as was done with the TT numbers.

Are sufficient details of methods and analysis provided to allow replication by others? Yes
If applicable, is the statistical analysis and its interpretation appropriate? Included samples 300 (38.8%) 92 (11.9%) 295 (38.1%) 87 (11.2%) We removed the repeated paragraph from the results.
We added the suggested Merbs et al citation to the discussion All instances of 'recurrence' in the paper have been changed to 'post-operative trichiasis'.
We agree there is a wider range of reported outcomes. This has been adjusted to reflect this. We believe the most suitable references for this is the review conducted by Rajak et al (Trachomatous Trichiasis and its management in endemic countries) and Cochrane review by Burton et al (Interventions for trachomatous trichiasis). We have revised the sentence and added these two references.

We have removed the incorrect citation in the introduction and corrected patients
We have added figures regarding ECAs as requested.
Competing Interests: No competing interests were disclosed.