- Open Access
The TGFβ-signaling pathway and colorectal cancer: associations between dysregulated genes and miRNAs
Journal of Translational Medicinevolume 16, Article number: 191 (2018)
The TGFβ-signaling pathway plays an important role in the pathogenesis of colorectal cancer (CRC). Loss of function of several genes within this pathway, such as bone morphogenetic proteins (BMPs) have been seen as key events in CRC progression.
In this study we comprehensively evaluate differential gene expression (RNASeq) of 81 genes in the TGFβ-signaling pathway and evaluate how dysregulated genes are associated with miRNA expression (Agilent Human miRNA Microarray V19.0). We utilize paired carcinoma and normal tissue from 217 CRC cases. We evaluate the associations between differentially expressed genes and miRNAs and sex, age, disease stage, and survival months.
Thirteen genes were significantly downregulated and 14 were significantly upregulated after considering fold change (FC) of > 1.50 or < 0.67 and multiple comparison adjustment. Bone morphogenetic protein genes BMP5, BMP6, and BMP2 and growth differentiation factor GDF7 were downregulated. BMP4, BMP7, INHBA (Inhibin beta A), TGFBR1, TGFB2, TGIF1, TGIF2, and TFDP1 were upregulated. In general, genes with the greatest dysregulation, such as BMP5 (FC 0.17, BMP6 (FC 0.25), BMP2 (FC 0.32), CDKN2B (FC 0.32), MYC (FC 3.70), BMP7 (FC 4.17), and INHBA (FC 9.34) showed dysregulation in the majority of the population (84.3, 77.4, 81.1, 80.2, 82.0, 51.2, and 75.1% respectively). Four genes, TGFBR2, ID4, ID1, and PITX2, were un-associated or slightly upregulated in microsatellite-stable (MSS) tumors while downregulated in microsatellite-unstable (MSI) tumors. Eight dysregulated genes were associated with miRNA differential expression. E2F5 and THBS1 were associated with one or two miRNAs; RBL1, TGFBR1, TGIF2, and INHBA were associated with seven or more miRNAs with multiple seed-region matches. Evaluation of the joint effects of mRNA:miRNA identified interactions that were stronger in more advanced disease stages and varied by survival months.
These data support an interaction between miRNAs and genes in the TGFβ-signaling pathway in association with CRC risk. These interactions are associated with unique clinical characteristics that may provide targets for further investigations.
The TGFβ-signaling pathway is important in the tumorigenesis of colorectal cancer (CRC) . This pathway is a regulator of cellular proliferation, differentiation, apoptosis, and extracellular matrix remodeling, and also is involved in angiogenesis and inflammation [2, 3]. It has been previously reported that components of the TGFβ-signaling pathway are mutated in 27% of non-hypermutated tumors and 87% of hypermutated tumors , and that inactivation of the pathway is a common event in CRC tumorigenesis . The TGFβ family of cytokine genes includes three isoforms of TGF-β, TGF-β1, TGF-β2, and TGF-β3, the type I receptors (TβR1 and ALK1) and a type II receptor (TβRII) . Other components of the TGFβ-signaling pathway include SMAD genes, key intracellular mediators of the transcriptional responses to TGF-β  and bone morphogenetic proteins (BMP), which trigger a SMAD-signaling cascade that is linked to cell proliferation and cellular growth [7, 8]. Growth differentiation factors (GDF)  and their receptors are components of the TGFβ superfamily as are Activin/inhibin and their receptors . BMP ligands bind to type 1 [BMPR1A, BMPR1B, Activin A receptor type 1 (ACVR1), and Activin receptor-like kinase 1 (ACVRL1)] and type 2 receptors [BMPR2, Activin A receptor type IIA (ACVR2A) and type IIB (ACVR2B)] . While both type 1 and 2 receptors are needed for BMP signaling, type I receptors bind with a higher affinity than type II receptors. BMPR1A has been reported as being inactivated in some studies focusing on familial syndromes such as mixed polyposis syndrome and familial juvenile polyposis [11,12,13]. Other factors that regulate TGFβ and its receptors are BAMBI (BMP and activing membrane bound inhibitor) , THBS1 (thrombospondin 1 also known as TSP1) , LEFTY (left–right determination factor, Factor 2 also known as TGFβ4) , and FST (Follistatin) .
It has been hypothesized that miRNAs, small, non-protein-coding RNA molecules involved in the regulation of gene expression either by post-transcriptionally suppressing mRNA translation or by causing mRNA degradation [18,19,20,21,22,23], may work with the TGFβ-signaling pathway to mediate cell growth and promote tumorigenesis . MiRNAs have been linked to TGFβ-signaling pathway in a variety of diseases, usually in studies which only examined few miRNAs and genes within the signaling pathway. For instance, miR-181a has been shown to have its expression altered by TGFβ [24, 25], miR-494 and miR-126-5p have been linked to BMP4 in the regulation of angiogenesis ; miR-98 with TBSP1 expression in asthma ; miR-590-5p with downregulation of TGFβ signaling in cardiosphere-derived stem cells ; the miR-17-92 cluster with cell proliferation in mouse mesenchymal cells ; and miR-140-5p with tumor growth . MiRNAs have also been linked with various disease processes in CRC. For example, miR-200a-3p regulates epithelial mesenchymal transition in CRC ; and miR-21-5p is upregulated in CRC, while its inhibition can hinder CRC tumor growth . Although the TGFβ-signaling pathway appears to be important for CRC, few studies have looked at this pathway with miRNAs. We have previously examined miRNA expression with SNPs in 21 genes associated with the TGFβ-signaling pathway and found that expression of several miRNAs varied by SNPs in TGFβ1 in normal mucosa . Additionally, other groups have found that specific miRNAs play a significant role in CRC via the TGFβ pathway. MiR-193b has been shown to play an important role in CRC by promoting cellular proliferation via SMAD3 ; miRs-203, 181d and 182 regulate genes, including TGF-β2 and RUNX2, involved in CRC proliferation, differentiation, and invasion ; and miR-34a mediates oxaliplatin resistance in CRC via the TGFβ/SMAD4 pathway .
In this study, we comprehensively examined expression differences between carcinoma and normal mucosa with all genes identified in the TGFβ-signaling pathway. For those genes that had significant differences in expression when considering level of expression (> 1.50 fold change (FC) or < 0.67 FC), we assessed their association with expression levels of over 800 miRNAs commonly expressed in CRC. We evaluated associations for overall CRC as well as for specific tumor molecular phenotype, i.e. tumors with or without mismatch repair deficiency, and how these associations relate to clinical features such as disease stage and survival months. Our goal is to obtain a better understanding of how miRNAs relate to the TGFβ-signaling pathway in CRC carcinogenesis. We believe that our epidemiological approach, using population-based cases, is an initial step in identifying associations that can be further examined in targeted laboratory studies with the goal of identification of targeted therapeutics.
Participants come from two population-based case–control studies that included all incident colon and rectal cancer diagnosed between 30 and 79 years of age in Utah or at Kaiser Permanente in Northern California (KPNC). Diagnosis was verified by tumor registry data as a first primary adenocarcinoma of the colon or rectum with a diagnosis date between October 1991 and September 1994 (colon study) and between May 1997 and May 2001 (rectal study) . Participants were non-Hispanic white, Hispanic, or black for the colon cancer study; the rectal cancer study also included people who self-reported being Asian [37, 38]. The Institutional Review Boards at the University of Utah and at KPNC approved the study.
Ethics, consent, and permissions
Study participants signed informed consent. The Institutional Review Boards at the University of Utah and at KPNC approved this study. Individual-level participant data are not reported.
Formalin-fixed paraffin embedded tissue from either surgery or the initial biopsy was used to extract RNA. RNA was extracted, isolated and purified from carcinoma tissue and adjacent normal mucosa as previously described . Differences in RNA quality were not observed based on age of the tissue; all samples were of a high quality.
mRNA: RNA-Seq Sequencing Library Preparation and Data Processing
RNA from 245 colorectal carcinoma and normal mucosa pairs was chosen for sequencing based on availability of RNA and high quality miRNA data; 217 pairs passed quality control (QC) and were used in these analyses. 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 . Illumina TruSeq v3 single read flow cell and a 50 cycle single-read sequence run was performed on an Illumina HiSeq instrument. Reads were 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. Total gene counts were calculated for each exon and UTR of the genes using a list of gene coordinates obtained from http://genome.ucsc.edu. We disregarded genes that were not expressed in our RNA-Seq data or for which the expression was missing for the majority of samples .
The Agilent Human miRNA Microarray V19.0 was used. 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. Samples failing to meet quality standards were re-labeled, hybridized to arrays, and re-scanned. If a sample failed QC assessment a second time, the sample was excluded from analysis. The repeatability associated with this microarray was extremely high (r = 0.98) ; comparison of miRNA expression levels obtained from the Agilent microarray to those obtained from qPCR had an agreement of 100% in terms of directionality of findings and the FCs were almost identical . To normalize differences in miRNA expression that could be attributed to the array, amount of RNA, location on array, or factors that could erroneously influence miRNA expression levels, total gene signal was normalized by multiplying each sample by a scaling factor which was the median of the 75th percentile of all the samples divided by the individual 75th percentile of each sample .
TGFβ-signaling pathway genes
The Kyoto encyclopedia of genes and genomes (KEGG) (www.genome.jp/kegg-bin/show_pathway?hsa04350) pathway map for TGFβ-signaling was used to identify genes within this pathway. Using this pathway map, we identified 84 genes, 81 of which had sufficient expression in CRC tissue for inclusion in the study (Additional file 1: Table S1).
We utilized a negative binomial mixed effects model in SAS (accounting for paired carcinoma/normal status) to determine genes in the TGFβ-signaling pathway that had a significant difference in expression between individually paired colorectal carcinoma and normal mucosa (i.e. differentially expressed). In this test we offset the overall exposure as the log expression of all identified protein-coding genes (n = 17,461). The Benjamini and Hochberg  procedure was used to control the false discovery rate (FDR) using a value of 0.05 or less. We calculated FC using the population means to give an estimate of the magnitude of difference in expression to identify mRNA and miRNAs. Population means were used because in some instances there was no expression in either the tumor or normal, making individual FCs impossible to calculate. Utilization of the population means we believe gave the best estimate of FC in the population. We calculated level of expression of each gene by dividing the total expression for that gene in an individual by the total expression of all protein-coding genes per million transcripts (RPMPCG or reads per million protein-coding genes). We focused on those genes with FC of > 1.50 or < 0.67 for analysis with differential miRNA expression since these levels of FC may have a greater biological significance than FCs closer to one. A FC of greater than one indicates a positive differential expression (i.e. up-regulated in carcinoma), while a FC between zero and one indicates a negative differential expression (i.e. down-regulated in carcinoma). We also report at the population level the percentage of the population with FC > 1.50 or < 0.67 to provide an estimate of the prevalence of these changes in the population.
We evaluated clinical indicators of age, sex, AJCC disease stage, and survival months with those genes that were dysregulated using Spearman correlation coefficients. Additionally we evaluated miRNA and mRNA associations for specific categories of sex, age (< 55, 55–65, and > 65), AJCC disease stage, and survival months (24, 24–60, and > 60 months).
There were 814 miRNAs expressed in greater than 20% of normal colorectal mucosa that were analyzed; differential expression was calculated using subject-level paired data as the expression in the carcinoma tissue minus the expression in the normal mucosa. In these analyses, we fit a least squares linear regression model to the RPMPCG differential expression levels and miRNA differential expression levels. P-values were generated using the bootstrap method by creating a distribution of 10,000 F statistics derived by resampling the residuals from the null hypothesis model of no association between gene expression and miRNA expression using the boot package in R. Linear models were adjusted for age and sex. Multiplicity adjustments for gene/miRNA associations were made at the gene level using the FDR by Benjamini and Hochberg .
We analyzed miRNAs and targeted mRNAs for seed region matches. The mRNA 3′ UTR FASTA as well as the seed region sequence of the associated miRNA were analyzed to determine seed region pairings between miRNA and mRNA. MiRNA seed regions were calculated as described in our previous work ; we calculated and included seeds of six, seven, and eight nucleotides in length. Our hypothesis is that a seed-region match would increase the likelihood that identified genes associated with a specific miRNA were more likely to have a direct association (as indicated by a negative beta coefficient) given a higher propensity for binding and thus mRNA degradation. We used mRNA FASTA sequences generated from both GRCh37 and GRCh38 Homo sapiens alignments, using UCSC Table Browser (https://genome.ucsc.edu/cgi-bin/hgTables) . We downloaded FASTA sequences that matched our Ensembl IDs and had a consensus coding sequences available. Analyses were done using scripts in R 3.2.3 and in perl 5.018002.
Utah miRNA data are available through GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE115513). Other data will be shared in accord with the signed consent and approved IRB studies. Contact study authors for data access.
Colon cancer cases comprised 77.9% of the study population (Table 1). The majority of cases were male (54.4%), non-Hispanic white (74.2%) and were diagnosed with a microsatellite stable (MSS) tumor. The mean age at diagnosis was 64.8 years.
Gene expression results
Of the 81 genes analyzed, 27 showed statistically significant differences in expression between carcinoma and normal mucosa with a more meaningful FC when tumor sub-site was not considered (Table 2). Of these 27 genes, 13 were downregulated and 14 were up-regulated. Of the downregulated genes, three were part of the BMP family (BMP5, BMP6, and BMP2) and two were part of the growth differentiation factor family (GDF7 and GDF6). INHBA (Inhibin beta A), which encodes for a subunit of activin and inhibin, was strongly upregulated, while the gene encoding activin type 1 receptor, ACVR1C, was downregulated. Two BMPs (BMP4 and BMP7) were upregulated. TGFBR1, TGFB2, TGIF1, TGIF2, and TFDP1 were upregulated. Other downregulated genes were: IFNG and CDKN2B (FC 0.32), AMHR2 and LEFTY2 (FC 0.43), LEFTY1 (FC 0.51), FST (FC 0.59), and THBS1 (FC 0.66). Other upregulated genes were: E2F5 (FC 1.52), AMH (FC 1.74), BAMBI (FC 1.85), INHBB (FC) 1.91), RBL1 (FC 2.30), and MYC (FC 3.70). While not all individuals had upregulated or downregulated expression in expression, when considering the entire population these genes were significantly up or downregulated. The stronger the FC in terms of being downregulated or upregulated the greater likelihood that more individuals in the population had that tumor characteristic. However, there are exceptions to this. For instance, LEFTY2 and AMHR2 both had an overall FC of 0.43, yet only 34.1 and 18.4% of individuals in the population had a FC of < 0.67 for these genes respectively. This illustrates the variability of gene expression in the population, but also suggest that some genes such as BMP5 BMP6, and BMP2 may be better markers of CRC given that 84.3, 77.4, and 81.1% of the population have significant downregulated expression. Figure 1 shows the KEGG TGFβ-Signaling Pathway and those genes that are up and downregulated within the pathway.
A small number of genes had either slightly stronger associations (as indicated by FC) for the MSS (Additional file 1: Table S2) or MSI (Additional file 1: Table S3) phenotype. Associations were most similar for MSS tumors and overall tumors; one gene (THBS1) was associated overall (FC 0.66) and for MSS-specific tumors (FC 0.68) that would not have merited inclusion by our cutpoints for FC. AMH (anti-Mullerian Hormone) was upregulated in overall tumors (FC 1.74) to a greater effect than for MSS tumors (FC 1.48). There were considerably more differences observed between MSI tumors and overall tumors. Seven genes were downregulated in MSI tumors with a much greater effect than for MSS tumors. Two of these genes were much more strongly downregulated in MSI tumors (ID3 FCMSS 0.86 and FCMSII 0.52; THBS1 FCMSS 0.68 and FCMSI 0.58). The other five genes were slightly upregulated in MSS tumors while downregulated in MSI tumors (INHBE FCMSS 1.06 and FCMSI 0.52; TGFBR2 FCMSS 1.18 and FCMSI 0.61; ID4 FCMSS 1.26 and FCMSI 0.62; ID1 FCMSS 1.30 and FCMSI 0.38; PITX2 FCMSS 2.36 and FCMSI 0.52). AMH was much more strongly upregulated for MSI tumors (FC 3.41).
For those genes that had significant differential expression with a FC of < 0.67 or > 1.50 we further evaluated their associations with sex, age, disease stage and survival months (Table 3). There were few correlations with any of these factors that we consider strong, although correlations > 0.134 were statistically significant with p values < 0.05. The strongest correlations for these factors were: INHBB with sex (r = − 0.175) indicating slightly stronger CDKN2B expression levels among women; BMP2 with age (r = 0.224) indicating less expression with increasing age; CDKN2B with disease stage (r = 0.197) indicating greater correlation with more advanced disease stage; and BMP4 with survival (r = − 0.174) indicating worse survival with greater BMP4 expression.
mRNA and miRNA associations
Eleven of the dysregulated genes were associated with miRNA differential expression (Table 4). All of the associated genes were upregulated except for THBS1. While two of the genes were only associated with a few miRNAs, E2F5 with miR-17-5p and miR-20a-5p, and THBS1 with miR-145-5p, the others were associated with seven or more miRNAs. The miRNAs associated with E2F5 had seed-region matches with this gene. RBL1 was associated with 11 miRNAs (five with a seed-region match), TGFBR1 with 11 miRNAs (nine with a seed-region match), TGIF2 with 26 miRNAs (eight with a seed-region match), INHBA with seven miRNAs (two with a seed-region match), MYC with 12 miRNAs, and TFDP1 with seven miRNAs. Of the miRNAs and mRNAs with a seed-region match four had inverse associations as indicated by a negative beta coefficient (TGFBR1 with miR-2117 and miR-6071; TGIF2 with miR-375 and miR-4749-3p). The FC of the differentially expressed miRNAs varied less in the population than the differential expression of the mRNAs. Although variability in the population still existed, when a population FC of > 1.50 or < 0.67 was observed for the miRNA, a much larger percentage of the population had a significant up or downregulated miRNA. The stronger the population FC, the larger percentage of individuals had a similar FC. Figure 1 displays the miRNAs with their associated mRNAs in the KEGG TGFβ-signaling pathway. Figure 2 highlights key components of the pathway, including those genes acting as extracellular factors (Fig. 2a), those acting as membrane receptors (Fig. 2b) and those acting as nuclear factors (Fig. 2c) that were associated miRNAs. Seed region matches are shown with a –| between the miRNA and mRNA, while interactions without seed matches are shown with a line.
Several miRNAs were associated with multiple genes (Table 5). Both miR-17-5p and miR-20a-5p were associated with the same three genes (RPBL1, TGIF2, and E2F5) with seed-region matches. MiR-199a-5p and miR-21-5p both had seed region matches to TGFBR1 and INHBA. MiR-93-5p had seed-region matches to both RBL1 and TGIF2. Both miR-1246 and miR-3651 were associated with RBL1, TGIF2, and MYC without seed-region matches. MiR-19b-3p and miR-20b-5p were associated with TGIF2, MYC, and TFDP1 without seed-region matches. MiR-501-3p and miR-663a were associated with MYC. MiR-92-3p was associated with the most genes without a seed-region match, RBL1, TGIF2, MYC, and TFDP1.
The combined mRNA:miRNA associations differed somewhat by category of sex, age, disease stage and months of survival (Table 6). While differences in the associations between the association of miRNA and mRNA did not differ for many associations by sex, it is noteworthy that for most of the miRNAs associated with TGIF2, correlations were twice as strong among men as among women. Most associations between RBL1, TGIF2, and TFDP1 and miRNAs tended to be much stronger for individuals < 55 years of age. Looking at disease stage, more differences in association existed for disease stage II, however several mRNA:miRNA associations were strongest for stage IV. Associations that were strongest for stage IV disease were: RBL1 with miR-17-5p, miRp19b-3p, miR-20a-5p, miR-25-3p, and miR-3651; TGFBR1 with miR-1203, miR-2117, mir-21-5p, miR-23a-3p, miR-24-3p, and miR-331-3p; TGIF2 with miR-106b-5p, miR-1246, miR-1291, miR-130b-3p, miR-17-5p, miR-196a-5p, miR-19b-3p, miR-20a-5p, miR-20b,5p, miR-221-3p, miR-25-3p, miR-361-5p, miR-3651, miR-424-3p, miR-425-5p, miR-4749-3p, miR-501-3p, miR-583, miR-663b, miR-92a-3p, and miR-93-5p; and TDFP1 with miR-145-5p, miR-17-5p, miR-19b-3p, miR-20a-5p, and miR-20b-5p. Several mRNA:miRNA associations were strongest among those with the fewest survival months including TGFBR1 with miR-331-3p and miR-21-5p; TGIF2 with miR-1291, mir-130b-3p, miR-17-5p, miR-196a-5p, miR-19b-3p, miR-20a-5p, miR-20b-5p, miR-221-3p, miR-25-3p, miR-361-5p, miR-424-3p, miR-483-3p, miR-501-3p, miR-583, miR-663b, miR-92a-3p, and miR-93-5p; and TFDP1 with miR-17-5p, miR-20a-5p, and miR-92a-3p. On the other hand, some mRNA:miRNA associations were strongest among those with better survival, such as RBL1 with miR-1246, miR-196b-5p, and miR-3651; and INHBA with miR-199a-3p, miR-199a-5p, miR-199b-5p, miR-214-3p, and miR-934.
The TGFβ pathway is one of the most important pathways in the development of CRC . In normal cells the TGFβ pathway plays an important role in suppressing growth and tumorigenesis; however, as cancer progresses the TGFβ pathway promotes epithelial-mesenchymal transition, invasion and metastasis . In this study, of the 27 dysregulated genes in this pathway, 13 genes were significantly downregulated. The genes that were downregulated included BMP5, BMP6 and BMP2, in addition to GDF7 and GDF6. We have previously shown an association between BMP2 SNPs and CRC risk . BMP5 and BMP6 play an important role in various cancers; both have been shown to induce apoptosis . BMP6 has previously been shown to be downregulated in metastatic breast cancer . Our data suggests these genes play an important inhibitory role in CRC tumorigenesis as well, and therefore, are downregulated similar to other cancers.
Growth differentiation factors are also known to have inhibitory effects on growth in human cancers, and both GDF9a and GDF9b have been associated with breast cancer . We found that both GDF6 and GDF7 were downregulated in CRC tissues. While we could find no earlier documentation of these genes having an association with CRC risk, GDF7 is a ligand in the BMP pathway that has been shown to regulate the Hedgehog and Wnt pathways and has been associated with the pro-inflammatory phenotypes in many diseases including Barrett’s Esophagus  and esophageal adenocarcinoma . Additionally, the role of Wnt in CRC has been established .
The genes that were upregulated in our study included BMP4, BMP7, TGFBR1, TGFB2, TGIF1, TGIF2, and TFDP1. BMP4 is known to regulate the SMAD 1/5 signal transduction pathway and plays an important role in angiogenesis . Our data suggest that in CRC, multiple BMPs have altered expression levels. However, it is interesting to note that downregulated BMP5, BMP6, and BMP2 were some of the most consistently downregulated genes in the population, with approximately 80% of cases having expression of these genes downregulated in tumors; BMP4 and BMP7 were upregulated in over 50% of CRC. BMPs are members of a pleiotropic family of growth factors whose signaling has been shown to protect against colonic polyposis with known mutations that can contribute to the development of the tumor microenvironment . Loss of BMP signals has been cited as one of the two main genetic alterations leading to CRC; disrupted BMP signaling allows tumor growth and expansion . BMPs have been shown to upregulate various cytokines important in cancer development and metastasis, and BMP expression levels have been associated with survival and prognosis in various cancers [54, 55]. Our data support these earlier findings of the importance of BMPs in the development and progression of CRC, in that BMP7, BMP5, BMP4, BMP2, and BMP6 were all significantly differentially expressed in CRC tumors and BMP4 differential expression was inversely associated with survival months. Moreover, they suggest that BMPs play very diverse roles in regulating the TGFβ pathway. Meaning that through their diverse group of downstream ligands and targets, some members of the BMP family can play an important inhibitory effect in tumorigenesis, and therefore are downregulated in CRC, while others play a stimulatory effect in tumorigenesis, and therefore must be upregulated in CRC.
In addition to evaluating the role of the TGFβ pathway in overall CRC risk, we evaluated gene expression levels for associations with mismatch repair proficiency and deficiency. Our data suggest that several genes in this pathway have significant associations with these phenotypes. Notably, THBS1 was associated with MSI-specific tumors. Members of the thrombospondin family previously have been associated with differential expression in different CRC phenotypes; THBS4 has been associated with CIMP high tumors, of which a subset are mismatch repair deficient . THBS1 has previously been associated with the development of the tumor micro-environment and angiogenesis  and metastasis . AMH was upregulated in all tumors, but had lower expression levels in MSS tumors and the highest expression levels in MSI tumors. AMH is selectively expressed in epithelial cells versus mesenchymal cells and a loss of AMH is known to induce EMT in various cancers . AMH signaling induces downstream phosphorylation of SMAD 1/5/8 . Notably SMADs 1/5/8 also are known to be downstream elements in the JNK and ERK pathways . TGFβ is a known activator of the ERK, JNK and p38 MAPK pathways , and our data suggests that AMH is a possible mediator of this cross-talk between pathways. We have previously shown that both ERK and JNK dependent pathways are significantly associated with CIMP, MSI and TP53 phenotypes in colon and rectal cancer . Therefore, it is possible that cross-talk between the TGFβ-signaling pathway and other pro-inflammatory pathways plays an important role in the determination of tumor molecular phenotype.
Our data also showed that eight genes in the TGFβ-signaling pathway are associated with differential miRNA expression. MiRNAs have previously been recognized as important in many pathways that regulate normal versus pathological function, and are dysregulated in most, if not all, cancers . In particular, miRNAs have been suggested as potential modulators of the TGFβ-signaling pathways’ seemingly paradoxical effects in regular colonic tissue and CRC development [24, 63]. Additionally, the canonical BMP pathways have been shown to alter gene expression levels through miRNAs . We have previously shown that miRNAs are extensively dysregulated in CRC  and that miRNAs can be used to differentiate between normal tissue and colonic carcinoma and rectal carcinoma at a molecular level .
In this study, we show that miRNAs play a role in the TGFβ’s function in CRC development, supporting earlier findings. Most earlier studies have focused on one, or several miRNAs such as miR-21 and miR-155 , miR-590 , miR-17-92  and miR-494 and miR-126-5p  or else focus on select genes in the TGFβ pathway, such as THBS1 . Here we assessed the associations between every gene in the KEGG TGFβ-signaling pathway and the expression levels of over 800 miRNAs. This allowed us to confirm earlier findings of associations with the miR-17-92 cluster (miR-17, 19, 20 and 92), in addition to detecting previously un-reported associations between a diverse array of miRNAs and genes in the TGFβ-signaling pathway. Moreover, we had previously reported that miR-17-5p, miR-20a-5p, miR-145-5p, miR-3651 and miR-92-3p were associated with both colon and rectal cancer, while miR-4749-3p was associated with rectal cancer only and miR-663a was associated with enriched biologic processes in normal colonic mucosa . Here we expand upon our previous findings and suggest that one of the mechanisms by which these miRNAs are associated with colorectal cancer may involve the dysregulation of the TGFβ pathway.
Our results suggest that the component of the TGFβ-signaling pathway most influenced by miRNAs is that part of the pathway that leads to apoptosis and cell cycle control (see Figs. 1, 2c). Dysregulated genes associated with miRNAs were TGFβR1, TGIF2, E2F4/E2F5, TFDP1, RBL1, and MYC. Within these associations it is clear that many genes are associated with several miRNAs and that the same miRNA could influence the expression, either directly or indirectly, of several genes. Notably in this pathway, miR-17-5p, miR-199a-5p, miR-20a-5p, miR-21-5p, and miR-93-5p have seed region matches with multiple genes as well as associations with other genes without a seed-region match, suggesting more of an indirect association. While it is beyond the scope of this paper to further assess these associations in laboratory-based functionality studies, these findings provide information that can be targeted in such studies to validate these findings. Upon validation, further studies to explore potential for targeting therapeutics based on these results can be undertaken.
Variations in clinical features of tumors have been noted by tumor subtype. Guinney and colleagues evaluated four consensus molecular subtypes (CMS), with CMS4 being characterized by activation of TGFβ-signaling . They observed that CMS4 tumors tended to be diagnosed at more advanced disease stages (stage III and IV) and had slightly worse survival. We noted variation within the TGFβ-signaling pathway in regard to both disease stage and survival months. However, the strongest differences by disease stage and survival months were observed when examining the interaction of miRNAs with mRNAs rather than mRNA expression. We believe that these findings may be important in identifying specific markers for further research and as possible treatment or screening modalities
A challenge in studies such as this, is determining which associations are important for follow-up and hold the potential for incorporating observations into clinical practice in the future. This is illustrated in our data where at the population level several genes were upregulated or downregulated, however a large percentage of cases did not demonstrate these characteristics. Two aspects of this observation are worth noting; first, mRNAs that had the greatest fold changes tended to have a greater percentage of the population expressing those changes. This also has implications for determining cutpoints when conducting statistical analysis. We focused our assessment of miRNAs with those genes that were statistically significant and also had a FC of > 1.50 or < 0.67, hoping to identify more biologically important genes. However, our data suggest that at that level fewer individuals expressed changes but those that did had greater changes. Secondly, we also observed less variability in the population for miRNA expression than for mRNA expression, in that miRNAs that were dysregulated tended to be dysregulated in the majority of the population, which was not the case for mRNAs.
This study has several strengths as well as some limitations. First, while our sample size is small, it is one of the largest available samples with paired tumor/normal data. When evaluating miRNA expression with mRNA expression, we could miss important associations since miRNAs have their impact post-transcriptionally. However, much of the current information on miRNA target genes comes from gene expression data, and any association observed may have important biological meaning, but must be acknowledged as being potentially incomplete [67, 68]. Our epidemiological approach should be seen as the first step in identifying associations that can be further examined in targeted laboratory studies to validate the identified associations that can then lead to targeted therapeutics. As such, this work should in part be viewed as a discovery of miRNAs that specifically are associated with genes in the TGFβ-signaling pathway that can be used to spearhead additional research in this field.
These data provide support for the role of miRNAs in CRC in part through their association with gene regulation in the TGFβ-signaling pathway. Through our comprehensive evaluation of this signaling pathway we have been able to identify genes within the pathway that are most dysregulated in CRC and how dysregulated miRNAs directly and indirectly influence those genes. While validation of study results is needed, these findings provide a basis for further laboratory evaluation and subsequent targeting of key genes and miRNAs.
false discovery rate
Kaiser Permanente of Northern California
Kyoto encyclopedia of genes and genomics
surveillance epidemiology and end results
Slattery ML, Wolff RK, Lundgreen A. A pathway approach to evaluating the association between the CHIEF pathway and risk of colorectal cancer. Carcinogenesis. 2015;36(1):49–59.
Gordon KJ, Blobe GC. Role of transforming growth factor-beta superfamily signaling pathways in human disease. Biochim Biophys Acta. 2008;1782(4):197–228.
Hong S, Lee C, Kim SJ. Smad7 sensitizes tumor necrosis factor induced apoptosis through the inhibition of antiapoptotic gene expression by suppressing activation of the nuclear factor-kappaB pathway. Cancer Res. 2007;67(19):9577–83.
Cancer Genome Atlas N. Comprehensive molecular characterization of human colon and rectal cancer. Nature. 2012;487(7407):330–7.
Heldin CH, Moustakas A. Signaling receptors for TGF-beta family members. Cold Spring Harb Perspect Biol. 2016;8(8):22053.
Yang G, Yang X. Smad4-mediated TGF-beta signaling in tumorigenesis. Int J Biol Sci. 2010;6(1):1–8.
Piccirillo SG, Vescovi AL. Bone morphogenetic proteins regulate tumorigenicity in human glioblastoma stem cells. Ernst Schering Found Symp Proc. 2006;5:59–81.
Piccirillo SG, Reynolds BA, Zanetti N, Lamorte G, Binda E, Broggi G, et al. Bone morphogenetic proteins inhibit the tumorigenic potential of human brain tumour-initiating cells. Nature. 2006;444(7120):761–5.
Alarmo EL, Kallioniemi A. Bone morphogenetic proteins in breast cancer: dual role in tumourigenesis? Endocr Relat Cancer. 2010;17(2):R123–39.
Kaneda H, Arao T, Matsumoto K, De Velasco MA, Tamura D, Aomatsu K, et al. Activin A inhibits vascular endothelial cell growth and suppresses tumour angiogenesis in gastric cancer. Br J Cancer. 2011;105(8):1210–7.
Nieminen TT, Abdel-Rahman WM, Ristimaki A, Lappalainen M, Lahermo P, Mecklin JP, et al. BMPR1A mutations in hereditary nonpolyposis colorectal cancer without mismatch repair deficiency. Gastroenterology. 2011;141(1):e23–6.
Rustgi AK. The genetics of hereditary colon cancer. Genes Dev. 2007;21(20):2525–38.
Hardwick JC, Kodach LL, Offerhaus GJ, van den Brink GR. Bone morphogenetic protein signalling in colorectal cancer. Nat Rev Cancer. 2008;8(10):806–12.
Sekiya T, Oda T, Matsuura K, Akiyama T. Transcriptional regulation of the TGF-beta pseudoreceptor BAMBI by TGF-beta signaling. Biochem Biophys Res Commun. 2004;320(3):680–4.
Teraoku H, Morine Y, Ikemoto T, Saito Y, Yamada S, Yoshikawa M, et al. Role of thrombospondin-1 expression in colorectal liver metastasis and its molecular mechanism. J Hepatobiliary Pancreat Sci. 2016;23(9):565–73.
Miyata N, Azuma T, Hozawa S, Higuchi H, Yokoyama A, Kabashima A, et al. Transforming growth factor beta and Ras/MEK/ERK signaling regulate the expression level of a novel tumor suppressor Lefty. Pancreas. 2012;41(5):745–52.
Pentek J, Parker L, Wu A, Arora K. Follistatin preferentially antagonizes activin rather than BMP signaling in Drosophila. Genesis. 2009;47(4):261–73.
Ambros V. The functions of animal microRNAs. Nature. 2004;431(7006):350–5.
Murray BS, Choe SE, Woods M, Ryan TE, Liu W. An in silico analysis of microRNAs: mining the miRNAome. Mol BioSyst. 2010;6(10):1853–62.
Arora S, Rana R, Chhabra A, Jaiswal A, Rani V. miRNA-transcription factor interactions: a combinatorial regulation of gene expression. Mol Genet Genom. 2013;288(3–4):77–87.
Gartel AL, Kandel ES. miRNAs: little known mediators of oncogenesis. Semin Cancer Biol. 2008;18(2):103–10.
Nam S, Li M, Choi K, Balch C, Kim S, Nephew KP. MicroRNA and mRNA integrated analysis (MMIA): a web tool for examining biological functions of microRNA expression. Nucleic Acids Res. 2009;37:W356–62.
Drusco A, Nuovo GJ, Zanesi N, Di Leva G, Pichiorri F, Volinia S, et al. MicroRNA profiles discriminate among colon cancer metastasis. PLoS ONE. 2014;9(6):e96670.
Welch DR, Hurst DR. Unraveling the ‘TGF-beta paradox’ one metastamir at a time. Breast Cancer Res. 2013;15(1):305.
Taylor MA, Sossey-Alaoui K, Thompson CL, Danielpour D, Schiemann WP. TGF-beta upregulates miR-181a expression to promote breast cancer metastasis. J Clin Invest. 2013;123(1):150–63.
Esser JS, Saretzki E, Pankratz F, Engert B, Grundmann S, Bode C, et al. Bone morphogenetic protein 4 regulates microRNAs miR-494 and miR-126-5p in control of endothelial cell function in angiogenesis. Thromb Haemost. 2017;117(4):734–49.
Chen L, Xu J, Chu X, Ju C. Micro RNA-98 interferes with Thrombospondin 1 expression in peripheral B cells of patients with asthma. Biosci Rep. 2017;37(4):BSR20170149.
Ekhteraei-Tousi S, Mohammad-Soltani B, Sadeghizadeh M, Mowla SJ, Parsi S, Soleimani M. Inhibitory effect of hsa-miR-590-5p on cardiosphere-derived stem cells differentiation through downregulation of TGFB signaling. J Cell Biochem. 2015;116(1):179–91.
Li L, Shi JY, Zhu GQ, Shi B. MiR-17-92 cluster regulates cell proliferation and collagen synthesis by targeting TGFB pathway in mouse palatal mesenchymal cells. J Cell Biochem. 2012;113(4):1235–44.
Yang H, Fang F, Chang R, Yang L. MicroRNA-140-5p suppresses tumor growth and metastasis by targeting transforming growth factor beta receptor 1 and fibroblast growth factor 9 in hepatocellular carcinoma. Hepatology. 2013;58(1):205–17.
Xu P, Wang J, Sun B, Xiao Z. Integrated analysis of miRNA and mRNA expression data identifies multiple miRNAs regulatory networks for the tumorigenesis of colorectal cancer. Gene. 2018;659:44–51.
Slattery ML, Trivellas A, Pellatt AJ, Mullany LE, Stevens JR, Wolff RK, et al. Genetic variants in the TGFbeta-signaling pathway influence expression of miRNAs in colon and rectal normal mucosa and tumor tissue. Oncotarget. 2017;8(10):16765.
Wu K, Zhao Z, Ma J, Chen J, Peng J, Yang S, et al. Deregulation of miR-193b affects the growth of colon cancer cells via transforming growth factor-beta and regulation of the SMAD3 pathway. Oncol Lett. 2017;13(4):2557–62.
Sells E, Pandey R, Chen H, Skovan BA, Cui H, Ignatenko NA. Specific microRNA-mRNA regulatory network of colon cancer invasion mediated by tissue kallikrein-related peptidase 6. Neoplasia. 2017;19(5):396–411.
Sun C, Wang FJ, Zhang HG, Xu XZ, Jia RC, Yao L, et al. miR-34a mediates oxaliplatin resistance of colorectal cancer cells by inhibiting macroautophagy via transforming growth factor-beta/Smad4 pathway. World J Gastroenterol. 2017;23(10):1816–27.
Slattery ML, Herrick JS, Pellatt DF, Stevens JR, Mullany LE, Wolff E, et al. MicroRNA profiles in colorectal carcinomas, adenomas and normal colonic mucosa: variations in miRNA expression and disease progression. Carcinogenesis. 2016;37(3):245–61.
Slattery ML, Potter J, Caan B, Edwards S, Coates A, Ma KN, et al. Energy balance and colon cancer–beyond physical activity. Cancer Res. 1997;57(1):75–80.
Slattery ML, Caan BJ, Benson J, Murtaugh M. Energy balance and rectal cancer: an evaluation of energy intake, energy expenditure, and body mass index. Nutr Cancer. 2003;46(2):166–71.
Slattery ML, Herrick JS, Mullany LE, Valeri N, Stevens J, Caan BJ, et al. An evaluation and replication of miRNAs with disease stage and colorectal cancer-specific mortality. Int J Cancer. 2015;137(2):428–38.
Slattery ML, Pellatt DF, Mullany LE, Wolff RK, Herrick JS. Gene expression in colon cancer: a focus on tumor site and molecular phenotype. Genes Chromosom Cancer. 2015;54(9):527–41.
Pellatt DF, Stevens JR, Wolff RK, Mullany LE, Herrick JS, Samowitz W, et al. Expression profiles of miRNA subsets distinguish human colorectal carcinoma and normal colonic mucosa. Clin Transl Gastroenterol. 2016;7:e152.
Agilent Technologies I. Agilent GeneSpring user manual. Santa Clara: Aglient Technologies Inc; 2013.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc. 1995;57:289–300.
Mullany LE, Herrick JS, Wolff RK, Slattery ML. MicroRNA seed region length impact on target messenger RNA expression and survival in colorectal cancer. PLoS ONE. 2016;11(4):e0154177.
Karolchik D, Hinrichs AS, Furey TS, Roskin KM, Sugnet CW, Haussler D, et al. The UCSC Table Browser data retrieval tool. Nucleic Acids Res. 2004;32:D493–6.
Slattery ML, Lundgreen A, Wolff RK, Herrick JS, Caan BJ. Genetic variation in the transforming growth factor-beta-signaling pathway, lifestyle factors, and risk of colon or rectal cancer. Dis Colon Rectum. 2012;55(5):532–40.
Singh A, Morris RJ. The Yin and Yang of bone morphogenetic proteins in cancer. Cytokine Growth Factor Rev. 2010;21(4):299–313.
Zhang L, Ye Y, Long X, Xiao P, Ren X, Yu J. BMP signaling and its paradoxical effects in tumorigenesis and dissemination. Oncotarget. 2016;7(47):78206–18.
Palles C, Chegwidden L, Li X, Findlay JM, Farnham G, Castro Giner F, et al. Polymorphisms near TBX5 and GDF7 are associated with increased risk for Barrett’s esophagus. Gastroenterology. 2015;148(2):367–78.
Becker J, May A, Gerges C, Anders M, Schmidt C, Veits L, et al. The Barrett-associated variants at GDF7 and TBX5 also increase esophageal adenocarcinoma risk. Cancer Med. 2016;5(5):888–91.
Grady WM, Pritchard CC. Molecular alterations and biomarkers in colorectal cancer. Toxicol Pathol. 2014;42(1):124–39.
Allaire JM, Roy SA, Ouellet C, Lemieux E, Jones C, Paquet M, et al. Bmp signaling in colonic mesenchyme regulates stromal microenvironment and protects from polyposis initiation. Int J Cancer. 2016;138(11):2700–12.
Whissell G, Montagni E, Martinelli P, Hernando-Momblona X, Sevillano M, Jung P, et al. The transcription factor GATA6 enables self-renewal of colon adenoma stem cells by repressing BMP gene expression. Nat Cell Biol. 2014;16(7):695–707.
Scarfi S. Use of bone morphogenetic proteins in mesenchymal stem cell stimulation of cartilage and bone repair. World J Stem Cells. 2016;8(1):1–12.
Chi Y, Yao L, Hu X, Huang S, Huang N, Li S, et al. The BMP inhibitor DAND5 in serum predicts poor survival in breast cancer. Oncotarget. 2016;7(12):14951–62.
Greco SA, Chia J, Inglis KJ, Cozzi SJ, Ramsnes I, Buttenshaw RL, et al. Thrombospondin-4 is a putative tumour-suppressor gene in colorectal cancer that exhibits age-related methylation. BMC Cancer. 2010;10:494.
Huang J, Soffer SZ, Kim ES, Yokoi A, Moore JT, McCrudden KW, et al. p53 accumulation in favorable-histology Wilms tumor is associated with angiogenesis and clinically aggressive disease. J Pediatr Surg. 2002;37(3):523–7.
Beck TN, Korobeynikov VA, Kudinov AE, Georgopoulos R, Solanki NR, Andrews-Hoke M, et al. Anti-mullerian hormone signaling regulates epithelial plasticity and chemoresistance in lung cancer. Cell Rep. 2016;16(3):657–71.
Kim BS, Kang HJ, Park JY, Lee J. Fucoidan promotes osteoblast differentiation via JNK- and ERK-dependent BMP2-Smad 1/5/8 signaling in human mesenchymal stem cells. Exp Mol Med. 2015;47:e128.
Sanchez-Capelo A. Dual role for TGF-beta1 in apoptosis. Cytokine Growth Factor Rev. 2005;16(1):15–34.
Slattery ML, Lundgreen A, Wolff RK. MAP kinase genes and colon and rectal cancer. Carcinogenesis. 2012;33(12):2398–408.
Svoronos AA, Engelman DM, Slack FJ. OncomiR or tumor suppressor? The duplicity of MicroRNAs in cancer. Cancer Res. 2016;76(13):3666–70.
Blahna MT, Hata A. Smad-mediated regulation of microRNA biosynthesis. FEBS Lett. 2012;586(14):1906–12.
Slattery ML, Herrick JS, Pellatt DF, Mullany LE, Stevens JR, Wolff E, et al. Site-specific associations between miRNA expression and survival in colorectal cancer cases. Oncotarget. 2016;7(37):60193.
Hata A, Davis BN. Control of microRNA biogenesis by TGFbeta signaling pathway-A novel role of Smads in the nucleus. Cytokine Growth Factor Rev. 2009;20(5–6):517–21.
Guinney J, Dienstmann R, Wang X, de Reynies A, Schlicker A, Soneson C, et al. The consensus molecular subtypes of colorectal cancer. Nat Med. 2015;21(11):1350–6.
Chou CH, Chang NW, Shrestha S, Hsu SD, Lin YL, Lee WH, et al. miRTarBase 2016: updates to the experimentally validated miRNA-target interactions database. Nucleic Acids Res. 2016;44(D1):D239–47.
Slattery ML, Herrick JS, Stevens JR, Wolff RK, Mullany LE. An assessment of database-validated microRNA target genes in normal colonic mucosa: implications for pathway analysis. Cancer Inform. 2017;16:1176935117716405.
AJP interpreted results and wrote the manuscript; MLS oversaw all aspects of study, obtained funding, helped write paper; LEM conducted bioinformatics analysis that included seed-region matches between mRNA and miRNA and helped write the manuscript; RKW oversaw laboratory genomic analysis; LCS contributed data and assisted in editing and writing the manuscript; WS provided pathology input; JSH managed data and conducted statistical analysis. All authors read and approved the final manuscript.
The contents of this manuscript are solely the responsibility of the authors and do not necessarily represent the official view of the National Cancer Institute. We acknowledge Sandra Edwards for data oversight and study management, and Michael Hoffman and Erica Wolff for miRNA analysis. We acknowledge Dr. Bette Caan and the staff at Kaiser Permanente Northern California for sample and data collection.
The authors declare that they have no competing interests.
Availability of data and materials
Data can be released in conjunction with study participant signed consent form. The datasets analyzed during the current study are not publicly available given the restrictions arising from the signed consent. To discuss feasibility of obtaining data contact Dr. Slattery.
Ethics approval and consent to participate
Study participants signed informed consent prior to study participation. The Institutional Review Boards at the University of Utah and at KPNC approved this study.
Consent for publication
This manuscript does not contain individual data presented in a manner that could identify them. All data are presented in aggregate.
This study was supported by NCI Grant CA163683.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.