LAMA4-Regulating miR-4274 and Its Host Gene SORCS2 Play a Role in IGFBP6-Dependent Effects on Phenotype of Basal-Like Breast Cancer

Specificity of RNAi to selected target is challenged by off-target effects, both canonical and non-canonical. Notably, more than half of all human microRNAs are co-expressed with hosting them proteincoding genes. Here we dissect regulatory subnetwork centered on IGFBP6 gene, which is associated with low proliferative state and high migratory activity of basal-like breast cancer. We inhibited expression of IGFBP6 gene in a model cell line for basal-like breast carcinoma MDA-MB-231, then traced secondary and tertiary effects of this knockdown to LAMA4, a laminin encoding gene that contributes to the phenotype of triple-negative breast cancer. LAMA4-regulating miRNA miR-4274 and its host gene SORCS2 were highlighted as intermediate regulators of the expression levels of LAMA4, which correlated in a basal-like breast carcinoma sample subset of TCGA to the levels of SORCS2 negatively. Overall, our study points that the secondary and tertiary layers of regulatory interactions are certainly underappreciated. As these types of molecular event may significantly contribute to the formation of the cell phenotypes after RNA interference based knockdowns, further studies of multilayered molecular networks affected by RNAi are warranted.

Specificity of RNAi to selected target is challenged by off-target effects, both canonical and non-canonical. Notably, more than half of all human microRNAs are co-expressed with hosting them proteincoding genes. Here we dissect regulatory subnetwork centered on IGFBP6 gene, which is associated with low proliferative state and high migratory activity of basal-like breast cancer. We inhibited expression of IGFBP6 gene in a model cell line for basal-like breast carcinoma MDA-MB-231, then traced secondary and tertiary effects of this knockdown to LAMA4, a laminin encoding gene that contributes to the phenotype of triple-negative breast cancer. LAMA4-regulating miRNA miR-4274 and its host gene SORCS2 were highlighted as intermediate regulators of the expression levels of LAMA4, which correlated in a basal-like breast carcinoma sample subset of TCGA to the levels of SORCS2 negatively. Overall, our study points that the secondary and tertiary layers of regulatory interactions are certainly underappreciated. As these types of molecular event may significantly contribute to the formation of the cell phenotypes after RNA interference based knockdowns, further studies of multilayered molecular networks affected by RNAi are warranted.

