Figures
Abstract
Purpose
Trachoma is a fibrotic disease of the conjunctiva initiated by Chlamydia trachomatis infection. This blinding disease affects over 40 million people worldwide yet the mechanisms underlying its pathogenesis remain poorly understood. We have investigated host microRNA (miR) expression in health (N) and disease (conjunctival scarring with (TSI) and without (TS) inflammation) to determine if these epigenetic differences are associated with pathology.
Methods
We collected two independent samples of human conjunctival swab specimens from individuals living in The Gambia (n = 63 & 194). miR was extracted, and we investigated the expression of 754 miR in the first sample of 63 specimens (23 N, 17 TS, 23 TSI) using Taqman qPCR array human miRNA genecards. Network and pathway analysis was performed on this dataset. Seven miR that were significantly differentially expressed between different phenotypic groups were then selected for validation by qPCR in the second sample of 194 specimens (93 N, 74 TS, 22 TSI).
Results
Array screening revealed differential expression of 82 miR between N, TS and TSI phenotypes (fold change >3, p<0.05). Predicted mRNA targets of these miR were enriched in pathways involved in fibrosis and epithelial cell differentiation. Two miR were confirmed as being differentially expressed upon validation by qPCR. miR-147b is significantly up-regulated in TSI versus N (fold change = 2.3, p = 0.03) and miR-1285 is up-regulated in TSI versus TS (fold change = 4.6, p = 0.005), which was consistent with the results of the qPCR array.
Author Summary
Trachoma is a debilitating disease that affects 40 million people worldwide. It can cause progressive fibrosis of the upper eyelid and blindness, yet the mechanism is poorly understood. We have investigated the expression of short sequences of genetic material (microRNA) that regulate gene expression. We screened for the expression of 754 microRNA sequences (miR) in genetic material isolated from conjunctival swab samples from individuals in trachoma-endemic communities in The Gambia. This sample included healthy controls, individuals with trachomatous scarring and individuals with trachomatous scarring in the presence of clinically significant inflammation. We found 82 miR that were differentially expressed. Computer simulations predict that these miR regulate genes in epithelial cell differentiation, inflammation and fibrosis pathways, all of which are involved in the scarring process. We then validated the expression of seven of these differentially expressed miR in a second larger biological sample set from The Gambia. We confirmed that miR-147b and miR-1285 have increased expression in individuals with trachomatous scarring in the presence of clinically significant inflammation. Further investigation into the functions of these miR will aid our understanding of this disease and present opportunities to develop treatments for ocular fibrotic diseases.
Citation: Derrick T, Roberts Ch, Rajasekhar M, Burr SE, Joof H, Makalo P, et al. (2013) Conjunctival MicroRNA Expression in Inflammatory Trachomatous Scarring. PLoS Negl Trop Dis 7(3): e2117. https://doi.org/10.1371/journal.pntd.0002117
Editor: Joseph M. Vinetz, University of California San Diego School of Medicine, United States of America
Received: December 11, 2012; Accepted: January 31, 2013; Published: March 14, 2013
Copyright: © 2013 Derrick et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This study was funded by grants from the Wellcome Trust (079246/Z/06/Z and GR079246MA) with additional support from Fight for Sight for TD in the form of a graduate studentship (http://www.wellcome.ac.uk; http://www.fightforsight.org.uk). The funders had no part in the study design; in the collection, analysis, and interpretation of data; in the writing of the report; and in the decision to submit the paper for publication.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Chlamydia trachomatis (Ct) is the causative agent of trachoma, the leading cause of blindness that results from infection. Forty million people have active trachoma and eight million people suffer with unoperated trichiasis [1]. Repeated infection of the conjunctiva by this intracellular bacterium during childhood causes a chronic inflammatory response, leading to progressive fibrosis and scarring in adult life. Scarring distorts the conjunctiva and the eyelashes are pulled inward to the extent that they scratch the cornea (trichiasis), causing pain and eventually blindness.
Chronic trachomatous inflammation is known to continue in the absence of current Ct infection and is believed to drive the continued scarring process, however the mechanisms by which this occurs are not completely understood [2]. Messenger RNA (mRNA) expression profiling of each clinical stage of trachoma has revealed many thousands of mRNAs that are differentially expressed [2], [3]. Key pathways that are differentially regulated in the conjunctiva are innate inflammatory pathways and extracellular matrix modifiers. MicroRNAs (miR) are known to have significant roles in the regulation of inflammation, fibrosis and cell differentiation [4]–[7] and can be dysregulated upon bacterial infection [8], [9].
miR are post-transcriptional regulators of gene expression. They are single-stranded RNA molecules typically 18–22 nucleotides in length. miR bind to complementary mRNA sequences in association with the RNA-induced silencing complex (RISC), causing transcriptional degradation of the transcript or repression of its translation [10]. The seed sequence (5′ nucleotides 2–7) of the miR guides target selection [11]. Complementary target sequence sites are usually, though not exclusively, found in the 3′ untranslated region (UTR) of mRNA transcripts. For a given miR complementary sequence sites may be present on a few or several hundred different mRNA targets, indicating the potential for a few miRs to regulate complete biological processes. Indeed, a relatively small total number of miR are thought to regulate over a third of all protein-coding genes [12]. miR have profound roles in the regulation of many biological processes and interest in their various functions in health and disease is growing. The number of known mature miR (http://www.mirbase.org/) is increasing rapidly as research in this area quickly unravels miR biology.
The ability of miR to regulate entire pathways offers investigators an opportunity to reduce the complexity of the trachoma transcriptome [2]. We suggest that the differential regulation of just a few miR in trachomatous disease may underlie the substantial differences in mRNA expression that characterize each phenotypic trachoma group. This reduction in complexity will enable more targeted research into the mechanisms of disease and may identify new potential therapeutic approaches.
Methods
Ethical permission, study participants and clinical assessment
The study was conducted in accordance with the tenets of the Declaration of Helsinki. The Ethics Committee of the Gambian Government/Medical Research Council Unit, and the ethics committee of the London School of Hygiene and Tropical Medicine approved the study. Samples were drawn from an archive built up under the MRC study numbers SCC729, SCC1177 and SCC1274 with specific approval for miR gene expression studies under SCC L2011.03. The samples described in this paper were collected from individuals recruited in trachoma-endemic communities across the whole geographic range of The Gambia, West Africa. Written informed consent was taken from individuals at the time of sample collection. For those participants aged <16 years that wished to take part in the study consent was obtained from a parent/guardian. All samples were anonymized. Cases of trachomatous scarring (TS) were identified from screening records, community ophthalmic nurse referral and opportunistic rapid screening. Control individuals with normal conjunctivae were selected by matching for age, sex, ethnicity and location. Clinical phenotypes were assessed in the field by experienced field supervisors trained and regularly assessed in trachoma grading. FPC scores [13] were assigned and grades were agreed by two experienced trachoma physicians using high-resolution photographic records taken in the field using a Nikon D3000 SLR camera with a VR AF-S micro Nikkor 105 mm 1:2.8G ED lens. Photographs were taken at the time of sample collection. Individuals were grouped into the following clinical phenotypes for analysis: individuals with trachomatous scarring (TS) had a C score between 1–3 (mild to severe scarring) and a P score of 0 or 1 (none or mild inflammation), individuals with trachomatous scarring in the presence of clinically significant inflammation (TSI) had a C score of 1–3 and a P score of 2 or 3 (moderate to severe inflammation), and control samples from individuals with normal healthy conjunctivae (N) had no conjunctival scarring (C0), papillary inflammation (P0) or follicles (F0).
Sample collection and processing
Swabs were taken from the upper tarsal conjunctiva using Dacron polyester-tipped swabs (Hardwood Products Company) and stored in 250 µl RNAlater (Ambion, Life Technologies) on ice blocks in the field and then archived at −20°C until processed.
Taqman microfluidic array miR genecards
A total of 63 specimens from the archive were selected for miR expression array profiling. Specimens were selected as representative examples of each phenotype group using the FPC scores. As control samples for these experiments, individuals with normal conjunctivae matched on age, sex, ethnicity and location were selected. MiR was extracted from swabs using the Qiagen Allprep DNA/RNA/protein kits with a modification to collect small non-coding RNAs. DNase1 digestion (Qiagen) was included. Total RNA purity was assessed by spectrophotometry using a nanodrop ND-1000 (Thermo Fisher Scientific). Reverse transcription and pre-amplification were performed using Megaplex human primer pools Av2.1 and Bv3.0 following the manufacturer's instructions (Taqman, Life Technologies). Quantitative PCR was performed using 72 µl of pre-amplified cDNA as template in the PCR master mix for the TaqMan Array Human MicroRNA genecards (Av2.0 and Bv3.0). Thermal cycling was performed on a 7900HT thermal cycler (Life Technologies). Plates were held at 50°C for 2 minutes, 94.5°C for 10 minutes, then underwent 40 cycles of 97°C for 30 seconds and one minute at 59.7°C. A total of 754 of the most well characterised unique human miR from Sanger miRBase V.14 (www.mirbase.org/) were screened. Sanger miRBase V.14 was the latest version of the miR database at the time of screening.
Array analysis
qPCR cycle threshold (CT) values were processed in SDS RQ manager (Life technologies); the threshold was set to 0.05 and baselines were detected automatically. Data from each array were uploaded and analysed using the High Throughput qPCR Package (HTqPCR) in Bioconductor R (www.bioconductor.org, www.r-project.org) [14]. Sample profiles were excluded from the analysis when the median CT value for the array was 40 since the majority of the CT values were either close to threshold or undetermined. Individual miR were retained in the analysis only when expressed (CT value<40) by at least five specimens. A and B genecard data were analysed separately due to differences in specimen performance on each card. Data were normalized to reduce technical bias in the analysis by a number of different standard methods (supplementary figure S1A & B). The coefficient of variation, standard deviation and correlation of raw against normalized data were used to evaluate the suitability of each method of normalisation, as described in Deo et al. (2011) [15]. The ‘Norm rank invariant’ method was chosen as the most effective normalisation strategy for both A and B cards (supplementary figure S1A & B). The distribution of the raw and Norm rank invariant normalized CT values are shown in supplementary figure S2 A–D. Differential expression was then assessed by empirical Bayes/moderated t-tests using HTqPCR in Bioconductor R. These data are deposited within the NCBI GEO public database (GSE37717) and, in line with MIQE guidelines [16], details are included as supplementary table S1.
MiR abundance
Relative abundance of miR in the conjunctiva was calculated from array CT data. An average was taken of the normalized CT values for all specimens (including all phenotypes) for each miR on A and B cards. The general equation for estimating relative differences in PCR was then applied to these values: 2−(CT_target−CT_calibrator) where the calibrator was the most abundant miR (miR-1274B). For each miR this value was then divided by the sum of all these values to create a relative abundance.
Network analysis
A network graph based on the specimen-to-specimen Pearson correlation was generated Biolayout express 3D v2.2 (www.biolayout.org/) [17]. The overall miR expression correlation matrix and graph were constructed from the rank invariant normalized raw CT data. The graph was arranged according to patient clinical classification (N, TS, TSI). Pearson correlation coefficients (r)>0.7 were retained and used as cut-offs in network construction. Nodes in the graph are individual miR linked by an edge if r>0.7. The graph was then clustered using a Markov Clustering algorithm using default inflation values. The partitioned clusters of expression contain sets of miR that exhibit a very strong degree of co-expression across the sample. Cluster content is independent of differential expression level. The co-expression clusters were then investigated in silico at the individual and pathway levels.
Pathway analysis
miR that were differentially expressed in the arrays at a significance level of p<0.05 and with a fold change (FC) over three (up or down-regulated) were entered into pathway analysis. DIANA-microT v4.0 (Beta version) target prediction was used in the DIANA mirPath software [18]. Multiple miRNA analysis was used for the significant miR within each comparison group.
Extraction of miR for quantitative real-time PCR
Total RNA including miR was extracted from swabs using a Qiagen miRNeasy kit, incorporating a DNase1 digestion step. Total RNA purity was assessed by spectrophotometry using a nanodrop ND-1000 (Thermo Fisher Scientific). miR was reverse transcribed using miScript II RT kit with the hiFlex buffer as per the manufacturer's instructions (Qiagen). qPCR was carried out using miScript Primer Assays and the miScript SYBR Green PCR kit (Qiagen) on a 7900HT thermal cycler (Life Technologies). Ten microlitres of RT product was diluted in 100 µl H20 and 0.5 µl was used as template in each qPCR assay, with 0.5 µl specific miR forward primer assay, 1 µl water, 0.5 µl universal reverse primer and 2.5 µl SYBR green master mix, in a total reaction volume of 5 µl. Each assay was performed in quadruplicate, including no template controls for each miR and for each specimen. Cycling conditions were as follows: 15 minutes at 95°C, followed by 70 cycles of 15 seconds at 94°C, 30 seconds at 55°C, and 30 seconds at 70°C. Data were collected at 94°C and 70°C. qPCR was run for 70 cycles to minimize the number of undetermined values. Fifteen percent of all miR assays had a CT value over 40.
qPCR analysis
CT values were derived in SDS RQ manager (ABI, Life technologies), with a threshold of 0.05 and an automatic baseline. Four replicate tests were used to calculate the geometric mean after outliers were removed. Analysis was done in R. Specimens were removed from the analysis if the endogenous control snoU6 (U6) CT values were ≥2× the standard deviation (s.d.) of the mean of all U6 CTs in the sample. For analysis purposes, any assay that did not amplify by 70 cycles was assigned a CT value of 80 (19 assays out of a total 1552). Target CT values for each specimen were normalized to U6 using Δ CT = CT target−CT U6. For each miR within each phenotype group, the Shapiro-Wilk method was used to test for normality of distribution in the raw-data. For each comparison (N v TS, N v TSI, TS v TSI) we calculated the fold change in miR expression between the phenotypic groups, using 2−ΔCT(median phenotype 1) −ΔCT(median phenotype 2). We used the Wilcoxon rank sum test (with continuity correction) to test for differences in the expression of each miR between phenotypic groups as the majority of the data were not normally distributed. Details of the qPCR and analysis in line with MIQE are included in supplementary table S1.
Results
Array samples
Sixty-three specimens were tested by miRNA array cards. Twenty-three specimens were excluded from A genecard profiles and 34 specimens from B genecard profiles. Basic demographic and clinical phenotype data are shown in table 1 both before and after filtering and these show that there was no systematic loss of any specific sample type based on clinical phenotype, age or sex as a result of the filtering process.
Basic characteristics of miR expression in the conjunctiva: Abundance and co-expression analysis
Following the filtering procedures described, 506/754 miRs were included in the final analysis. Relative abundance of all of 506 expressed miR in the conjunctiva was calculated from the CT data. Of the miR that were tested, just six constitute 90% of total miR present in the conjunctival samples (figure 1). miR-1274B has the highest overall expression and miR-623 the lowest. Networks of co-expression, independent of differential expression, based on rank invariant CT values were explored in the entire data set of 506 miR using Biolayout express 3D. The undirected graph contained 126 miR connected by 215 edges. Markov clustering partitioned the network into 11 clusters of co-expressed miR. These clusters ranged in size from 23 to 4 co-expressed miR and accounted for 80 miR in the original network. Each of the 11 clusters is laid out in figure 2 and the specific miR content of each cluster is available in supplementary table S2. The major biology revealed by these co-expression clusters indicates that these miR target mRNA in four major pathways. These are repeatedly identified and are shown in figure 2. The MAPK signaling pathway and focal adhesion pathway contain the largest number of miR target genes and have the highest levels of enrichment (over-representation).
Abundance of all miR tested expressed relative to miR-1274B. Abundance was calculated from cycle threshold values irrespective of sample phenotype.
Clusters of co-expressed miR are shown, calculated from normalized array data irrespective of phenotype. Clusters are highlighted according to the pathway most enriched upon pathway analysis of miR in each cluster.
Differential expression analysis
A total of 82 miR are differentially expressed across the comparison groups (FC>3, p<0.05) (supplementary table S3). The number of up- and down-regulated miR in each comparison is shown in table 2. A miR that is up-regulated in the N v TS comparison has a lower CT value in TS relative to N. The same applies to N v TSI and TS v TSI, where the latter is up- or down-regulated relative to the former in each comparison. A larger number of miR are differentially expressed in comparisons with the TSI phenotype. Fewer miR are down-regulated than up-regulated, particularly in the N v TSI group.
Twenty miR are differentially expressed in both N v TSI and TS v TSI comparisons (figure 3) indicating they might be features of inflammation. In contrast, very few miR are shared with the N v TS miR gene set and there are none that overlap between all three groups. This indicates that TSI and TS phenotypes are distinct and have characteristic miR signatures.
Venn diagram showing the number of differentially expressed (FC>3, p<0.05) miR that are unique or that overlap between the different clinical phenotypes. Selected miR of interest are shown with arrows illustrating whether they are up- or down-regulated in the indicated comparison group. An upward facing arrow indicates up-regulation and a downward facing arrow indicates down-regulation. * miR-23b is up-regulated in N v TS but down-regulated in the TS v TSI comparison group.
Of the 103 miR found in the networks by Markov clustering, 15 had some evidence of differential regulation based on p-value alone. miR-492 and miR-548d were both up-regulated (4.9 and 3.2 times respectively) in TSI individuals compared to controls whilst 3 miR (miR-508, miR-509 and miR-664) were >3 fold down-regulated. The miR precursor let-7b showed modest evidence of differential up-regulation in TS individuals versus normal controls (p = 0.037).
Pathway analysis
Differentially expressed miR in each comparison group (FC>3 p<0.05) were entered into DIANA mirPath pathway analysis. The top ten most enriched pathways for each comparison are listed in table 3, where a greater −ln(p-value) reflects increasing enrichment of miR targets within a pathway. Many of the same pathways are enriched in each comparison group, despite little overlap in miR between the N v TS, N v TSI and TS v TSI groups (figure 3). Axon guidance, focal adhesion, and the TGF-β signaling pathway are all present in all three groups. A large number of genes in the TGF-β pathway are predicted targets of differentially expressed miR in the N v TSI comparison (figure 4). Within the TGF-β pathway, 53% of transcripts are differentially regulated based on differences found in a mRNA transcriptome array using Ethiopian conjunctival samples (GSE23705) from similar phenotypic groups [2]. Analysis was also carried out on each phenotypic comparison group split into up- and down-regulated gene sets (supplementary table S5). Interestingly, TGF-β is enriched in the down-regulated gene set for each comparison group. Given that this pathway is enriched for miR targets that would be silenced and these miR are down-regulated, this would suggest an up-regulation or release of the TGF-β signaling pathway. Analysis was also performed on the 20 miR differentially expressed in both N v TSI and TS v TSI comparisons, with enrichment again in MAPK, TGF-β and Wnt pathways (supplementary table S6), supporting the hypothesis that these miR are characteristic of the major pathways under miR control in the conjunctiva and are perturbed by inflammation.
Genes highlighted in yellow are predicted targets of miR differentially regulated in N v TSI (35 miR) (FC>3 p<0.05). Genes highlighted in green are not predicted targets.
Validation of differential expression of selected miR using an alternative RT-qPCR system
In the validation set it was not feasible to assay the expression of all the potentially differentially expressed miR and we selected for follow-up a small number of miR that exhibited a high FC, low p-value and homogeneity in the raw data. Seven miR were selected for follow-up (supplementary table S7) including three miR that were differentially regulated in the N v TS comparison (miR-30c, miR-32, miR-203) and four from the N v TSI comparison (miR-10a, miR-147b, miR-1285, miR-1305). Each candidate miR was tested for differential expression in a second sample of 194 independent archival Gambian clinical specimens, selected as representative examples of each phenotype group. In these experiments, small nucleolar (sno) U6 RNA was used as the calibrator snoU6 CT values were not different between the phenotypic groups, which implied that it was a stably expressed reference gene (supplementary information figure S3). Five specimens were excluded because they had outlying snoU6 values (average CT>2 s.d. of the sample mean), leaving a total of 189 specimens to be tested for statistical differences in expression levels. Summary statistics for these 189 specimens are shown in table 4.
Data were tested for differential expression between the three phenotypic groups as is presented for the analysis of the array data (N, TS and TSI). Only miR-1285 and miR-147b showed a significant difference between the different phenotypic groups in this validation set (table 5). MiR-147b was up-regulated 2.3 fold in individuals with TSI relative to N (p = 0.0332). This is consistent with the array results in which miR-147b was up-regulated 9.6 fold in individuals with TSI versus N. MiR-1285 was up-regulated 4.6 fold in TSI relative to TS (p = 0.005). This is also consistent with the array results in which miR-1285 was up-regulated 16 fold in TSI versus TS.
Discussion
Array analysis revealed that a large number of miR are potentially differentially regulated between different disease states and healthy controls. Trachomatous scarring with inflammation (TSI) has a distinct miR signature compared to scarring trachoma (TS). TS may be a less active disease process than TSI as fewer miR were differentially regulated. On validation, two miR remained significantly differentially regulated. MiR-147b was up-regulated in individuals with TSI compared to N and miR-1285 was up-regulated in people with TSI compared to those with TS alone. In a transcriptome analysis of similar phenotypic comparison groups in Ethiopians with scarring trachoma [2], 25% of predicted targets of miR-1285 and 52% of predicted targets of miR-147b predicted targets (TargetScan v6.2) were differentially regulated (adjusted p<0.05).
MiR-1285 directly targets the 3′UTR of p53 mRNA in HEK 293T cells [19]. The loss of p53 is associated with many cancers via disruption of the normal function of p53 in the initiation of apoptosis and growth arrest. In contrast, Hidaka and colleagues [20] find miR-1285 to be a tumor suppressor. Expression levels of miR-1285 were reduced in clinical samples of renal cell carcinoma (RCC) compared to healthy mucosa and miR-1285 transfection in RCC cell lines in vitro led to inhibition of cell proliferation, migration and invasion [20]. The authors verified transglutaminase 2 (TGM2) as a target of miR-1285. Interestingly, TGM2 is linked to several cancers and the process of epithelial-mesenchymal transition (EMT) [21]. EMT can be initiated by chronic inflammation [22] and is implicated in the pathology of many fibrotic diseases [23], so could have a role in trachomatous disease. The conflicting conclusions of these studies are likely due to the different cell types used. In vitro studies in primary conjunctival epithelia will be required to understand the function of miR-1285 in chlamydial infection and trachoma.
There is limited literature on the clinical significance of changes in miR-147b expression, but it is known to be down-regulated in rectal cancer [24]. This is consistent with the ability of miR-147b to induce apoptosis four days post-transfection in A549 cells [25]. Mmu-miR-147, a functional homologue of human miR-147b, is induced by multiple toll-like-receptor (TLR) signals and negatively regulates inflammation in murine macrophages [26]. It acts in a negative feedback loop to prevent excessive inflammation. Bertero and colleagues [25] showed that LPS (bacterial lipopolysaccharide) and TNFα strongly induced hsa-miR-147b expression in A549 cells in vitro; supporting the hypothesis that miR-147b has a homologous role in the regulation of inflammation in humans. The up-regulation of miR-147b seen in the TSI phenotypic group may reflect ongoing and uncontrolled inflammation. Understanding the role of miR-147b in TSI should also be aided by further in vitro functional studies.
MiR-23b-5p is up-regulated in TS relative to both N (FC = 3.6 p = 0.026) and TSI (FC = 5.8 p = 0.008) (figure 3), indicating that it might be a feature of scarring in the absence of inflammation. MiR-23b is a member of the miR-23b cluster, which includes miR-27b and miR-24-1. This cluster is known to target members of the TGF-β signaling pathway [27]. An up-regulation of miR-23b, as is seen in scarred individuals, would lead to a decrease in TGF-β signaling. Although this may seem counterintuitive, any dysregulation of the TGF-β signaling pathway could lead to aberrant wound healing [28]. MiR-23b is also anti-inflammatory through inhibition of NFkB activation [29], preventing up-regulation of inflammatory cytokines such as IL-17. In turn, IL-17 inhibits miR-23b, leading to inflammation. An up-regulation of IL-17A is associated with active trachomatous disease [30]. The relative down-regulation of miR-23b in TSI compared to TS and N conjunctival samples could reflect a down-regulation of this miR's expression by IL-17 in individuals with TSI.
MiR-30c was up-regulated 15 fold in N v TS (p = 0.01) and 11 fold in N v TSI (p = 0.04) in the microarray experiments. This miR is thought to regulate fibrinolysis and collagen production through targeting serine protease inhibitor SERPINE1 and connective tissue growth factor (CTGF) respectively [31]. Over expression of this miR is known to inhibit the proliferative and migratory properties of endometrial cancer cells [32], however qPCR validation found no association of this miR with disease.
Pathway analysis revealed that many of the same pathways are enriched amongst predicted targets of differentially expressed miR in each comparison group. This is surprising due to the minimal overlap of differentially expressed miR in N v TS compared with the other groups. Many of these pathways are characteristic of epithelial cell and fibroblast communication, differentiation and fibrosis. In particular, the Wnt pathway has been implicated in the disruption of epithelial cell homeostasis and the pathology of C. trachomatis infection [33]. Importantly, epithelial cell differentiation, development and cytoskeleton remodeling pathways (including TGF-β and Wnt) are enriched in gene sets from a differential expression analysis of transcriptome data of similar phenotypic groups in Ethiopians [2]. Of all the members of the KEGG defined TGF-β pathway, 53% are differentially regulated in this Ethiopian transcriptome [2]. We theorise that the miR that are differentially expressed in this study are at least partly responsible for the observed changes in the normal function of the Wnt and TGF-β pathways in trachoma patients.
Investigation of miR abundance in the conjunctiva shows that most miR have very low levels of expression, with just a few making up the vast majority of the population. MiR-1274B is the most highly abundant miR in the conjunctiva, but the nature of this miR has been called in to question. A previous study [34] found evidence that miR-1274B may not be a canonical miR, rather that it is a tRNA-derived small RNA (tsRNA). tsRNAs are thought to be abundant in the genome [26] and may have a role similar to miR in regulating gene functions [35]. It may be the case that the high abundance of miR-1274B can be explained simply by its origin in tRNA and the generally high abundance of tRNA in cells. Regardless of this uncertainty over its designation as either a miR or a tsRNA, a role for miR-1274B in the regulation of gene expression cannot be excluded at this time.
The discrepancy between the array results and qPCR validation could be due to a high number of false positives accepted in the initial analysis or as a result of a number of technical differences leading to differential miR isolation, extraction, amplification bias and normalisation [36]. Different methods of miRNA isolation and qPCR were employed in this study in the screening and validation stages which could introduce technical variation [37]. In addition, a non-proscriptive filtering process was used in the identification of potentially differentially expressed miR. Acceptance of a high number of potential false positive associations at this initial filtering stage was considered acceptable based on the arguments presented by Rothman [38] on principles of p-value adjustment and enabled us to explore the wider biology of miR conjunctival expression (supplementary table S4). As a result, many of the miR chosen for follow up in the validation clinical samples appear not to be differentially regulated. This highlights the need to verify array-profiling results even when the method of choice is the apparently robust gold standard technique to examine differential expression [39]–[41]. Furthermore the design of the B genecards, which cover less abundant miR, adds difficultly to analysis pipeline. Even with pre-amplification the amount of genetic material that can be obtained from a conjunctival swab is very small. Specimens were judged to perform less well on B genecards leading to the exclusion of a larger number of specimens, resulting in a reduction of the number of biological replicates in each phenotypic group and therefore some loss of statistical power. The raw and analysed data is publicly available in the NCBI GEO.
This is the first description and initial identification of specific miR expression in an ocular fibrotic disease of humans. It is possible that these miR play a role in other ocular surface inflammatory diseases. We highlight the major pathways under the control of miR that are expressed in the conjunctival epithelium and suggest that it is plausible that dysregulation of expression in these miR leads to the release of the TGF-β signaling pathway in trachomatous disease. Two miR with significantly increased expression in trachomatous scarring and inflammation were identified (miR-147b and miR-1285). In order to understand the mechanisms by which these miR may contribute to health and disease of the conjunctiva the associations of miR-147b and miR-1285 with trachomatous disease now requires further study. A combination of in vitro experimentation in model systems and in vivo application in animal models will also facilitate our understanding of this association and whether these miR are reflective or causative effectors of disease. Research in this area of RNA biology is a rapidly evolving field that is only now beginning to realize its potential. We hope that its application to trachomatous disease may lead to the development of therapeutics or biomarkers for the diagnosis and treatment of trachoma and other fibrotic ocular pathologies.
Supporting Information
Figure S1.
Supplementary figure S1A & B. Coefficient of variation (CV) for uncorrected and normalized array data for A and B card genesets. Various methods of normalization were tested (Left to right: uncorrected data, geometric mean, quantile, delta-CT using RQ manager, delta-CT using geNorm, norm rank invariance, scale rank invariance).
https://doi.org/10.1371/journal.pntd.0002117.s001
(TIFF)
Figure S2.
Supplementary figure S2A–D. Boxplots of CT distribution for each sample in A and B geneset groups after quality control filtering before normalization (A, B), and after normalization (C, D).
https://doi.org/10.1371/journal.pntd.0002117.s002
(TIFF)
Figure S3.
Boxplots of endogenous control snoRNA U6 raw CT distribution in each phenotype group. These distributions are not significantly different (Kruskal-Wallis p = 0.5469).
https://doi.org/10.1371/journal.pntd.0002117.s003
(TIFF)
Table S1.
MIQE guidelines for array and qPCR processing.
https://doi.org/10.1371/journal.pntd.0002117.s004
(XLS)
Table S2.
MiR in the 11 clusters derived from Markov chain clustering algorithm of array data. Footnote: The most enriched pathway among predicted targets of miR in each cluster is shown with values and FC for each miR from the array analysis.
https://doi.org/10.1371/journal.pntd.0002117.s005
(XLSX)
Table S3.
Differentially regulated miR (FC>3 p<0.05) in each phenotypic comparison group in the array analysis.
https://doi.org/10.1371/journal.pntd.0002117.s006
(XLSX)
Table S4.
Number and percentage of false positives in trachoma miR array study. Footnote: With no inflation of p-values the chance of finding one or more significant differences in 1518 tests = 100%. Bonferroni's adjustment indicates critical p-value = 3.38×10−5 and assuming outcomes are moderately correlated (r = 0.5) then a critical p-value = 0.001 would be required.
https://doi.org/10.1371/journal.pntd.0002117.s007
(XLSX)
Table S5.
Pathway analysis of up- and down-regulated miR in the array analysis (FC>3 p<0.05).
https://doi.org/10.1371/journal.pntd.0002117.s008
(XLSX)
Table S6.
Pathway analysis on 20 miR differentially expressed (FC>3 p<0.05) in both N v TSI and TS v TSI comparisons.
https://doi.org/10.1371/journal.pntd.0002117.s009
(XLSX)
Table S7.
MiScript primer assays used in qPCR validation: accession number and mature miR sequence.
https://doi.org/10.1371/journal.pntd.0002117.s010
(XLSX)
Acknowledgments
We are very grateful to the study participants and the support staff at MRC Unit, The Gambia. We would also like to acknowledge the helpfulness and co-operation of the Gambian National Eye Health Programme staff and community ophthalmic nurses in the identification of individuals recruited in this study.
Author Contributions
Obtained funding: DCWM RLB MJH MJB. Conceived and designed the experiments: MJH DCWM. Performed the experiments: TD MJH MR ChR. Analyzed the data: TD ChR MJH RLB MJB. Contributed reagents/materials/analysis tools: SEB PM HJ TD MJH. Wrote the paper: TD MJH.
References
- 1. Mariotti SP, Pascolini D, Rose-Nussbaumer J (2009) Trachoma: global magnitude of a preventable cause of blindness. Br J Ophthalmol 93: 563–568 .
- 2. Burton MJ, Rajak SN, Bauer J, Weiss HA, Tolbert SB, et al. (2011) Conjunctival transcriptome in scarring trachoma. Infect Immun 79: 499–511 .
- 3. Natividad A, Freeman TC, Jeffries D, Burton MJ, Mabey DCW, et al. (2010) Human conjunctival transcriptome analysis reveals the prominence of innate defense in Chlamydia trachomatis infection. Infect Immun 78: 4895–4911 .
- 4. Oertli M, Engler DB, Kohler E, Koch M, Meyer TF, et al. (2011) MicroRNA-155 Is Essential for the T Cell-Mediated Control of Helicobacter pylori Infection and for the Induction of Chronic Gastritis and Colitis. J Immunol 187: 3578–3586 .
- 5. Gregory PA, Bert AG, Paterson EL, Barry SC, Tsykin A, et al. (2008) The miR-200 family and miR-205 regulate epithelial to mesenchymal transition by targeting ZEB1 and SIP1. Nat Cell Biol 10: 593–601 .
- 6. Cushing L, Kuang PP, Qian J, Shao F, Wu J, et al. (2011) miR-29 is a major regulator of genes associated with pulmonary fibrosis. Am J Respir Cell Mol Biol 45: 287–294 .
- 7. Lu TX, Hartner J, Lim E-J, Fabry V, Mingler MK, et al. (2011) MicroRNA-21 limits in vivo immune response-mediated activation of the IL-12/IFN-gamma pathway, Th1 polarization, and the severity of delayed-type hypersensitivity. J Immunol 187: 3362–3373 .
- 8. Zhu Y, Jiang Q, Lou X, Ji X, Wen Z, et al. (2012) MicroRNAs up-regulated by CagA of Helicobacter pylori induce intestinal metaplasia of gastric epithelial cells. PLoS One 7: e35147 .
- 9. Izar B, Mannala GK, Mraheil MA, Chakraborty T, Hain T (2012) microRNA Response to Listeria monocytogenes Infection in Epithelial Cells. International Journal of Molecular Sciences 13: 1173–1185 .
- 10. Martinez J, Patkaniowska A, Urlaub H, Lührmann R, Tuschl T (2002) Single-Stranded Antisense siRNAs Guide Target RNA Cleavage in RNAi. Cell 110: 563–574 .
- 11. Vejnar CE, Zdobnov EM (2012) miRmap: Comprehensive prediction of microRNA target repression strength. Nucleic Acids Res 40: 11673–83 .
- 12. Lewis BP, Burge CB, Bartel DP (2005) Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell 120: 15–20 .
- 13.
Dawson CR, Jones BR TM (1981) Guide to Trachoma Control. Geneva: World Health Organization.
- 14. Dvinge H, Bertone P (2009) HTqPCR: high-throughput analysis and visualization of quantitative real-time PCR data in R. Bioinformatics 25: 3325–3326 .
- 15. Deo A, Carlsson J, Lindlöf A (2011) How to choose a normalization strategy for miRNA quantitative real-time (qPCR) arrays. J Bioinform Comput Biol 9: 795–812.
- 16. Bustin SA, Beaulieu J-F, Huggett J, Jaggi R, Kibenge FSB, et al. (2010) MIQE précis: Practical implementation of minimum standard guidelines for fluorescence-based quantitative real-time PCR experiments. BMC Mol Biol 11: 74 .
- 17. Goldovsky L, Cases I, Enright AJ, Ouzounis CA (2005) BioLayout(Java): versatile network visualisation of structural and functional relationships. Appl Bioinformatics 4: 71–74.
- 18. Papadopoulos GL, Alexiou P, Maragkakis M, Reczko M, Hatzigeorgiou AG (2009) DIANA-mirPath: Integrating human and mouse microRNAs in pathways. Bioinformatics 25 (15) 1991–1993 .
- 19. Tian S, Huang S, Wu S, Guo W, Li J, et al. (2010) MicroRNA-1285 inhibits the expression of p53 by directly targeting its 3′ untranslated region. Biochem Biophys Res Commun 396: 435–439 .
- 20. Hidaka H, Seki N, Yoshino H, Yamasaki T, Yamada Y, et al. (2012) Tumor suppressive microRNA-1285 regulates novel molecular targets: aberrant expression and functional significance in renal cell carcinoma. Oncotarget 3: 44–57.
- 21. Shao M, Cao L, Shen C, Satpathy M, Chelladurai B, et al. (2009) Epithelial-to-mesenchymal transition and ovarian tumor progression induced by tissue transglutaminase. Cancer Res 69: 9192–9201 .
- 22. Kalluri R, Weinberg RA (2009) Review series The basics of epithelial-mesenchymal transition. J Clin Invest 119 (6) 1420–8 .
- 23. Kalluri R, Neilson EG (2003) Epithelial-mesenchymal transition and its implications for fibrosis. J Clin Invest 112: 1776–1784 .
- 24. Gaedcke J, Grade M, Camps J, Søkilde R, Kaczkowski B, et al. (2012) The Rectal Cancer microRNAome - microRNA Expression in Rectal Cancer and Matched Normal Mucosa. Clin Cancer Res 18: 4919–4930 .
- 25. Bertero T, Grosso S, Robbe-Sermesant K, Lebrigand K, Hénaoui I-S, et al. (2012) “Seed-Milarity” Confers to hsa-miR-210 and hsa-miR-147b Similar Functional Activity. PLoS One 7: e44919 .
- 26. Liu G, Friggeri A, Yang Y, Park Y-J, Tsuruta Y, et al. (2009) miR-147, a microRNA that is induced upon Toll-like receptor stimulation, regulates murine macrophage inflammatory responses. Proc Natl Acad Sci U S A 106: 15819–15824 .
- 27. Rogler CE, Levoci L, Ader T, Massimi A, Tchaikovskaya T, et al. (2009) MicroRNA-23b cluster microRNAs regulate transforming growth factor-beta/bone morphogenetic protein signaling and liver stem cell differentiation by targeting Smads. Hepatology 50: 575–584 .
- 28. Bavan L, Midwood K, Nanchahal J (2011) A New Avenue for Wound Healing Research. BioDrugs 25: 27–41 .
- 29. Zhu S, Pan W, Song X, Liu Y, Shao X, et al. (2012) The microRNA miR-23b suppresses IL-17-associated autoimmune inflammation by targeting TAB2, TAB3 and IKK-α. Nat Med 18: 1077–1086 .
- 30. Burton MJ, Ramadhani A, Weiss HA, Hu V, Massae P, et al. (2011) Active trachoma is associated with increased conjunctival expression of IL17A and profibrotic cytokines. Infect Immun 79: 4977–4983 .
- 31. Duisters RF, Tijsen AJ, Schroen B, Leenders JJ, Lentink V, et al. (2009) miR-133 and miR-30 regulate connective tissue growth factor: implications for a role of microRNAs in myocardial matrix remodeling. Circ Res 104: 170–8 .
- 32. Zhou H, Xu X, Xun Q, Yu D, Ling J, et al. (2012) microRNA-30c negatively regulates endometrial cancer cells by targeting metastasis-associated gene-1. Oncol Rep 27: 807–812 .
- 33. Kessler M, Zielecki J, Thieck O, Mollenkopf H-J, Fotopoulou C, et al. (2011) Chlamydia Trachomatis Disturbs Epithelial Tissue Homeostasis in Fallopian Tubes via Paracrine Wnt Signaling. Am J Pathol 180: 198–186 .
- 34. Schopman NCT, Heynen S, Haasnoot J, Berkhout B (2010) A miRNA-tRNA mix-up: tRNA origin of proposed miRNA. RNA Biol 7: 573–576 .
- 35. Cole C, Sobala A, Lu C, Thatcher SR, Bowman A, et al. (2009) Filtering of deep sequencing data reveals the existence of abundant Dicer-dependent small RNAs derived from tRNAs. RNA 15: 2147–2160 .
- 36. Thomas ME, Bueno M, Wirapati P, Lefort K, Dotto GP, et al. (2009) Impact of normalization on miRNA microarray expression profiling. RNA 15 (3) 493–501 .
- 37. Ach RA, Wang H, Curry B (2008) Measuring microRNAs: comparisons of microarray and quantitative PCR measurements, and of different total RNA prep methods. BMC Biotechnol 8: 69 .
- 38.
Rothman KJ (1986) Modern epidemiology. 2nd edition. Little Brown & Co.
- 39. Canales RD, Luo Y, Willey JC, Austermiller B, Barbacioru CC, et al. (2006) Evaluation of DNA microarray results with quantitative gene expression platforms. Nat Biotechnol 24: 1115–1122 .
- 40.
Git A, Dvinge H, Osborne M, Kutter C, Hadfield J, et al.. (2010) Systematic comparison of microarray profiling , real-time PCR, and next-generation sequencing technologies for measuring differential microRNA expression: 991–1006. doi:10.1261/rna.1947110.one.
- 41. VanGuilder HD, Vrana KE, Freeman WM (2008) Twenty-five years of quantitative PCR for gene expression analysis. Biotechniques 44: 619–626 .