Telomere Length, TERT, and miRNA Expression

It has been proposed that miRNAs are involved in the control of telomeres. We test that hypothesis by examining the association between miRNAs and telomere length (TL). Additionally, we evaluate if genetic variation in telomerase reverse transcriptase (TERT) is associated with miRNA expression levels. We use data from a population-based study of colorectal cancer (CRC), where we have previously shown associations between TL and TERT and CRC, to test associations between TL and miRNA expression and TERT and miRNA expression. To gain insight into functions of miRNAs associated with TERT we tested linear associations between miRNAs and their targeted gene mRNAs. An Agilent platform that contained information on over 2000 miRNAs was used. TL was measured using a multiplexed quantitative PCR (qPCR). RNAseq was used to assess gene expression. Our sample consisted of 1152 individuals with SNP data and miRNA expression data; 363 individuals with both TL and miRNA; and 148 individuals with miRNA and mRNA data. Thirty-three miRNAs were directly associated with TL after adjusting for age and sex (false discovery rate (FDR) of 0.05). TERT rs2736118 was associated with differences in miRNA expression between carcinoma and normal colonic mucosa for 75 miRNAs (FDR <0.05). Genes regulated by these miRNAs, as indicated by mRNA/miRNA associations, were associated with major signaling pathways beyond their TL-related functions, including PTEN, and PI3K/AKT signaling. Our data support a direct association between miRNAs and TL; differences in miRNA expression levels by TERT genotype were observed. Based on miRNA and targeted mRNA associations our data suggest that TERT is involved in non-TL-related functions by acting through altered miRNA expression.


