Genome-Wide Uncovering of STAT3-Mediated miRNA Expression Profiles in Colorectal Cancer Cell Lines

Colorectal cancer (CRC) is one of the most common malignancies resulting in high mortality worldwide. Signal transducer and activator of transcription 3 (STAT3) is an oncogenic transcription factor which is frequently activated and aberrantly expressed in CRC. MicroRNAs (miRNAs) are a class of small noncoding RNAs which play important roles in many cancers. However, little is known about the global miRNA profiles mediated by STAT3 in CRC cells. In the present study, we applied RNA interference to inhibit STAT3 expression and profiled the miRNA expression levels regulated by STAT3 in CRC cell lines with deep sequencing. We found that 26 and 21 known miRNAs were significantly overexpressed and downexpressed, respectively, in the STAT3-knockdown CRC cell line SW480 (SW480/STAT3-siRNA) compared to SW480 transfected with scrambled siRNAs (SW480/siRNA-control). The miRNA expression profiling was then validated by quantitative real-time PCR for selected known miRNAs. We further predicted the putative target genes for the dysregulated miRNAs and carried out functional annotation including GO enrichment and KEGG pathway analysis for selected miRNA targets. This study directly depicts STAT3-mediated miRNA profiles in CRC cells, which provides a possible way to discover biomarkers for CRC therapy.