INTRODUCTION
MicroRNAs constitute one of the largest families of non-coding RNAs and are known to regulate the proliferative and migratory activity of tumor cells (Makarova et al., 2015). According to the current estimates, microRNAs regulate expression of more than 60% of protein-coding genes (Makarova et al., 2016). An increase in levels of certain microRNA may lead to a decrease in mRNA levels of their target genes and vice versa. In turn, this may lead to significant alterations of cell phenotypes (Galatenko et al., 2018a). RNA interference (RNAi) is a widely used tool for specific downregulation of expression of a certain target gene. Its central assumption, namely, its specificity to selected target, has been challenged by commonly detected off-target effects, driven both by "seed" sequence based annealing of shRNA to the genes encoding related proteins of same family (Jackson et al., 2006) or even to unrelated mRNA sequences, which may be targeted by either canonical or noncanonical interactions (Birmingham et al., 2006;Seok et al., 2018).
Notably, microRNA encoding genes may be arranged not only as independent transcription units but also within the introns and exons of both protein-coding and non-coding RNA genes (Makarova et al., 2016). More than half of all human microRNAs are co-expressed with the proteincoding genes and are located within their introns (Hinske et al., 2010(Hinske et al., , 2014. Previously, we have demonstrated the regulatory interactions between the protein-coding genes which host certain miRNAs and the protein-coding target genes of these miRNA (Galatenko et al., 2018a).
In this paper we concentrate on a dissection of a particular regulatory subnetwork centered on IGFBP6 gene, which is associated with low proliferative state and high migratory activity of basal-like breast cancer (Poloznikov et al., 2019). Previous studies suggest that RNAi-guided downregulation of IGFBP6 gene expression leads to an elevation in the expression of genes encoding miR-100 and miR-let-7a (Poloznikov et al., 2019). Surprisingly, known targets of the above-mentioned microRNAs include IGF2, IGF1R, and INSR genes, known as drivers of the process of metastasis (Nikulin et al., 2018). In other words, suppression of anti-migratory gene IGFBP6 may take place on two different levels, direct and indirect, with the latter being realized through miRNA based suppression of INSR/IGF1R 1 axis of pro-metastatic signaling (Poloznikov et al., 2019).
To explore putative miRNA-based co-regulation followed by the suppression of the IGFBP6 gene, we inhibited expression of this gene in MDA-MB-231 cells, a model cell line for basallike breast carcinoma, by shRNA, then we have traced secondary and tertiary effects of this knockdown to LAMA4, a laminin encoding gene that contributes to the phenotype of triplenegative breast cancer.

For
RNAi-mediated knockdown of the IGFBP6 gene, the MDA-MB-231 cells were transduced with lentiviral expression vector pLVX-shRNA1 (Clontech Laboratories, USA) carrying the following shRNAencoding sequence: 5 ′ -gatccGCCCAATTGTGACCATCGA TCGATTCAAGAGATCGATCGTCACAATGGGGCTTTTTTT ACGCGTg-3 ′ ). The control cell line was transduced with the same vector encoding shRNA targeting firefly (Photinus pyralis) luciferase gene. The production of lentiviral particles, transduction, and selection procedures were performed according to the protocol described earlier (Knyazev et al., 2015). Phenotypic characteristics of resultant MDA-MB-231 cell line with stable knockdown of IGFBP6 were described in our earlier work (Nikulin et al., 2018). The levels of endogenous expression of IGFBP6 in this stable cell line are at 22% of that in luciferase expression control.

RNA Extraction and Analysis of Expression
Before RNA isolation, cells were suspended in 2.5 ml of culture media, plated onto 6-wells plates (2 × 10 5 cells/ml) and incubated for 48 h at 37 • C in 5% CO 2 . Afterward, the cells were washed three times with cold FSB solution (PanEco) and lysed in Qiazol Reagent (Qiagen, Germany). RNA was isolated with miRNeasy Mini Kit (Qiagen) according to the manufacturer's protocol. The concentrations and the quality of the purified RNA were assessed using both Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific, USA) and Agilent Bioanalyzer 2100 (Agilent Genomics, USA) (Fomicheva et al., 2017). A260/A280 ratios were in range of 1.92.2. Each RNA sample was isolated from independently grown flask, with three flasks per cell type undergoing mRNA extraction. RNA quality controls were performed as described in Krainova et al. (2013). The RNA integrity number (RIN) values were higher than 9.0.
Procedures for the preparation of total RNA samples for hybridization on Human Transcriptome Arrays 2.0 microarrays (Affymetrix, USA) were carried out according to the manufacture protocol as described in Khaustova et al. (2017) and Kudriaeva et al. (2017). cDNA synthesis was performed using 500 ng total RNA as the starting material. Hybridization of fragmented and labeled cDNA on microarrays, array washing, staining, and scanning were performed as described in Sakharov et al. (2012). Experiment has been repeated in three biological replicates. Scans of the microarrays were converted into CEL files using the scanner software, and were then processed in Transcriptome Analysis Console 4.0.1.36 according to the summation method Gene + Exon -SST-RMA (Analysis Version 2) and the ebayes statistics; respective p-values were calculated after Benjamin-Hochberg correction for multiple tests. The fold changes in gene expression were shown using a log2 scale.
GSEA analyses were performed the Pathway Studio environment (www.pathwaystudio.com, accessed in Sep 2019).

Bioinformatics Analysis
The bioinformatics data analysis was carried out in R statistical environment according to the following algorithm: 1. Genes, which have changed their expression 2-fold or larger (with FDR p ≤ 0.01) after IGFBP6 knockdown were selected; 2. Among those genes, the genes encoding microRNAs were further determined (Hinske et al., 2014); 3. For each intragenic miRNA identified at Step 2, miRTarBase 7.0 (Chou et al., 2017) was utilized to compile the list of its target genes; 4. Target lists were analyzed to identify any target genes with expression changed of more than four-fold (at FDR p ≤ 0.01); 5. The directionality of the expression changes of each host gene to each target gene was explored in triple-negative breast cancer sample subset of TCGA database (N = 190), accessed on July 28th, 2019.

RESULTS
Transcriptional knockdown of IGFBP6 gene in MDA-MB-231 cells resulted in inhibition of cell migration with concomitant increase in the rate of proliferation (Nikulin et al., 2018). A comparison of IGFBP6 knockdown and the control population of cells revealed a significant change in expression levels of 380 genes (FDR p ≤ 0.01, fold change ≥ 2). Among those, 183 genes were downregulated, while the expression of other 197 genes was increased. Enrichment analysis revealed that "Mesenchymal Cell Differentiation" as the most altered Molecular Function, along with a predominance of "Cell Adhesion" and "Cell Communication" centered gene lists among the most altered biological processes (Supplementary Table 1). When analyses of upregulated and downregulated genes were performed separately, cellular and molecular functions performed by upregulated components were minimally perturbed when compared to that of entire list of genes with altered expression (Supplementary Table 2). Surprisingly, Top cellular component perturbed by downregulation was "vesicle, " while the Top of the list of downregulated biological processes was dominated by various aspects of tissue morphogenesis (Supplementary Table 3). Importantly, a total of 29 upregulated genes were co-located with various microRNAs ( Table 1). These co-located intragenic miRNA loci encode a total of 35 pre-miRNAs and 41 different mature microRNAs.
Mining of miRTarBase 7.0 database (Chou et al., 2017) for possible interactions of these 41 microRNAs resulted in a total of 3,429 miRNA-target gene interactions cataloged. For each of these target genes, gene expression changes were assessed by applying the threshold of FDR p ≤ 0.01 and at least 4fold change of both miRNA host and miRNA target genes. According to these criteria, there was only one host/target gene pair with significant decrease of miRNA host and concomitant increase of respective target gene expression, namely, SORCS2 as a host (a decrease of 4.32-fold) and LAMA4 as a target (an increase of 5.46 fold). The regulator connecting LAMA4 and other target genes was hsa-miR-4274 encoded by the second intron of SORCS2. Figure 1 shows a scheme of SORCS2 gene with miR-4274 location.
Taking into account that interactions between hsa-miR-4274 and LAMA4 gene has been detected in course of PAR-CLIP (Whisnant et al., 2013) experiment only, secondary analysis of the host and target genes expression was performed in the basallike breast cancer dataset of 190 samples retrieved from TCGA and annotated using the pam50 classifier (geneFu package).
The results of SORCS2 and LAMA4 genes co-expression analysis are shown in Figure 2. Notably, in basal-like breast carcinoma subdataset of TCGA, expression levels of these two genes negatively correlate with each other (Spearman's R = −0.57, p < 2.2 × 10 −16 ).

DISCUSSION
miRNAs are well-known modulators of gene expression, playing their roles both at the level of the transcripts and through epigenetic gene silencing. Although significant progress has been made in understanding functions of individual miRNA, mechanisms regulating their own transcription remain elusive. Speaking generally, intragenic miRNAs are expected to be coexpressed with their host genes, displaying either synergistic or antagonistic correlation patterns. Because of that, shRNA based suppression of the activity of individual gene, which hosts miRNA in its intron, is expected to change expression pattern for both the host gene and the miRNA it harbors.
IGFBP6 gene, a well-known suppressor of IGF-II signaling, suppresses proliferation and stimulates apoptosis in a wide range of cells. On the other hand, in some tumor cells, including those of rhabdomyosarcoma and basal types of breast cell carcinoma (Poloznikov et al., 2019), it induces migration, likely by an IGF-independent mechanism. Respectively, depending on particular cell type, effects of IGFBP6 may be classified as that of a tumor suppressor, or a facilitator of the metastatic spread. Thus, understanding the physiological role of IGFBP6 protein may be context-dependent and requires thorough evaluation of the relative balance of its IGF-dependent and IGF-independent actions.
In a previous study, shRNA based suppression of IGFBP6 in basal-like breast cancer cell line MDA-MB-231 led to a notable inhibition of cell migration along with an increase in its rate of proliferation (Poloznikov et al., 2019). GSEA analysis showed that knockdown of IGFBP6 profoundly changes expression programs related to cell adhesion and cell communication, while sparing machinery involved in propagation of the cell cycle. Among pathways most enriched by downregulated mRNAs were ones directly involved in tissue morphogenesis, and in the formation of vesicles, including exosomes. These effects were accompanied by the elevation in the expression of a metastasissuppressing miR-100 and miR-let-7a (Poloznikov et al., 2019). These two observations may be tied together by the hypothesis of intracellular accumulation of miRNAs due to suppression of their transport into intercellular space and/or circulation.
Further dissection of the transcription program enacted in response to the treatment with IGFBP6 shRNA revealed concerted changes in the levels of many transcripts lacking canonical complementarity to its intended target. Many previous works have previously acknowledged detrimental miRNA-like effects of shRNA (Anderson et al., 2008;Seok et al., 2018),  including unwanted toxicities, and their clear potential for misleading an interpretation of resultant phenotypes. In this paper, we show that off-target effects of a transcriptional knockdown of IGFBP6 are multilayered. After knockdown of IGFBP6, a total of 380 genes change their expression significantly, representing the first layer of the off-target complexity. As many of these genes (N = 29) hosted various intragenic pre-microRNA (Table 1), we attempted to uncover a second regulatory layer by embarking on a further analysis of the putative targets of these miRNA. Indeed, analysis of coordinated changes in expression of the host genes along with the targets for miRNAs embed in these genes revealed a number of candidate pairs. Expression levels of the strongest co-regulated gene pair, SORCS2 as a host and LAMA4 as a target, correlated negatively (R = −0.57, p < 2.2 × 10 −16 ), when analyzed in a basal-like breast carcinoma sample subset of TCGA (N = 190) (Figure 2). Interestingly, gene variant located in a vicinity of SORCS2 (rs4234798) demonstrated genome-wide significant associations with circulating concentrations of IGFBP-3 and IGF-I (Kaplan et al., 2011;Teumer et al., 2016), previously described as a prognostic marker in basal-like triple-negative breast cancer (Marzec et al., 2015;Julovi et al., 2018), and as a part of 19-gene expression signature which substantially outperforms Dukes' classification predicting the survival of patients with colorectal cancer (Aziz et al., 2016). Protein encoded by SORCS2 gene coordinates intracellular transport of the glutamate and cysteine transporter EAAT3 (Malik et al., 2019) in the neurons, where it provides a degree of protection from oxidative stress and epilepsy-induced cell death (Shkurnikov et al., 2013). Given that expression of SORCS2 mRNA is highly restricted to the nervous system and the kidneys, it is tempting to speculate that the function of relatively high levels of SORCS2 expression observed in basal-like breast carcinomas should be attributed to intragenic miRNA hsa-miR-4274.
As many other miRNAs, miR-4274 is capable of binding to many different mRNAs. One of its confirmed targets, the LAMA4 mRNA, encodes an α4-chain of laminins, the secreted heterotrimeric glycoproteins involved in the induction of the epithelial-mesenchymal transition. Elevation of LAMA4/LAMA5 protein ratio correlates with an increase in the permeability of the basal membrane as well as the metastatic and invasive potential of colorectal cancer (Galatenko et al., 2018b). Elevated level of LAMA4 expression as well as two other chains of laminin 411 is observed in metastatic lymph nodes of prostate cancer patients in comparison with metastasis-free lymph nodes (Shkurnikov et al., 2018). In breast carcinomas, an increase in LAMA4 expression marks the transition from pre-malignant lesions to malignant carcinomas, while in well-developed tumors an overexpression of LAMA4 mRNA is associated with shorter relapse-free survival (Ross et al., 2015). Specifically, in basal-like triple negative breast carcinoma cells, the levels of both LAMA4 mRNA and laminin-α4 are elevated when compared to that in adjacent normal tissue samples. Moreover, knockdown of LAMA4 inhibits proliferation, migration and invasion of these cancer cells in vitro (Maltseva and Rodin, 2018;Yang et al., 2018). Taken together, observations concerning LAMA4, its regulating miRNA miR-4274 and its host gene SORCS2, point at a substantial role played by secondary and tertiary effects produced by a transcriptional knockdown of IGFBP6 in basal-like breast carcinoma cells.
Overall, our study shows that the secondary and tertiary layers of regulatory interactions are certainly underappreciated. This study is limited in its scope, as it studies results of IGFBP6 knockdown in only one model cell line rather than in a panel of cells of various tissue origins. Our findings were validated only in a specific type of human tumor, basallike triple negative breast carcinoma and, therefore, cannot be generalized to all cancers. Moreover, preliminary data point that negative correlation of expression levels for genes SORCS2 and LAMA4 is a phenomenon which is, indeed, specific to only this tumor type. While it is tempting to speculate that strict cell type specificity of observed phenomenon may inform future therapy for triple negative breast cancer, detailed studies of SORC2/LAMA4 relationship are warranted. Similarly, further studies of multilayered molecular networks to be uncovered by deconvoluting of miRNA mediated correlations in gene expression patterns should be performed in order to improve our understanding of the formation of cell phenotypes after RNA interference based knockdowns.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://www.cancer.gov/tcga.

AUTHOR CONTRIBUTIONS
MS, AP, AB, and AT contributed conception and design of the study. MS and SNe performed the statistical analysis. SZ performed enrichment analysis. MS wrote the first draft of the manuscript. SNi, AB, US, DW, and AT wrote sections of the manuscript. All authors contributed to manuscript revision, read and approved the submitted version.