Introduction
MicroRNAs (miRNA) are non-coding RNAs that function as gene regulators [1]. As such, they provide a mechanism whereby genes involved in various signaling pathways can be regulated simultaneously [2]. Both aging and cellular senescence are thought to be regulated in part by miRNAs [3,4]. It also has been suggested that miRNAs control the function of telomeres in cancer [5]. Telomeres, located at the end of chromosomes, shorten with every cell division; when they become too short the cell is unable to divide, leading to the induction of senescence and apoptosis, which are key features of normal cell function, but absent in tumor cells. Telomere length (TL) has been associated with colorectal cancer, with higher risk of disease for both long and short telomeres [6,7]. Telomerase, a ribonucleoprotein enzyme that is essential for the replication of chromosome ends, is low or non-expressing in normal human somatic cells, but has been shown to be expressed in the majority of tumor cells [8]. Telomerase reverse transcriptase (TERT) is essential for maintenance of telomere length by protecting chromosome ends. TERT has been proposed to regulate miRNAs by regulation of miRNA biogenesis, with diminished primary miRNA expression in TERT-reduced cells [9], and has been shown to be associated with both TL and cancer [6,[10][11][12][13][14]. TERT is thought to be one of the most critical, if not the rate-limiting step, in production of telomerase activity [8,15]. We have shown previously that TERT rs2853676 was associated with TL and that TL was associated with colorectal cancer [11]. Additionally TERT-CLPTM1L rs2853668 and TERT rs2736118 were associated with colorectal cancer risk in this study population [14].
In this paper we examine the associations between TL and miRNA expression in colorectal tissue. Additionally we test the hypothesis that variation in TERT influences miRNA expression levels. We focus on associations between miRNAs and TERT SNPs previously associated with colorectal cancer and examine miRNA expression in normal colonic mucosa tissue as well as the differential expression of miRNAs between carcinoma and normal colonic mucosa tissue.

Study population
The study population consisted of individuals previously enrolled in a study of Diet, Lifestyle, and Colon cancer at the University of Utah and the Kaiser Permanente Medical Research Program (KPMRP) [16] for whom both miRNA from carcinoma and adjacent normal colonic mucosa tissue were available. Study subjects included incident cases of colon cancer between the ages of 30 and 79 who were non-Hispanic white, Hispanic, or African American, and were able to provide a signed informed consent prior to participation in the study. The study was approved by the University of Utah Institutional Review Board for Human Subjects. miRNA processing RNA (miRNA) was extracted from formalin-fixed paraffin embedded tissues. We assessed slides and tumor blocks that were prepared over the duration of the study prior to the time of miRNA isolation to determine their suitability. The study pathologist reviewed slides to delineate carcinoma and normal colonic mucosa tissue. Cells were dissected from 1-4 sequential sections on aniline blue stained slides using an H&E slide for reference. Total RNA containing miRNA was extracted, isolated, and purified using the RecoverAll Total Nucleic Acid isolation kit (Ambion); RNA yields were determined using a NanoDrop spectrophotometer. 100 ng total RNA was labeled with Cy3 and hybridized to Agilent Human miRNA Microarray V19.0 and were scanned on an Agilent SureScan microarray scanner model G2600D. The Agilent Human microarray was generated using known miRNA sequence information compiled in the Sanger miRBase database v19.0. The microarray contains probes for 2006 unique human miR-NAs. The miRNA array contains 60,000 unique human sequences and averages 30 replicates per probe sequence. Data were extracted from the scanned image using Agilent Feature Extract software v.11.5.1.1. Data were required to pass stringent QC parameters established by Agilent that included tests for excessive background fluorescence, excessive variation among probe sequence replicates on the array, and measures of the total gene signal on the array to assess low signal. If samples failed to meet quality standards for any of these parameters, the sample was re-labeled, hybridized to arrays, and scanned. If a sample failed QC assessment a second time the sample was deemed to be of poor quality and the individual was excluded from downstream analysis. A 75 th percentile scaling was performed to normalize across arrays was done using preprocessCore [17] (www.bioconductor.org) to minimize differences that could be attributed to the array, amount of RNA, location on array, or other factors that could erroneously influence expression measurements. Testing of both reliability and comparative validity were done [18]. We showed a 98% repeatability in the Agilent platform between repeat samples. Additionally, the Agilent platform had fairly good agreement with Nanostring results and 100% agreement in expression values and fold change with qPCR results [19].

RNA-Seq Sequencing Library Preparation and Data Processing
Total RNA was available from 197 carcinoma and normal mucosa pairs. These samples were taken from the study subjects used for miRNA analysis and were extracted, isolated and purified in the same manner as previously described [20]. RNA library construction was done with the Illumina TruSeq Stranded Total RNA Sample Preparation Kit with Ribo-Zero. The samples were then fragmented and primed for cDNA synthesis, adapters were then ligated onto the cDNA, and the resulting samples were then amplified using PCR; the amplified library was then purified using Agencount AMPure XP beads. A more detailed description of the methods can be found in our previous work [21].
Sequencing was done using an Illumina TruSeq v3 single read flow cell and a 50 cycle single-read sequence run was performed on an Illumina HiSeq instrument. Reads were then aligned to a sequence database containing the human genome (build GRCh37/hg19, February 2009 from genome.ucsc.edu) and alignment was performed using novoalign v2.08.01. Python and a pysam library were used to calculate counts for each exon and UTR of the genes using a list of gene coordinates obtained from http://genome.ucsc.edu. We dropped features that were not expressed in our data or for which the expression was missing for the majority of samples. A more detailed description of the methods can be found in our previous work [21].
Telomere Length (TL) Telomere-related Gene Analysis TL was measured using a multiplexed quantitative PCR (qPCR) method previously described by Cawthon [22]. This method modified earlier qPCR methods in which telomere signals were measured separately from single copy gene signals in order to produce a T/S ratio corresponding to TL, which is proportional to the average telomere length in a cell. Therefore a larger T/S ratio represents a longer telomere length. The multiplexed PCR analysis uses a single dye and measures both the telomere signals and single copy gene signals in the same well. This is achieved by using CG clamps to stabilize the single copy gene giving it a higher melting point. The unique sequence is then amplified during a different cycle than the telomeric sequence, allowing for a single qPCR to determine the T/S ratio [6].
DNA was obtained from immortalized cell lines from study participants for the colon cancer study and from whole blood for the rectal cancer study. To assure the quality of telomere assay for these samples from different DNA origins we evaluated associations with TL separately for cases and controls from each study and as shown previous, results were similar for both study groups and we concluded that DNA source did not significantly alter TL length in this study. Two samples yielded a T/S ratio greater than three and were excluded. These samples were excluded because we concluded, through gel electrophoresis, that the DNA was degraded and therefore the T/S ratio was unreliable. We focused on Single Nucleotide Polymorphisms (SNPs) that were previously associated with TL and/or colon cancer [6,14]. We analyzed TERT-CLPTM1L rs2853668 which is 5.1kb upstream of TERT and has been associated with colon cancer in GWAS [23]. We also test TERT rs2736118 and TERT rs2853676 with miRNA expression.

Statistical Analysis
Our sample consisted of 1152 individuals with both SNP data and miRNA expression data; 363 individuals with both TL and miRNA expression data in the normal colonic mucosa; and 148 individuals who had both miRNA and mRNA data, and 176 cases in the colon cancer study and 184 cases in the rectal cancer study had both TL and SNP data.
Three major components of the analysis were conducted. First, we evaluated the linear associations between TL and those miRNAs for which at least 20% of the population had expression (N = 817 miRNA) by fitting a linear regression model to the log base 2 transformed expression levels; we adjusted for age and sex. P-values for each regression analysis were generated from the bootstrap method by creating a distribution of 1,000 β coefficients for each miRNA and evaluating H 0 : β = 0 vs. H 1 : β6 ¼0 using the boot package in R. We adjusted for multiple comparisons using a false discovery rate (FDR) level of 0.05 [24]. We standardized the slopes generated from the overall dataset in order to compare the results across the miRNA.
Second, we assigned inheritance models for TERT SNPs based on previous findings for colorectal cancer risk using either dominant or recessive models [14]. We compared log base 2 transformed expression levels across the selected genotype models using the significance analysis of microarrays (SAM) technique implemented in the R package siggenes; p-values were based on 1,000 permutations. We utilized an FDR level of 0.05 to determine which TERT SNP/ miRNA associations were significant after adjustment for multiple comparisons.
For both the TL/miRNA analysis and the TERT SNP/miRNA analysis, we evaluated miRNA expression in both normal colonic mucosa and for differentially expressed miRNAs (i.e. difference between carcinoma tissue and normal colonic mucosa). Looking at normal colonic mucosa allowed us to evaluate the effect of TL or SNP on overall miRNA expression. However, since an association with normal colonic mucosa could influence baseline level of the miRNA expressed and could potentially miss important changes in miRNA expression between carcinoma and normal colonic mucosa, we also evaluated miRNA differential expression. This allowed us to gain insight into regulation of miRNAs that could be more directly associated with cancer through their influence on carcinoma miRNA expression; by comparing to normal expression we were able to control for differences in tumor subsite.
One of the goals of the study was to identify non-telomere related mechanisms whereby TERT may be operating through altered miRNA expression levels. To accomplish this, we undertook the third part of the analysis, where we used miRTarBase (www.http://mirtarbase. mbc.nctu.edu.tw/) v6.0 (as of 09/15/2015) [25] to identify validated target genes for the significant miRNAs identified in our TERT SNP/miRNA analysis. We calculated the differential mRNA expression as the difference of the log base 2 of the RPKM (Reads per Kilobase per Million) for the carcinoma and normal mucosa tissues for the 4047 target genes identified by miR-TarBase. We then evaluated the 6309 combinations of miRNAs with their respective mRNA target genes. We used the bootstrap method as described above to assess any linear relationship between the differential mRNA and miRNA expression levels, adjusting for age and sex.
QIAGEN's Ingenuity Pathway Analysis (IPA) (QIAGEN Redwood City, www.qiagen.com/ ingenuity) [26] was then used to evaluate the genes that were dysregulated in conjunction with the miRNAs. We considered direct relationships from the IPA knowledge base of genes only, limited to experimentally verified and mammalian results only, and considered all mutations, data sources, and tissues. We report canonical pathways that were significant enriched with these genes after adjustment for multiple comparisons.

Results
The average age of study participants from the colon cancer study was slightly older than from the rectal cancer study ( Table 1). The distribution of age, sex, and AJCC stage was similar for those individuals included in the smaller TL dataset as with the larger SNP dataset.
Assessment of the linear relationship between TL and miRNA expression levels in for the 363 individuals with both TL and miRNA from normal colorectal mucosa showed that 34 miR-NAs were associated with TL after adjustment for age and sex ( Table 2). Thirty-three of the miRNAs were directly associated with TL, in that TL increased as miRNA expression level increased. Only miR-487b was inversely associated with TL, showing lower levels of miR-487b expression with longer TL length.
Evaluation of differences in miRNA expression levels between carcinoma tissue and normal colorectal mucosa showed that TERT rs2736118 was associated with differential miRNA expression based on genotype (Table 3). TERT rs2736118 was associated with 75 miRNAs with an FDR of <0.05. Individuals with the AA (homozygous dominant) genotype were more likely to have greater miRNA expression in the tumor than in the normal colorectal mucosa (61% of the time) while individuals with the AG/GG genotypes were more likely to have greater miRNA expression in the normal colorectal mucosa (84% of the time). Neither TERT-CLPT-MIL rs2853668 nor TERT rs2853676 significantly altered miRNA expression in normal colonic mucosa or miRNA expression between carcinoma and normal colorectal mucosa.
The 75 miRNAs dysregulated by TERT rs2736118 were associated with 4047 genes and 6309 miRNA/mRNA associations according to miRTarBase. To determine which of these genes might be associated with mRNA regulation in colonic tissue, we assessed each miRNA with its mRNA gene target in colonic carcinoma tissue and colonic normal mucosa using a linear regression. We identified 573 miRNA/mRNA associations p-value of <0.05 prior to adjustment for multiple comparisons. Of these, 150 significant miRNA/mRNA target associations remained when the FDR was 0.05 and 270 miRNA/mRNA associations remained with the FDR of <0.1 ( Table 4). The miRNAs and their targeted genes shown in Table 4 also indicate the direct of the association based on the beta coefficient. Of note is that the associations were limited to 13 miRNAs. All but five of these associations showed inverse associations between the miRNA and the mRNA. Using IPA to construct pathways enriched with the 270 genes with an FDR of <0.1, we identified 14 pathways (Fig 1). Although less than 10% of genes in each pathway were included hsa-miR-939-5p 9.06 0.14 < .0001 < .0001 Table 3. Associations between TERT rs2736118 and miRNA differential expression 1 with an FDR of <0.05.

Discussion
These data suggest a direct association between miRNAs and TL as well as between TERT rs2736118 and miRNA differential expression between carcinoma tissue and normal colonic mucosa. While TERT expression is low in most somatic cells, endothelial cells represent an established exception and associations observed for differential expression between carcinoma and normal mucosa is a likely extension of that observation. Our results add support to the hypothesis that miRNAs are associated with TL and that genetic variation in telomere-related genes influence miRNA expression levels. Given the number of miRNAs altered by TERT rs2736118, our ability to show associations between altered miRNAs and their targeted mRNA, and the pathways represented by these target genes, we believe that our data support the hypothesis that TERT has functions beyond that of regulating TL. Focusing only on those miR-NAs that were statistically associated with target mRNA expression in colon tissue, we identified 14 canonical pathways that were associated significantly with these TERT-regulated miRNAs after adjustment for multiple comparison, 13 of which were not associated with telomerase.
There have been few reports of miRNAs associated with TL that have looked at actual miRNA expression levels in conjunction with TL; our sample size is by far the largest study to date to examine miRNA expression with TL. One study with only two samples compared shortened telomeres to intact telomeres and showed that 47 miRNAs had over a twofold difference in expression levels [27]. MiR-155, miR-138, and miR-34a,b,c also have been hypothesized as being associated with TL [3,5]; none of these miRNAs were associated in a linear fashion with TL in our study. The miRNAs associated with TL in our data for the most part increased TL with increased miRNA expression. Differences in association between our study and others may in part be from the methods of analysis. We assessed a linear association, while   other studies have often used categorical variables and looked at either short or long TL, which could account for some of the differences observed. A limitation of the study is that we are unable to confirm our observed associations between TL and miRNAs in bio-functional studies; we encourage others to do so. However, support for our findings come from several sources. First, we have previously shown that 26 of the 34 miR-NAs associated with TL were differentially expressed between colorectal carcinoma and normal colonic mucosa and 28 of the 34 were dysregulated between carcinoma and adenomatous polyp tissue; three miRNAs were dysregulated between normal colonic mucosa and adenomatous polyp tissue [18]. Interestingly, the only miRNA that was directly associated with TL, miR-487b, was the only miRNA that was not associated with dysregulation between carcinoma, adenoma, and/or normal mucosa. Additionally, using bioinformatics tools to identify the target genes for the dysregulated miRNAs associated with TL, we observed that 12 of these miRNAs associated with TL target genes (summary in Table 5) are components of the alternative lengthening of telomeres pathway (ALT) and Shelterin [28]. Shelterin is a conserved protein component on telomeres that serves as a functional framework of telomere chromatin pathways. ALT is a telomerase-independent pathway of lengthening telomeres in human cells. Given that there is no overlap between TERT SNP-associated miRNAs and TL-associated miRNAs, and that a large proportion of these miRNAs target genes in the ALT pathway, it is possible that miRNAs regulate TL through this alternative mechanism. As TL has been shown to be associated with CRC, the associations between TL and miRNA expression may unique clinical significance.
We also observed that genetic variation in TERT, which we have previously reported as being associated with TL and cancer [6,14], was associated with miRNA expression levels. Many of these miRNAs have been shown to be dysregulated in CRC [18] and 26 of these miR-NAs target ALT pathway genes ( Table 5). One of the pathways enriched by genes targeted by dysregulated miRNAs is 'telomerase signaling' . Since telomerase is generally not expressed in non-tumor somatic tissue but is expressed in the majority of tumors, it is logical that TERT rs2836118 was not associated with miRNA expression in normal colonic mucosa, but only in differential miRNA expression between carcinoma and normal colorectal mucosa. Telomerase is often reactivated by TERT upregulation, which gives tumor cells their immortal quality, and this is thought to be a key step in the adenoma-carcinoma transition in human tumorigenesis [29].
As stated previously, we identified miRNAs that were differentially expressed between normal colonic mucosa and adenomatous polyp tissues as well as between adenoma and carcinoma tissues, which show the changes in miRNA expression as tissue progresses from normal to adenomatous to cancerous [18]. Both hsa-miR-520e and hsa-miR-583 were differentially expressed between normal colonic mucosa and adenomatous polyp tissue [18]; in this study we show that they also regulate mRNA expression. These two miRNAs regulated 40 genes, including ITGA2, BMP8B, and PDPK1, which were found in numerous pathways in IPA, including PTEN Signaling, BMP signaling pathway, and Basal Cell Carcinoma Signaling. Hsa-miR-520e was associated with CBX5, FANCA and NR2F2, and hsa-miR-583 was associated with CBX5 and CDKN1A, all of which are part of the ALT pathway. We also identified six miRNAs that were downregulated with the AG/GG (heterozygous-homozygous recessive) genotype of rs2736118 that we previously found were dysregulated between adenomatous and carcinoma tissue; these miRNAs were: hsa-miR-206, hsa-miR-4510, hsa-miR-4674, hsa-miR-639, hsa-miR-6718-5p, and hsa-miR-708-5p. This supports the hypothesis that TERT regulates miRNA expression and this influences the formation of adenomatous tissue and contributes to tumorigenesis. The association of TERT rs2736118 genotype with miRNA expression levels may provide additional insight into pathways in which TERT is involved. We therefore examined genes and pathways regulated by the miRNAs that were influenced by TERT rs2736118. We observed that genetic variation in TERT was associated with dysregulation of 75 miRNAs, which in turn could influence translation and function of over 6000 genes. While many genes are regulated by these miRNAs, it is not clear which of these are related to colorectal cancer. Evaluation of over 6000 genes leads to a very non-specific assessment of potential functionality. To better focus on functionality we attempted to cluster genes together into what could be potentially a more common function and we also assessed seed regions. Using a cluster analysis approach, we still had over fifty pathways that also lacked specificity to colorectal cancer. Focusing on seed regions, we were left with 73 seed regions and did not minimize the potential number of genes that could be associated with TERT. Thus, we took an approach that focused only on those genes where the gene targeted by the miRNA was linearly associated with the corresponding mRNA expression in colorectal cancer. This reduced the number of genes and pathways considerably and was specific to colorectal cancer. However, while we gained specificity, we undoubtedly were limited in the target genes we assessed. We believe that the pathway information obtained is meaningful, but most certainly not complete. These data do however support the hypothesis that TERT has non-TL-related pathways [29,30] that may be acting through altered miRNA expression.
In addition to its established role in telomere maintenance, it has been proposed that TERT operates through other pathways. Further examples of TERT involvement in non-TL related pathways include both NFκB complex and the Wnt/β-catenin pathways, which have been hypothesized as being telomerase-targeted pathways as both are pathways related to inflammation [5]; these pathways were both represented by our target genes in IPA although they are not shown in Fig 1. Most notably, the NFκB-signaling pathway activation has been suggested as being related to TERT [31]. TERT levels also have been associated with differential expression of genes involved in development/morphogenesis, signal transduction, and cell adhesion-signaling pathways [32,33]. Our data also suggest that associations between TERT and colorectal cancer, possibly mediated by miRNAs, involve numerous pathways other than telomerase. These pathways involve tumor suppressor genes (such as PTEN signaling), apoptosis (PI3K/AKT signaling, NGF signaling, and FLT3 signaling), angiogenesis (Erythropoietin signaling), immune regulation and response (Prolactin signaling and Aryl Hydrocarbon Receptor), cell migration and proliferation (Neuregulin signaling), and growth factor activity (NGF signaling and EIF4 and p70S6K signaling). Another enriched pathway, Tight Junction Signaling, is important in cellular communication and facilities crosstalk across several of these signaling pathways including PI3K/AKT and PTEN [34].
Additionally, specific genes such as STAT5A, MELK, MAPK1 and MAP3K9, RELA (part of the NFKB1 complex), CDK1, tumor necrosis factor genes, IRF1, and TGFβ-related genes were associated with miRNAs differentially expressed by TERT genotype. A number of translocations associated with cancer involve the Ig heavy chain locus on chromosome 14. For example t(8;14) brings c-MYC from chromosome 8 to the IgH gene locus causing an overexpression of c-MYC in Burkitt lymphoma. A t(18;14) translocation moving Bcl-2 from chromosome 18 to the IgH locus on 14 is involved in follicular lymphoma. Also CCND1 is involved in a t (11;14) translocation bringing cyclin d1 to the IgH resulting in cell cycle dysregulation in Mantle cell lymphoma. Studies suggest that these cancers are associated with telomere length [35][36][37][38] as well as with miRNA levels [39]. We show that TERT regulation of miRNAs alters Bcl-2 and an enriched pathway associated with the dysregulated genes is chronic myeloid leukemia, which has overlapping genes with some of the lymphomas previously studied. Our data support the hypothesis that TERT is involved in disease processes that include non-TL-related pathways.
One of the major strengths of our study was the availability of miRNA expression data, SNP data, and mRNA data on a large population. We have previously analyzed these SNPs with both TL and colon cancer, which helps facilitate the interpretation of our data. Additionally, we applied a FDR to guide our interpretation of meaningful results. Larger samples would have enabled us to split the sample and replicate our findings with a second data set. While we are unable to replicate these findings in a second data set, we encourage others with similar data to undertake such analysis. There are many available tools to assist in identifying verified mRNA targets for miRNA. In this study we used several bioinformatics tools to better understand function of these miRNAs, however all bioinformatics tools have limitations as to the extent of the knowledge and the specificity of the information to specific miRNAs. To have a more focused assessment of pathways, we restricted our bioinformatics assessment to only those targeted genes significantly associated with miRNA. Although functionality data of specific miRNA and TERT associations are not available or feasible from data available to us, we believe that the information provided gives direction for more laboratory-based studies that can directly test these associations.

Conclusions
In conclusion, our data suggest that miRNAs are involved in regulating TL. TERT appears to influence carcinoma/normal mucosa differential expression. Given the number of miRNAs associated with TL, TERT rs2736118 and ALT genes, our data also support the hypothesis that telomere-related genes, in addition to affecting TL, impact non-TL-related functions that can importantly influence cancer risk.
Dr. Bette Caan and the Kaiser Permanente Medical Research Program for sample contributions, Brett Milash and the Bioinformatics Shared Resource of the Huntsman Cancer Institute and University of Utah for miRNA and mRNA bioinformatics data processing, Sandra Edwards for the data collection and management, Wade Samowitz for slide review, and Erica Wolff and Michael Hoffman for miRNA processing.