Introduction
MiRNAs are a class of endogenously expressed noncoding small RNAs. Over the past few years, miRNAs have been identified as abundant regulators of gene expression at posttranscriptional level across various biological processes, such as cancers and autoimmunity [1][2][3][4]. As is similar to proteincoding genes, the transcription of miRNAs is also regulated by an important class of gene regulators-transcription factors (TFs) which act at the transcriptional level. MiRNAs and TFs can cooperate to tune gene expression and form feedback or feed forward loops [5]. Therefore, the integration of TFs, miRNAs, and their target genes is required for transcriptional and posttranscriptional regulatory networks. However, the transcriptional networks in which TFs control miRNA expression (TFs-miRNA) have received little attention [6].
Signal transducers and activators of transcription (STAT), a family of latent cytoplasmic TFs which become activated by phosphorylation on a single tyrosine, can convey extracellular signals from kinds of cytokines, growth factors, or peptides to the nucleus [7]. They are critical signaling components which can drive tumorigenesis, cell proliferation, cell survival, and angiogenesis [8,9]. Among the STAT family members, STAT3 frequently responds to a variety of signals such as growth factors, cytokines, and oncogenes and is involved in diverse signaling pathways [10]. Numerous reports have taken STAT3 as a critical link between tumor cells and their microenvironments by regulating tumor growth and tumor-associated inflammation [7,11]. Constitutively activated or overexpressed STAT3 can be detected in multiple tumor-derived cell lines as well as samples from common malignant tumors, such as colon, 2 BioMed Research International skin, gastric, breast, and lung [12][13][14][15]. Conversely, conditional deletion of STAT3 results in reduced tumor development in mice [16]. In vitro and in situ studies show a frequent upregulation and activation of STAT3 protein in colorectal cancer (CRC) which is one of the most common carcinomas resulting in high mortality in western countries [17][18][19][20][21].
There is a body of evidence demonstrating the strong reciprocal regulations between STATs and miRNAs. STAT1/2 is upregulated under the transcriptional control of INF-alpha signaling after repression of miR-221/222 in glioblastoma U251 cells [22]. Interferon-can induce STAT1-dependent upregulation of the tumor suppressing miR-29 family in melanoma cells [23]. Inflammatory cytokines can increase miR-155 expression in human retinal pigment epithelial cells by activating STAT1 and enhancing putative STAT1 protein binding to the promoter region of miR-155 [24]. The most extensive studied STAT3/miR interaction is the STAT3/miR-21 pathway [25][26][27][28]. It is well accepted that STAT3 can directly activate miR-21, one of the miRNAs that promote cancer cell survival and proliferation [26,28]. MiR-21 is upregulated in many types of malignant tumors and has been identified as an antiapoptotic miRNA which can directly target programmed cell death protein 4 (PDCD4) and phosphatase and tensin homolog (PTEN) [29][30][31]. As a downstream effector of IL-6, STAT3 can activate transcription of miR-21 and miR-181b-1 which directly target PTEN and cylindromatosis (CYLD) tumor suppressor genes linking inflammation to cancer [25]. In contrast, miR-124 is reported as a CRC suppressor. It can program tumor cell apoptosis and suppress growth of the tumor by targeting and reducing STAT3 action [32].
Although a large number of reports have elaborated the STAT3-miRNAs relationship [10,[25][26][27][28][29][30][31], there is still a lack of thoroughly detailed study on the temporal dynamic regulation of miRNA expression library due to the action of STAT3 in CRC. In this paper, we used direct sequencing to document STAT3-mediated miRNA expression profiles in CRC cell line SW480. We identified a panel of differentially expressed known and novel miRNAs, which contribute to better understanding of STAT3-miRNAs interaction and involved signaling pathways in CRC cells.

Suppression of STAT3 by siRNA in SW480 Cells.
Knockdown of STAT3 in SW480 cells has been described in our previous study [32]. In brief, SW480 cells were transfected with 50 nM STAT3-siRNA or siRNA-control using siPORTNeoFX (Ambion). After 48 h of transfection, the protein levels of STAT3 were measured by western blotting.

RNA Extraction.
Total RNA was extracted from SW480 cells transfected with STAT3-siRNA or siRNA-control using TRIZOL reagent (Invitrogen) in accordance with the manufacture's protocol. RNA samples then passed the RNA quality control for deep sequencing.

Analysis of Sequencing Data.
High throughput sequencing was performed on the Illumina Cluster Station and Genome Analyzer II (Illumina Inc., USA). The raw sequences went through data cleaning process, including getting rid of the low quality sequences and several kinds of contamination formed by adaptor-adaptor ligation. Length distribution of clean sequences was then summarized to reveal the compositions of small RNA samples. The small RNA sequences of 18 to 30 nt were retained for further analyses. The standard bioinformatics analysis to annotate the clean sequences has been previously described [33]: (1) map all the small RNA sequences that pass filters to the reference human genome by Short Oligonucleotide Alignment Program (SOAP 2.0) and analyze their expression and distribution on the genome; (2) annotate the small RNA sequences with rRNA, small cytosol (sc)RNA, small nucleolar (sno)RNA, small nuclear (sn)RNA, and tRNA using Rfam 10.1 (http://rfam.sanger.ac.uk/) and Genbank (http://www.ncbi.nlm.nih.gov/genbank/) databases to get rid of matched sequences from unannotated sequences; (3) align the small RNA sequences to the miRNA precursor or mature miRNA of human species in miRBase18 (http://www.mirbase.org/) to get (a) the known miRNA count, (b) base bias on the first position among identified miRNAs with fixed length (18-30 nt), and (c) base bias on each position among all identified miRNAs; (4) align the small RNA sequences to repeated associated RNA to find matched sequences in the samples; (5) align the small RNA sequences to exons and introns of mRNA and match the small RNA to their original sites in genome; (6) if the clean sequences could not be annotated to match any category, we took them to predict the novel miRNAs. To make every specific small RNA mapped to only one annotation, we obeyed the following priority rule: rRNAetc (in which Genbank > Rfam) > known miRNA > repeat > exon > intron. The total rRNA ratio of less than 40% was a mark for sample quality check.
We used prediction software Mireap (http://sourceforge .net/projects/mireap/) to predict novel candidate miRNAs by detecting the secondary hairpin structure, the Dicer cleavage site, and the minimum free energy of the unannotated small RNA sequences which could be mapped to genome. The following criteria were conducted for defining highconfidence miRNA candidates: (1) the characteristic stable hairpin structure with low free energy (<-20 kcal/mol); (2) miRNA candidates expressed in both two samples at detectable levels (1 TPM, one transcript per million tags).
We compared the miRNA expression levels between two samples to detect the differentially expressed miRNAs.
The expression levels of miRNAs in two samples were first normalized to get the expression of TPM. The fold change and values were then calculated from the normalized expression level. In general, if the adjusted values were <0.01 based on the Benjamini and Hochberg multiple testing correction and there was at least a 2-fold change ((SW480/STAT3-siRNA)/(SW480/siRNA-control)) in the normalized expression, one could consider the miRNAs as significantly differentially expressed [34].
We applied RNAHybrid (http://bibiserv.techfak.unibielefeld.de/rnahybrid/) to predict potential targets of differentially expressed miRNAs by detecting the minimum free energy hybridization of the small RNA sequences and mRNAs. The parameters were set like this: -c -d 1.9, 0.28 -t cel-hbl-1.fasta -q cel-let-7.fasta. The target genes of dysregulated miRNAs were annotated by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses to predict their potential functions.

Validation of miRNA Expression by Quantitative RT-PCR.
Assays to quantify the known and novel miRNAs were done by using miScript PCR System (Qiagen) according to the manufacturer's instruction. RT reactions with miScript II RT Kit (Qiagen) contained 1.0 g total RNA, 4 L 5 × miScriptHiSpec buffer, 2 L 10 × miScriptnucleics mix, and 2 L miScript reverse transcriptase mix in each reaction (20 L). The RT reaction was conducted under the following conditions: 37 ∘ C for 60 min and then 95 ∘ C for 5 min. After that, the cDNA products from RT reaction were diluted 15 times. PCR was carried out with 1.5 L of the diluted products in 20 L PCR reaction containing 10 L 2 × QuantiTect SYBR Green PCR master mix, 2 L 10 × miScript universal primer, and 2 L 10 × miScript primer assay. Amplification was performed as follows: 95 ∘ C for 15 min, followed by 40 cycles at 94 ∘ C for 15 s, 55 ∘ C for 30 s, and 70 ∘ C for 30 s. All reactions were run in triplicate. Relative expression was calculated using the comparative CT method and normalized to the expression of RNU6B.

Sequencing Data Description.
Suppression of STAT3 by STAT3-siRNA in SW480 cells was described by our previous study [32]. We obtained ∼14M and ∼10M high-quality clean reads from the raw sequences after removing contaminants in SW480/siRNA-control and SW480/STAT3-siRNA samples, respectively (Table S1 in the Supplementary Material available online at http://dx.doi.org/10.1155/2014/187105). The more clean reads in SW480/STAT3-siRNA samples might result from the small interfering RNA. We then summarized the length distribution of these clean reads, which is helpful to discover the compositions of small RNA samples. The most abundant group in both samples was 22 nt in length which is consistent with numerous studies of miRNA size distribution reported in human and animals ( Figure S1). In both SW480/siRNA-control and SW480/STAT3-siRNA samples, most 22-nt small RNAs began with the base "U. " The 23-nt small RNAs displayed a bias to "A" and "U" at first base in SW480/STAT3-siRNA while only "U" in SW480/siRNAcontrol. Both libraries exhibited similar compositions of four bases and most sites kept a major base, indicating that small RNA sequences from libraries were relatively conserved ( Figure S2).
A total of 829,863 unique tags were obtained after removing repeats from 24,669,643 total tags in both SW480/siRNAcontrol and SW480/STAT3-siRNA samples. The two samples shared 19.17% of unique common tags and 96.42% of total common tags. The very few unique common tags indicated that SW480/STAT3-siRNA presented a distinctive small RNA profile compared to SW480/siRNA-control. The high percentage of total common tags might suggest that miRNAs have been successively enriched from both libraries. SW480/STAT3-siRNA had more specific small RNAs than SW480/siRNA-control which might be also due to the small interfering RNA experiments (unique: 51.39% and 29.43%; total: 2.35% and 1.23% for SW480/STAT3-siRNA and SW480/siRNA-control, resp.) ( Figure S3).
The unique tags and total tags were both mapped onto human genomes, with a total of 29.15% (117,573 tags) in SW480/siRNA-control and 38.56% (225,833 tags) in SW480/STAT3-siRNA for unique tags, while about 74.90% (7,930,277 tags) in SW480/siRNA-control and 76.66% (10,795,100 tags) in SW480/STAT3-siRNA for total tags. As shown in Figure S4, small RNAs were unevenly distributed across chromosomes between sense and antisense chains. More small RNAs were mapped in the sense chains. For example, there were 27,276 unique tags and 983 unique tags in exon-sense and exon-antisense regions in SW480/siRNA-control, respectively, with ratio of 28 : 1; while in SW480/STAT3-siRNA, the ratio was 68 : 1 (100,478 : 1,482) (Figures 1(a) and 1(c)). The numbers of total tags were 34,158 and 1,494 in exon-sense and exon-antisense chains in SW480/siRNA-control, while 136,505 and 3,528 in SW480/STAT3-siRNA (Figures 1(b) and 1(d)).
The remaining unannotated small RNAs were used for novel miRNA prediction, accounting for 68.05% and 57.87% of total reads in SW480/siRNA-control and SW480/STAT3-siRNA cells, respectively. We predicted the sequences with miRNA stem loop structure and dicer cleavage sites from unannotated sequences to be novel miRNAs. A total of 52,597 and 44,685 sequence tags were identified to be  58 and 70 novel miRNAs in SW480/siRNA-control and SW480/STAT3-siRNA, respectively, 30 of which were shared by both cells.
To investigate the differential expression levels of individual miRNA between SW480/siRNA-control and SW480/STAT3-siRNA, we compared their miRNA profiles and identified that 26 and 21 known miRNAs were overexpressed and downexpressed significantly (adjusted < 0.01, Table 1, Figure 2(a)) in SW480/STAT3-siRNA cells, respectively.

Targets and Pathway Prediction of Known miRNAs.
The dysregulated known miRNAs were predicted to target about 30,846 genes based on RNAHybrid prediction. The GO enrichment analysis for all the target genes showed that they enriched significantly in regulation of biological processes and gene expression, such as "cellular process, " "regulation of metabolic process, " "regulation of biological process, " and "regulation of gene expression" (see details in Table 3, adjusted < 0.05, Bonferroni correction). A single mRNA could be controlled posttranscriptionally by hundreds of miRNAs. As illustrated in Figure 3, one target can be affected by multiple miRNAs from both upand downregulated groups. Among the numerous targets, MICAL3, DDX58, and DOCK9 were inhibited potentially by 15 upregulated miRNAs and ANK3 was targeted by 13 downregulated miRNAs. MICAL3 is involved in zinc ion binding (GO: 0008270). DDX58 is ATP-dependent RNA helicase, associating with ATP binding (GO: 0005524), ATPdependent helicase activity (GO: 0008026), and nucleic acid binding (GO: 0003676). DOCK9 is a dedicator of cytokinesis protein, joining in small GTPase-mediated signal transduction (GO: 0007264) and guanyl-nucleotide exchange  factor activity (GO: 0005085). ANK3 takes part in biological process of axon guidance (GO: 0007411), protein targeting to plasma membrane (GO: 0019228), and synapse organization (GO: 0050808). Additionally, WDR35, DDX4, PDPK1, CLTC, WDR82, and DNAH17 were targeted by 14 overexpressed miRNAs, while TET2, SYT2, and SRRM4 were mediated by 12 downexpressed miRNAs.
To avoid the bias of speculative interference, only the target genes mediated by more than 3 up-or downexpressed miRNAs were selected to perform KEGG pathway analysis. Among the significantly enriched KEGG pathways, "Pathway in Cancer (hsa05200)" was the top one in both KEGG analyses for targets of up-or downregulated miRNAs (Table 4, adjusted < 0.05, Bonferroni correction). More than 200 target genes joined in this complex pathway. STAT3 was expected to participate in the signaling system by activating the JAK-STAT3 pathway. Obviously, the effects were predicted to extend from the JAK-STAT3 pathway to multiple networks in the cancer cells.
The targets of downregulated miRNAs were supposed to be overexpressed in STAT3-silence cells due to the decreased level of related miRNAs. We observed that, except "Pathway in Cancer, " there are still signaling related pathways enriched significantly in targets of downregulated miRNAs, such as focal adhesion (hsa04510), calcium signaling pathway (hsa04020), and colorectal cancer (hsa05210) and other pathways ( Table 4, adjusted < 0.05, Bonferroni correction). However, the targets of upregulated miRNAs were predicted to be inhibited by those increasing miRNAs in STAT3silence cells. These targets took part in signaling pathways such as Axon guidance (hsa04360), ErbB signaling pathway (hsa04012), focal adhesion (hsa04510), and insulin signaling pathway (hsa04910) ( Table 4, adjusted < 0.05, Bonferroni correction).

Targets and Pathway Prediction of Novel miRNAs.
The top five significant GO terms of novel miRNA targets were neurogenesis, cell development, transcription from RNA polymerase II promoter, regulation of transcription from RNA polymerase II promoter, and generation of neurons (Table S2, adjusted < 0.05, Bonferroni correction). KEGG analysis showed that the STAT3-mediated novel miRNAs influence the cancer-related signaling pathway significantly, including ErbB signaling pathway, MAPK signaling pathway, and Notch signaling pathway (Table S3, adjusted < 0.05, Bonferroni correction).

Discussion
Over the past years, the study of miRNA biogenesis and function has made remarkable progress. However, the upstream regulators of miRNAs as well as the mechanisms that miRNAs use to control gene expression are still mysterious. STAT3 is known as an oncogene and the consistent activation of STAT3 proteins has been characterized in tumor tissues or in cell lines derived from human tumors in accumulated reports [12][13][14][15]. At present we are still far from completely understanding the interaction of STAT3 and miRNAs. In this study, we identified the STAT3-mediated miRNA expression profiles through deep sequencing technique, which was followed by validation with qRT-PCR. We then conducted functional annotation for the predicted target genes of these differentially expressed miRNAs with GO and KEGG analyses.
In general, 26 known miRNAs were discovered to be overexpressed in SW480/STAT3-siRNA cells, of which miR-934 was the most upregulated. To our knowledge, this miRNA has seldom been studied before until Lui et al. reported its existence in 2007 [35]. The most abundant miRNA-miR-10a-5p, although not the top upregulated one in this study, has been demonstrated in a variety of cancers. For example, it is aberrantly overexpressed in Nucleophosmin 1-(NPM1-) mutated acute myeloid leukemia (AML) and interferes with the p53 pathway partly through MDM4 [36,37]. In CRC samples and corresponding normal tissues, miR-10a-5p has also been shown to be one of the most abundantly expressed miRNAs revealed by other deep sequencing studies [38,39]. In contrast, we found that 21 known miRNAs were significantly downexpressed in SW480/STAT3-siRNA cells. MiR-3656 was the most downregulated one in our study which had been identified to be significantly downregulated in colorectal carcinogenesis only in one literature [40]. MiR-25-5p was the most abundant miRNA among those significantly downregulated miRNAs in both samples. It was found that miR-25 was downregulated in human colon cancer tissues, contributing to promote cell proliferation and migration. Therefore, it might function as a tumor suppressor by targeting Smad7 in colon cancer [41]. MiR-215 was the second abundant miRNA in the significantly downregulated group, which was reported as a biomarker for colon cancer [42]. MiR-143 was frequently downregulated in colorectal carcinoma tissues, the restoration of which in colon cell lines might reduce tumor cell growth and soft-agar colony formation [43,44]. It was shown to target DNA methyltransferase 3A (DNMT3A) directly by luciferase reporter assay [43] and Hexokinase 2 (HK2) by microarray in combination with seed site enrichment analysis [44].
Recently, Rozovski et al. used microarray approach to show that STAT3 can directly and indirectly modulate miRNA expression levels in B-cell chronic lymphocytic leukemia (CLL) cells [45]. In their results, miR-21 and miR-181a were both downregulated by transfection with STAT3-shRNA. However, our results showed a different profile. MiR-181a was upregulated in SW480/STAT3-siRNA samples, which might be attributed to different cancer tissues. KEGG analysis of novel miRNA candidate target revealed that ErbB signaling pathway is the most significant. The ErbB family of receptor tyrosine kinases can be autophosphorylated on specific tyrosine residues for phosphotyrosine binding as well as cytoplasmic signaling molecules to activate numerous intracellular signaling pathways. The STAT proteins are one of the intermediates in the pathway cascade [46]. ErbB1 is able to phosphorylate STAT1 and STAT3 in vitro [47,48]. In A431 cells, STAT1, STAT3, and STAT5 were constitutively complexed with ErbB1 and rapidly phosphorylated on tyrosine in response to EGF [46].
Apparently, the development of miRNA microarrays, RT-PCR platforms, in situ hybridization, and next generation sequencing methodologies lead to a gradually growing number of miRNA profiling studies, which have paved the way to new approaches for biomarker discovery. People have reached a consensus that differential expression of miRNAs in cancers may have substantial diagnostic and prognostic values. Taken together, our data directly showed that STAT3 can modulate the miRNA expression levels in CRC cells, which contribute to better understanding of STAT3-miRNA interaction as well as the involved complex regulatory network.