**Transcription Factor and miRNA Interplays Can Manifest the Survival of ccRCC Patients**

### **Shijie Qin, Xuejia Shi, Canbiao Wang, Ping Jin \* and Fei Ma \***

Laboratory for Comparative Genomics and Bioinformatics, College of Life Science, Nanjing Normal University, Nanjing 210046, China; 191201005@stu.njnu.edu.cn (S.Q.); 181202063@stu.njnu.edu.cn (X.S.); 191202061@stu.njnu.edu.cn (C.W.)

**\*** Correspondence: jinping@njnu.edu.cn (P.J.); mafei01@tsinghua.org.cn (F.M.)

Received: 28 September 2019; Accepted: 24 October 2019; Published: 28 October 2019

**Abstract:** Clear cell renal cell carcinoma (ccRCC) still remains a higher mortality rate in worldwide. Obtaining promising biomakers is very crucial for improving the diagnosis and prognosis of ccRCC patients. Herein, we firstly identified eight potentially prognostic miRNAs (hsa-miR-144-5p, hsa-miR-223-3p, hsa-miR-365b-3p, hsa-miR-3613-5p, hsa-miR-9-5p, hsa-miR-183-5p, hsa-miR-335-3p, hsa-miR-1269a). Secondly, we found that a signature containing these eight miRNAs showed obviously superior to a single miRNA in the prognostic effect and credibility for predicting the survival of ccRCC patients. Thirdly, we discovered that twenty-two transcription factors (TFs) interact with these eight miRNAs, and a signature combining nine TFs (*TFAP2A*, *KLF5*, *IRF1*, *RUNX1*, *RARA*, *GATA3*, *IKZF1*, *POU2F2*, and *FOXM1*) could promote the prognosis of ccRCC patients. Finally, we further identified eleven genes (hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-1269a, hsa-miR-144-5p, hsa-miR-183-5p, hsa-miR-335-3p, *TFAP2A*, *KLF5*, *IRF1*, *MYC*, *IKZF1*) that could combine as a signature to improve the prognosis effect of ccRCC patients, which distinctly outperformed the eight-miRNA signature and the nine-TF signature. Overall, we identified several new prognosis factors for ccRCC, and revealed a potential mechanism that TFs and miRNAs interplay cooperatively or oppositely regulate a certain number of tumor suppressors, driver genes, and oncogenes to facilitate the survival of ccRCC patients.

**Keywords:** ccRCC; prognostic biomarker; miRNA; transcription factor; interplay

### **1. Introduction**

Clear cell renal cell carcinoma (ccRCC) is the most common malignant tumor subtype of kidney cancer [1], which still remains a higher mortality rate in worldwide [2]. At present, the main treatment method on ccRCC patients is early resection, but its curative effect and prognosis are not very good for these terminally ccRCC patients [3]. Currently, above 30–50% of ccRCC patients have missed the best surgical opportunity due to the lack of early clinical symptoms [3]. Therefore, acquiring new molecular biomarkers not only are urgently needed for establishing clinically stratifying system to improve the diagnostic efficiency of ccRCC patients, but are of great clinical value for effectively improving the management and treatment strategy of ccRCC patients.

MicroRNAs (miRNAs), a class of small noncoding RNAs with about 22 nucleotides, can negatively regulate gene expression at the post-transcriptional level by influencing the mRNA degradation and/or translational efficiency involved in manifold biological processes [4–7]. In recent years, many studies have revealed that miRNAs not only can act as oncogenes or tumor suppressor genes involve in tumorigenesis and progressions of various cancers [8], but can serve as valuable boimarkers for the detection and prognosis of cancer patients [9–15]. Thus, the current study has been focused on finding novel miRNAs as effective prognostic predictors for the overall survival of cancer patients [16–20].

At present, some studies have suggested that a lot of miRNAs could act as oncogenic miRNAs or tumor suppressor miRNAs to participate in the tumorigenesis and progressions of ccRCC, and they might serve as diagnosic and prognosic biomarkers for ccRCC patients [21–24]. However, several studies have shown that some miRNAs may play complete contradictor roles in diagnosic and prognosic effects for ccRCC, such as miR-99a [25,26], miR-106a [27,28], miR-125b [21,29], miR-144 [30,31], miR-203 [32,33], and miR-378 [34,35], which might greatly limit their applications into the clinical diagnosis and prognosis for ccRCC patients. Therefore, we should further study the functional roles of miRNAs as potentially diagnostic and prognostic biomarkers for ccRCC, whilst it is also very necessary for expanding the screening of new reliable miRNAs as diagnostic and prognostic biomarkers for ccRCC.

Nowadays, a number of studies are mainly focused on the fuctional mechanism of miRNAs that regulate their targets to cause the occurrence and development of ccRCC. However, how are miRNAs themselves regulated in the occurrence of ccRCC, which is still largely unknown. Intriguingly, currently, many studies have shown that transcription factors (TFs) can regulate miRNA expressions, and miRNAs may also regulate TF expressions in gene regulatory networks [36–39], and TFs and miRNAs interplay could precisely modulate gene expressions to maintain cell homeostasis [37,40,41]. Therefore, considering the complexity of TFs and miRNAs interplay mediating gene expression, substantial works should be further made in elucidating the mechanism that TFs and miRNAs interplay drives the occurrence and development of ccRCC, and finding TFs and miRNAs as novel diagnostic and prognostic biomarkers in the survival of ccRCC patients.

Considering that miRNAs often carry out their functions through fine-tuning the expression of their target genes, and multiple miRNAs can also synergistically or antagonistically regulate one or more target genes to control the strength and duration of cell response [42,43]. Here, we firstly screened eight potentially prognostic miRNAs based on RNA-seq and clinical information from the TCGA database. We next combined the eight miRNAs as an integrative prognostic predictor to evaluate the prognostic efficiency. This result showed that the prognostic efficiency and credibility of the eight-miRNA signature significantly outperformed a single miRNA, which implies that the synergistical regulation of miRNAs plays key roles in the tumorigenesis and progression of ccRCC. Subsequently, we utilized target prediction software and KEGG enrichment analysis to find that these eight miRNAs mainly control a certain number of oncogenic and onco-suppressive genes from some crucially cancer-related pathways improving the survival of ccRCC patients. Additionally, we found that twenty-two TFs could interact with eight miRNAs based on deepCAGE, TransmiR v2.0, and MirWalk3.0 database [44–46]. To further reveal the possible molecular mechanism that TFs and miRNAs interplay facilitates the survival of ccRCC patients, we further constructed a interplay network of TFs and miRNAs and the network analysis revealed that the interplay between twenty-two TFs and eight miRNAs could synergistically control the expression of oncogenes, driver genes, and tumor suppressor genes to involve in the regulation of tumorigenesis and progression of ccRCC. Finally, we performed cox regression analysis to identify eleven genes as a eleven-gene signature, including six miRNAs and five TFs and the prognostic effect and the credibility of the eleven-gene signature also obviously outperformed the eight-miRNA signature and the nine-TF signature. Taken together, our findings not only revealed a novel possible mechanism that TFs and miRNAs interplay could regulate cooperatively oncogenes, driver genes, and tumor suppressor genes to facilitate the survival of ccRCC patients, but also identified some new potential prognostic factors and therapeutic targets for ccRCC patients.

#### **2. Results**

#### *2.1. Identification of miRNAs as Potential Prognostic Biomarkers*

In this work, we found 110 differentially expressed miRNAs between 480 ccRCC tissues and 68 paracancerous tissues, including 50 up-regulated and 60 down-regulated miRNAs (Figure S1A). Here, to identify these prognostic miRNAs for predicting the overall survival in ccRCC patients, we grouped 480 patients with at least 90 days into high- and low-expression groups according to

the median expression level of each of 110 differentially expressed miRNAs. We next performed survival analysis to find that 37 miRNAs might be associated with the survival of ccRCC patients (*p*-value < 0.05) (Table S1). Next, we performed univariate analysis and screened these top 21 miRNAs with a *p*-value < 0.001 (Figure 1, Table S2) for multivariate stepwise cox regression analysis to further determine independently prognostic miRNA biomarkers for ccRCC patients. We finally obtained eight potential prognostic miRNAs, including four down-regulated miRNAs (hsa-miR-9-5p, hsa-miR-1269a, hsa-miR-183-5p, hsa-miR-335-3p) and four up-regulated miRNAs (hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-144-5p, hsa-miR-3613-5p), which might be involved in the survival of ccRCC patients. We divided 480 ccRCC patients into high-risk and low-risk groups according to the median univariate cox risk score of each of these eight miRNAs to further detect the association between these identified eight miRNAs and the overall survival of ccRCC patients. Kaplan-Meier survival analysis and log-rank test indicated that, in all eight independent miRNA cohorts, ccRCC patients with high-risk groups exhibited the overall survival more badly than low-risk groups (all cohorts *p*-value < 0.01) (Figure S2). Interestingly, the range of the AUC value for eight miRNAs was about 0.6~0.7 (Figure S3), which indicated that the established prognosis model has a very good prognosis effect. Overall, these findings suggested that any of eight miRNAs might act as a possible prognostic biomarker for the survival of ccRCC patients. median expression level of each of 110 differentially expressed miRNAs. We next performed survival analysis to find that 37 miRNAs might be associated with the survival of ccRCC patients (*p*-value < 0.05) (Table S1). Next, we performed univariate analysis and screened these top 21 miRNAs with a *p*value < 0.001 (Figure 1, Table S2) for multivariate stepwise cox regression analysis to further determine independently prognostic miRNA biomarkers for ccRCC patients. We finally obtained eight potential prognostic miRNAs, including four down-regulated miRNAs (hsa-miR-9-5p, hsamiR-1269a, hsa-miR-183-5p, hsa-miR-335-3p) and four up-regulated miRNAs (hsa-miR-365b-3p, hsamiR-223-3p, hsa-miR-144-5p, hsa-miR-3613-5p), which might be involved in the survival of ccRCC patients. We divided 480 ccRCC patients into high-risk and low-risk groups according to the median univariate cox risk score of each of these eight miRNAs to further detect the association between these identified eight miRNAs and the overall survival of ccRCC patients. Kaplan-Meier survival analysis and log-rank test indicated that, in all eight independent miRNA cohorts, ccRCC patients with highrisk groups exhibited the overall survival more badly than low-risk groups (all cohorts *p-*value < 0.01) (Figure S2). Interestingly, the range of the AUC value for eight miRNAs was about 0.6~0.7 (Figure S3), which indicated that the established prognosis model has a very good prognosis effect. Overall, these findings suggested that any of eight miRNAs might act as a possible prognostic biomarker for the survival of ccRCC patients.

*Cancers* **2019**, *11*, x FOR PEER REVIEW 3 of 22

**Figure 1. The forest plots of hazard ratios (HR) of top 21 most significant survival associated miRNAs in ccRCC (***p***-value < 0.001)**. Red represents the risk factors (HR > 1) and blue represents protective factors (HR < 1). **Figure 1. The forest plots of hazard ratios (HR) of top 21 most significant survival associated miRNAs in ccRCC (***p***-value** < **0.001).** Red represents the risk factors (HR > 1) and blue represents protective factors (HR < 1).

#### *2.2. The Exprssion Level of Eight Prognostic miRNAs is Associated with the Survival of ccRCC Patients 2.2. The Exprssion Level of Eight Prognostic miRNAs is Associated with the Survival of ccRCC Patients*

To explore whether the eight miRNAs could be used as potential diagnostic biomarkers for distinguishing patients with ccRCC from controls. Here, we used survival curves to assess the association between the expression level of these eight miRNAs and the overall survival of ccRCC patients. We divided ccRCC patients into high-expression and low-expression miRNA groups according to the median expression level of each of these eight miRNAs. This result indicated that, in all eight independent miRNA cohorts, ccRCC patients with high-expression miRNA groups exhibited a worse overall survival rate than the low-expression miRNA groups, except for hsa-miR-144-5p (Figure 2). These results seemed to indicate that these highly expressed hsa-miR-9-5p, hsamiR-1269a, hsa-miR-183-5p, hsa-miR-335-3p, hsa-miR-365b-3p, hsa-miR-223-3p, and hsa-miR-3613- 5p might be associated with a poor prognosis for ccRCC patients. It is noteworthy that the highly To explore whether the eight miRNAs could be used as potential diagnostic biomarkers for distinguishing patients with ccRCC from controls. Here, we used survival curves to assess the association between the expression level of these eight miRNAs and the overall survival of ccRCC patients. We divided ccRCC patients into high-expression and low-expression miRNA groups according to the median expression level of each of these eight miRNAs. This result indicated that, in all eight independent miRNA cohorts, ccRCC patients with high-expression miRNA groups exhibited a worse overall survival rate than the low-expression miRNA groups, except for hsa-miR-144-5p (Figure 2). These results seemed to indicate that these highly expressed hsa-miR-9-5p, hsa-miR-1269a, hsa-miR-183-5p, hsa-miR-335-3p, hsa-miR-365b-3p, hsa-miR-223-3p, and hsa-miR-3613-5p might be associated with a poor prognosis for ccRCC patients. It is noteworthy

that the highly expressed hsa-miR-144-5p was associated with better overall survival (Figure 2), which suggested that hsa-miR-144-5p should be a good prognostic factor for ccRCC patients. In addition, we further calculated the association between the expression level of eight miRNAs and patient's clinical diagnostic factors, respectively. Our results showed that eight prognostic miRNAs were significantly associated with T stage, M stage, G stage, and pathologic stage (Figure S4, Table S3), implying that these eight miRNAs might be involved in tumorigenesis and progression of ccRCC and could be served as prognostic biomarkers for the survival of ccRCC patients. *Cancers* **2019**, *11*, x FOR PEER REVIEW 4 of 22 calculated the association between the expression level of eight miRNAs and patient's clinical diagnostic factors, respectively. Our results showed that eight prognostic miRNAs were significantly associated with T stage, M stage, G stage, and pathologic stage (Figure S4, Table S3), implying that these eight miRNAs might be involved in tumorigenesis and progression of ccRCC and could be served as prognostic biomarkers for the survival of ccRCC patients.

**Figure 2. Kaplan Meier survival based on the expression level of eight miRNAs.** Overall survival curves for high-expression and low-expression ccRCC patient cohorts. (**A**): hsa-miR-1269a; (**B**): hsamiR-183-5p; (**C**): hsa-miR-9-5p; (**D**): hsa-miR-335-3p; (**E**): hsa-miR-3613-5p; (**F**): hsa-miR-223-3p; (**G**): hsa-miR-365b-3p; and, (**H**): hsa-miR-144-5p. **Figure 2. Kaplan Meier survival based on the expression level of eight miRNAs.** Overall survival curves for high-expression and low-expression ccRCC patient cohorts. (**A**): hsa-miR-1269a; (**B**): hsa-miR-183-5p; (**C**): hsa-miR-9-5p; (**D**): hsa-miR-335-3p; (**E**): hsa-miR-3613-5p; (**F**): hsa-miR-223-3p; (**G**): hsa-miR-365b-3p; and, (**H**): hsa-miR-144-5p.

#### *2.3. Prognostic Value of Combined Eight miRNAs as a Signature in ccRCC Patients Cancers* **2019**, *11*, x FOR PEER REVIEW 5 of 22

Considering multiple miRNAs could synergistically or antagonistically regulate one or more target genes to control cell fate. Herein, we further combined these eight miRNAs as an integrative prognostic predictor. The 480 patients were divided into low-risk group and high-risk group and then subjected to survival analysis. Our results showed that there was a significant difference in the overall survival between the two risk groups, and the high-risk group had more worse overall survival than the low-risk group (*p* < 0.0001, Figure 3A). In addition, the ROC curve based on the eight-miRNA signature also, respectively, showed an average 3, 5, and 10 year AUC for 0.762, 0.747, and 0.746 (Figure 3B). Interestingly, the concordance index (0.7305) of the combined prognostic model of the eight-miRNA signature was higher than that of all single miRNA (Table S4), whereas the Akaike information criterion (1622.2491) of the combined prognostic model of the eight-miRNA signature was lower than that of all single miRNA (Table S4), which indicated that the prognostic effect and credibility of the eight-miRNA signature were clearly superior to all single miRNA. These findings suggested that the eight-miRNA signature could act as a prognostic biomarker for promoting the survival of ccRCC patients. Considering multiple miRNAs could synergistically or antagonistically regulate one or more target genes to control cell fate. Herein, we further combined these eight miRNAs as an integrative prognostic predictor. The 480 patients were divided into low-risk group and high-risk group and then subjected to survival analysis. Our results showed that there was a significant difference in the overall survival between the two risk groups, and the high-risk group had more worse overall survival than the low-risk group (*p* < 0.0001, Figure 3A). In addition, the ROC curve based on the eight-miRNA signature also, respectively, showed an average 3, 5, and 10 year AUC for 0.762, 0.747, and 0.746 (Figure 3B). Interestingly, the concordance index (0.7305) of the combined prognostic model of the eight-miRNA signature was higher than that of all single miRNA (Table S4), whereas the Akaike information criterion (1622.2491) of the combined prognostic model of the eight-miRNA signature was lower than that of all single miRNA (Table S4), which indicated that the prognostic effect and credibility of the eight-miRNA signature were clearly superior to all single miRNA. These findings suggested that the eight-miRNA signature could act as a prognostic biomarker for promoting the survival of ccRCC patients.

Cox proportional hazard regression analysis was further used to characterize the impact of various clinical factors on the overall survival of ccRCC patients (Table 1). Age, gender, tumor size, metastasis, pathologic stage, neoplasm histologic grade, and the combined eight miRNAs signature were coded as continuous variables. As shown in Table 1, the univariate analysis showed that all factors, except for gender, might act as prognostic indicators for ccRCC patients. However, the multivariate analysis indicated that only age and the eight-miRNA signature can be used as independent prognostic indicators for ccRCC patients. This result revealed that the eight-miRNA signature could not only could serve as an independent prognostic factor for overall survival of ccRCC patients, but also act as an effective risk stratification indicator for ccRCC patient diagnosis. Cox proportional hazard regression analysis was further used to characterize the impact of various clinical factors on the overall survival of ccRCC patients (Table 1). Age, gender, tumor size, metastasis, pathologic stage, neoplasm histologic grade, and the combined eight miRNAs signature were coded as continuous variables. As shown in Table 1, the univariate analysis showed that all factors, except for gender, might act as prognostic indicators for ccRCC patients. However, the multivariate analysis indicated that only age and the eight-miRNA signature can be used as independent prognostic indicators for ccRCC patients. This result revealed that the eight-miRNA signature could not only could serve as an independent prognostic factor for overall survival of ccRCC patients, but also act as an effective risk stratification indicator for ccRCC patient diagnosis.

**Figure 3. Kaplan Meier survival and receiver operating characteristic (ROC) curves based on the riskscore of the eight-miRNA signature.** (**A**): Overall survival curves of high-risk and low-risk based on the eight-miRNA signature model. (**B**): Receiver operating characteristic (ROC) curves for high **Figure 3. Kaplan Meier survival and receiver operating characteristic (ROC) curves based on the riskscore of the eight-miRNA signature.** (**A**): Overall survival curves of high-risk and low-risk based on the eight-miRNA signature model. (**B**): Receiver operating characteristic (ROC) curves for high and low risk from the eight-miRNA signature model.

and low risk from the eight-miRNA signature model.


**Table 1.** Univariate and multivariate Cox regression analysis of the overall survival for clinical factors and risk of the combined eight prognostic miRNAs as a signature.

Age, gender, tumor stage, metastasis pathologic, pathologic stage, histologic grade, and the eight-miRNA signature were coded as continuous variable. Specifically, pathologic stage was coded as I = 1, II = 2, III = 3, IV = 4. Tumor stage was coded as T1 = 1, T2 = 2, T3 = 3, T4 = 4. Histologic grade was coded as G1 = 1, G2 = 2, G3 = 3, G4 = 4.

#### *2.4. Function Roles of Eight Prognostic miRNAs in ccRCC*

To reveal the functional role of eight prognostic miRNAs in ccRCC, we first identified 3672 differentially expressed genes between 480 ccRCC tissues and 68 paracancerous tissues, including 2057 up-regulated genes (of which 116 transcription factors) and 1012 down-regulated genes (of which 61 transcription factors) (Figure S1B). Next, we obtained 357 down-regulated targets of 4 up-regulated miRNAs (hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-144-5p, hsa-miR-3613-5p) and 1012 up-regulated targets of 4 down-regulated miRNAs (hsa-miR-9-5p, hsa-miR-1269a, hsa-miR-183-5p, hsa-miR-335-3p).

To elucidate the function of these target genes of these eight miRNAs in ccRCC, we further performed KEGG pathway analysis using clusterProfiler R package. We found that these significantly up-regulated target genes of four down-regulated miRNAs were widely involved in cancer-related signaling pathways, such as MAPK, Ras, NF-kappa B, Chemokine and Cytokine-cytokine receptor (Figure S5A). These results suggested that the interplay among multiple signaling pathways might synergistically mediate the occurrence and development of ccRCC.

We further picked out these up-regulated target genes of four down-regulated miRNAs from these main cancer signaling pathways to construct the miRNA-gene regulation network (Figure 4A). As shown in Figure 4A, some up-regulated target genes were regulated by more than one down-regulated miRNA, implying that the cooperative regulation of multiple miRNAs might play key roles in the initiation and progression of ccRCC. In addition, the protein-protein interaction (PPI) network was also constructed for these up-regulated target genes of four down-regulated miRNAs using the STRING database, which demonstrated a close interaction within these target genes (Figure 4B). Here, a node with ≥ 20 degrees is defined as a hub gene, thus we found 30 hub genes (Table S5). We next used multivariate cox regression analysis for 30 hub genes to further select prognosis-related genes for ccRCC patients. We found ten potential prognostic genes, including eight tumor suppressor and/or driver genes (*VEGFA*, *CCND1*, *BAX*, *IL7R*, *SHC1*, *FLT1*, *IL7,* and *JAK3*) [47–51], as well as two chemokines (*CXCL9* and *CXCL10*) [50]. Additionally, we combined the above ten hub genes as a signature to perform survival analysis. The result indicated that the low-risk group has a better overall survival rate than the high-risk group (*p* < 0.0001, Figure 4C), whilst the ROC curve also demonstrated that the ten-hub gene signature had a better prognostic effect and credibility for ccRCC patients with an average 3, 5 and 10 year AUC for 0.728, 0.751 and 0.796, respectively (Figure 4D). Taken together, our present findings suggested a possible molecular mechanism that the down-regulated expressed hsa-miR-9-5p, hsa-miR-1269a, hsa-miR-183-5p, and hsa-miR-335-3p might cooperatively up-regulate the expression level of numerous tumor suppressor and/or driver genes from some cancer-related pathways to improve the overall survival of ccRCC patients.

*Cancers* **2019**, *11*, x FOR PEER REVIEW 7 of 22

**Figure 4. The cancer signaling pathway and protein network analysis of these up-regulated targets and the prognostic model of the ten-hub gene signature.** (**A**): Target genes of four down-regulated miRNAs involve in the network map of cancer-related signaling pathways. The green represents the down-regulated expression and the red represents up-regulated expression; (**B**): Target protein interaction network of four down-regulated miRNAs. The blue line means low credibility and the purple line means high credibility; (**C**): Overall survival curves for high-risk and low-risk groups based on the ten-hub gene model; and, (**D**): ROC curves for high-risk and low-risk based on the ten-**Figure 4. The cancer signaling pathway and protein network analysis of these up-regulated targets and the prognostic model of the ten-hub gene signature**. (**A**): Target genes of four down-regulated miRNAs involve in the network map of cancer-related signaling pathways. The green represents the down-regulated expression and the red represents up-regulated expression; (**B**): Target protein interaction network of four down-regulated miRNAs. The blue line means low credibility and the purple line means high credibility; (**C**): Overall survival curves for high-risk and low-risk groups based on the ten-hub gene model; and, (**D**): ROC curves for high-risk and low-risk based on the ten-hub gene model.

hub gene model. Compared to four down-regulated miRNAs, these down-regulated target genes of four up-Compared to four down-regulated miRNAs, these down-regulated target genes of four up-regulated miRNAs were less enriched in cancer-associated signaling pathways (Figure S5B). However, we

regulated miRNAs were less enriched in cancer-associated signaling pathways (Figure S5B).

further analyzed these target genes of three up-regulated miRNAs (hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-3613-5p), finding that many genes of their target genes are involved in some cancer-associated signaling pathways (Figure 5A). Therefore, here, we reused multivariate cox regression analysis to screen potential prognostic genes for ccRCC patients. The result showed that five genes (*PRKCA*, *ADORA1*, *PPARGC1A*, *KL*, *GNG7*) could be identified as potentially prognostic factors, and they have also been suggested as tumor suppressor genes [49]. Interestingly, the survival analysis indicated that the five-gene signature could significantly stratify ccRCC patients into a high- and low-risk group (*p* < 0.0001, Figure 5B), and the AUC value of an average 3, 5 and 10 year is 0.696, 0.698, and 0.708, respectively (Figure 5C). Overall, we proposed a possibly functional mechanism that three highly expressed miRNAs (hsa-miR-365b-3p, hsa-miR-223-3p, and hsa-miR-3613-5p) might synergistically down-regulate the expression of many tumor suppressor genes to decrease the survival of ccRCC patients. However, we further analyzed these target genes of three up-regulated miRNAs (hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-3613-5p), finding that many genes of their target genes are involved in some cancer-associated signaling pathways (Figure 5A). Therefore, here, we reused multivariate cox regression analysis to screen potential prognostic genes for ccRCC patients. The result showed that five genes (*PRKCA*, *ADORA1*, *PPARGC1A*, *KL*, *GNG7*) could be identified as potentially prognostic factors, and they have also been suggested as tumor suppressor genes [49]. Interestingly, the survival analysis indicated that the five-gene signature could significantly stratify ccRCC patients into a highand low-risk group (*p* < 0.0001, Figure 5B), and the AUC value of an average 3, 5 and 10 year is 0.696, 0.698, and 0.708, respectively (Figure 5C). Overall, we proposed a possibly functional mechanism that three highly expressed miRNAs (hsa-miR-365b-3p, hsa-miR-223-3p, and hsa-miR-3613-5p) might synergistically down-regulate the expression of many tumor suppressor genes to decrease the survival of ccRCC patients.

*Cancers* **2019**, *11*, x FOR PEER REVIEW 8 of 22

**Figure 5. The heat map of down-regulated targets of three up-regulated miRNAs involve in cancerrelated signaling pathway and the prognostic model of the five-gene signature.** (**A**): The heat map of these down-regulated targets regulated by has-miR-3613-5p, has-miR-223-3p, and has-miR-365b-3p involved in cancer-related signaling pathways; (**B**): Overall survival curves for high-risk and lowrisk groups based on the prognostic model of a five-gene signature; (**C**): ROC curves for high and low risk from the prognostic model of a five-gene signature. **Figure 5. The heat map of down-regulated targets of three up-regulated miRNAs involve in cancer-related signaling pathway and the prognostic model of the five-gene signature**. (**A**): The heat map of these down-regulated targets regulated by has-miR-3613-5p, has-miR-223-3p, and has-miR-365b-3p involved in cancer-related signaling pathways; (**B**): Overall survival curves for high-risk and low-risk groups based on the prognostic model of a five-gene signature; (**C**): ROC curves for high and low risk from the prognostic model of a five-gene signature.

Interestingly, our above results have indicated that the highly expressed hsa-miR-144-5p could act as a good prognostic factor for ccRCC patients. How does the highly expressed hsa-miR-144-5p facilitate the survival of ccRCC patients? Thus, we carried out a in-depth analysis for these target genes of hsa-miR-144-5p. Intriguingly, we found that hsa-miR-144-5p could regulate 55 genes, and, in particular, twelve genes of them have been reported as oncogenes or driver genes (Table S6) [49]. Thus, we further used multivariate cox regression analysis for twelve oncogenes and driver genes to screen potential prognostic factors for ccRCC patients. Consequently, these five genes (*MAGI3*, *CDKL1*, *CDH1*, *PPM1K*, *MSI2*) could be served as potential prognostic factors. We further combined the five genes as an integrative prognostic predictor for survival analysis. These results showed that the high-risk group had a worse overall survival than the low-risk group (*p* < 0.0001, Figure 6A), and the AUC value of an average 3, 5 and 10 year is 0.659, 0.676 and 0.759, respectively, based on the ROC curve (Figure 6B). Collectively, our present findings implied that the highly expressed hsa-miR-144-5p might facilitate the overall survival of ccRCC patients through down-regulating the expression level of some certain oncogenes and/or driver genes. Interestingly, our above results have indicated that the highly expressed hsa-miR-144-5p could act as a good prognostic factor for ccRCC patients. How does the highly expressed hsa-miR-144-5p facilitate the survival of ccRCC patients? Thus, we carried out a in-depth analysis for these target genes of hsa-miR-144-5p. Intriguingly, we found that hsa-miR-144-5p could regulate 55 genes, and, in particular, twelve genes of them have been reported as oncogenes or driver genes (Table S6) [49]. Thus, we further used multivariate cox regression analysis for twelve oncogenes and driver genes to screen potential prognostic factors for ccRCC patients. Consequently, these five genes (*MAGI3*, *CDKL1*, *CDH1*, *PPM1K*, *MSI2*) could be served as potential prognostic factors. We further combined the five genes as an integrative prognostic predictor for survival analysis. These results showed that the high-risk group had a worse overall survival than the low-risk group (*p* < 0.0001, Figure 6A), and the AUC value of an average 3, 5 and 10 year is 0.659, 0.676 and 0.759, respectively, based on the ROC curve (Figure 6B). Collectively, our present findings implied that the highly expressed hsa-miR-144- 5p might facilitate the overall survival of ccRCC patients through down-regulating the expression level of some certain oncogenes and/or driver genes.

**Figure 6. Kaplan Meier survival and ROC curves based on the risk score of five prognostic targets of has-miR-144-5p.** (**A**): Overall survival curves for high-risk and low-risk groups based on the prognostic model of a five-gene signature; (**B**): ROC curves for high and low risk from the prognostic model of a five-gene signature. **Figure 6. Kaplan Meier survival and ROC curves based on the risk score of five prognostic targets of has-miR-144-5p.** (**A**): Overall survival curves for high-risk and low-risk groups based on the prognostic model of a five-gene signature; (**B**): ROC curves for high and low risk from the prognostic model of a five-gene signature.

#### *2.5. The Interplay Network Between Twenty-Two TFs and Eight miRNAs*

*2.5. The Interplay Network Between Twenty-Two TFs and Eight miRNAs*  To explore the interplay between TFs and eight prognostic miRNAs, we further used TransmiR2.0 to predict these TFs that may regulate the eight prognostic miRNAs. Here, we found six down-regulated TFs (KLF5, SREBF2, TFAP2A, HIF1A, GATA3, and GATA2), which could regulate three down-regulated miRNAs (hsa-miR-9-5p, hsa-miR-183-5p, and hsa-miR-335-3p), but no TF was found for regulating hsa-miR-1269a (Figure S6A). Whilst, 16 up-regulated TFs (CEBPA, E2F1, FOXM1, FLI1, HEY1, IKZF1, IRF1, MEF2C, MYC, POU2F2, POU5F1, PRDM1, RARA, RUNX1, RUNX3, and TCF4) were also found to regulate four up-regulated miRNAs (hsa-miR-365b-3p, hsamiR-223-3p, hsa-miR-144-5p, hsa-miR-3613-5p) (Figure S6B). Furthermore, we also predicted these twenty-two TFs regulated by the eight prognostic miRNAs. These detailed TF-miRNA regulatory pairs were outlined in the Table S7. Based on the TF-miRNA pairs, we further constructed a TFmiRNA interplay network (Figure 7), which showed that these up-regulated TFs, such as MYC, IKZF1, and IRF1 could up-regulate the expression level of hsa-miR-223-3p and hsa-miR-365b-3p to To explore the interplay between TFs and eight prognostic miRNAs, we further used TransmiR2.0 to predict these TFs that may regulate the eight prognostic miRNAs. Here, we found six down-regulated TFs (KLF5, SREBF2, TFAP2A, HIF1A, GATA3, and GATA2), which could regulate three down-regulated miRNAs (hsa-miR-9-5p, hsa-miR-183-5p, and hsa-miR-335-3p), but no TF was found for regulating hsa-miR-1269a (Figure S6A). Whilst, 16 up-regulated TFs (CEBPA, E2F1, FOXM1, FLI1, HEY1, IKZF1, IRF1, MEF2C, MYC, POU2F2, POU5F1, PRDM1, RARA, RUNX1, RUNX3, and TCF4) were also found to regulate four up-regulated miRNAs (hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-144-5p, hsa-miR-3613-5p) (Figure S6B). Furthermore, we also predicted these twenty-two TFs regulated by the eight prognostic miRNAs. These detailed TF-miRNA regulatory pairs were outlined in the Table S7. Based on the TF-miRNA pairs, we further constructed a TF-miRNA interplay network (Figure 7), which showed that these up-regulated TFs, such as MYC, IKZF1, and IRF1 could up-regulate the expression level of hsa-miR-223-3p and hsa-miR-365b-3p to inhibit the expression of some TFs, such as GATA2 and SREBF2, then reducing the expression of hsa-miR-183-5p and hsa-miR-335-3p to

inhibit the expression of some TFs, such as GATA2 and SREBF2, then reducing the expression of hsa-

up-regulate tumor suppressor gene expression (Figure 7). Whilst these down-regulated TFs, such as KLF5, SREBF2, TFAP2A, HIF1A, GATA3, and GATA2, could also down-regulate the expression level of hsa-miR-9-5p, hsa-miR-183-5p, and hsa-miR-335-3p to up-regulate their target TF expression, such as E2F1, RUNX1, and RUNX3, and then up-regulated the expression of hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-144-5p, and hsa-miR-3613-5p to down-regulate the expression of some tumor suppressor genes or oncogenes (Figure 7). Specially, TF and miRNA carry out opposing functions. Therefore, our study seemed to imply that the interplay between twenty-two TFs and eight prognostic miRNAs might precisely control the expression of oncogenes, driver genes, and tumor suppressor genes to facilitate the survival of ccRCC patients. these down-regulated TFs, such as KLF5, SREBF2, TFAP2A, HIF1A, GATA3, and GATA2, could also down-regulate the expression level of hsa-miR-9-5p, hsa-miR-183-5p, and hsa-miR-335-3p to upregulate their target TF expression, such as E2F1, RUNX1, and RUNX3, and then up-regulated the expression of hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-144-5p, and hsa-miR-3613-5p to downregulate the expression of some tumor suppressor genes or oncogenes (Figure 7). Specially, TF and miRNA carry out opposing functions. Therefore, our study seemed to imply that the interplay between twenty-two TFs and eight prognostic miRNAs might precisely control the expression of oncogenes, driver genes, and tumor suppressor genes to facilitate the survival of ccRCC patients.

*Cancers* **2019**, *11*, x FOR PEER REVIEW 10 of 22

**Figure 7. The interplay network between twenty-two transcription factors and eight miRNAs.** The circle represents the miRNA and the diamond represents the transcription factor. The blue represents a down-regulation expression and the purple represents an up-regulation expression. The sharp **Figure 7. The interplay network between twenty-two transcription factors and eight miRNAs.** The circle represents the miRNA and the diamond represents the transcription factor. The blue represents a down-regulation expression and the purple represents an up-regulation expression. The sharp arrow represents activation and the flat arrow represents inhibition.

#### arrow represents activation and the flat arrow represents inhibition. *2.6. Prognostic Value of the Combined Nine TFs as a Signature in ccRCC Patients*

the survival of ccRCC patients.

*2.6. Prognostic Value of the Combined Nine TFs as a Signature in ccRCC Patients*  Based on these above results, we further used multivariate cox regression analysis for these above twenty-two TFs to determine independent prognostic TFs for ccRCC patients. Herein, we identified nine potential prognostic TFs (*TFAP2A*, *KLF5*, *IRF1*, *RUNX1*, *RARA*, *GATA3*, *IKZF1*, *POU2F2,* and *FOXM1*) that could act as prognostic factors for ccRCC patients. We next combined these nine TFs as a signature to detect the prognostic effect for ccRCC patients. The 480 patients were divided into low-risk group and high-risk group and subjected to survival analysis. Our results showed that the high-risk group had worse overall survival than the low-risk group (*p* < 0.0001, Figure 8A). In addition, the ROC curve based on the nine-TF pool also, respectively, showed an average 3, 5, and 10 year AUC for 0.721, 0.748, and 0.780 (Figure 8B), indicating that the nine-TF signature had very good prognosis effect and credibility for the survival of ccRCC patients. These Based on these above results, we further used multivariate cox regression analysis for these above twenty-two TFs to determine independent prognostic TFs for ccRCC patients. Herein, we identified nine potential prognostic TFs (*TFAP2A*, *KLF5*, *IRF1*, *RUNX1*, *RARA*, *GATA3*, *IKZF1*, *POU2F2,* and *FOXM1*) that could act as prognostic factors for ccRCC patients. We next combined these nine TFs as a signature to detect the prognostic effect for ccRCC patients. The 480 patients were divided into low-risk group and high-risk group and subjected to survival analysis. Our results showed that the high-risk group had worse overall survival than the low-risk group (*<sup>p</sup>* <sup>&</sup>lt; 0.0001, Figure 8A). In addition,the ROC curve based on the nine-TF pool also, respectively, showed an average 3, 5, and 10 year AUC for 0.721, 0.748, and 0.780 (Figure 8B), indicating that the nine-TF signature had very good prognosis effect and credibility for the survival of ccRCC patients. These findings suggested that the nine-TF signature could serve as a prognostic biomarker for improving the survival of ccRCC patients.

findings suggested that the nine-TF signature could serve as a prognostic biomarker for improving

**Figure 8. Kaplan Meier survival and ROC curves based on the risk score of the combined nine transcription factors as a signature.** (**A**): Overall survival curves for high-risk and low-risk groups based on the prognostic model of a nine-TF signature; (**B**): ROC curves for high and low risk from the prognostic model of a nine-TF signature. **Figure 8. Kaplan Meier survival and ROC curves based on the risk score of the combined nine transcription factors as a signature.** (**A**): Overall survival curves for high-risk and low-risk groups based on the prognostic model of a nine-TF signature; (**B**): ROC curves for high and low risk from the prognostic model of a nine-TF signature.

#### *2.7. Clinical Value of TFs and miRNAs Interplay as a Prognostic Signature in ccRCC Patients*

*2.7. Clinical Value of TFs and miRNAs Interplay as a Prognostic Signature in ccRCC Patients*  To further reveal the role of the interaction between TFs and miRNAs in the overall survival of ccRCC patients, we hypothesized that the interaction between TFs and miRNAs is likely to improve prognosis. Thus, we next performed multivariate analysis for above twenty-two TFs and eight miRNAs to identify prognostic TFs and miRNAs for ccRCC patients. We ultimately screened eleven potential prognostic factors, including hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-1269a, hsa-miR-144-5p, hsa-miR-183-5p, hsa-miR-335-3p, *TFAP2A*, *KLF5*, *IRF1*, *MYC*, *IKZF1*. We further combined the above eleven genes as a signature for survival analysis. We found that the high-risk group had worse overall survival than the low-risk group (*p* < 0.0001, Figure 9A). Whilst the ROC curve based on the eleven-gene signature also, respectively, showed an average 3, 5, and 10 year AUC values for 0.777, 0.771, and 0.785 (Figure 10B). Strikingly, the concordance index (0.7552) of the combined prognostic model of the eleven-gene signature was more higher than that of the eight-miRNA signature (0.7305) and the nine-TF signature (0.7281), and the Akaike information criterion (1606.0516) of the combined prognostic model of the eleven-gene signature was lower than that of the eight-miRNA signature (1622.2941) and the nine-TF signature (1632.5955) (Table S8), which indicated that the prognostic effect and the credibility of the eleven-gene signature were better than both the eight-miRNA signature and the nine-TF signature did. These findings suggested that the To further reveal the role of the interaction between TFs and miRNAs in the overall survival of ccRCC patients, we hypothesized that the interaction between TFs and miRNAs is likely to improve prognosis. Thus, we next performed multivariate analysis for above twenty-two TFs and eight miRNAs to identify prognostic TFs and miRNAs for ccRCC patients. We ultimately screened eleven potential prognostic factors, including hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-1269a, hsa-miR-144-5p, hsa-miR-183-5p, hsa-miR-335-3p, *TFAP2A*, *KLF5*, *IRF1*, *MYC*, *IKZF1*. We further combined the above eleven genes as a signature for survival analysis. We found that the high-risk group had worse overall survival than the low-risk group (*p* < 0.0001, Figure 9A). Whilst the ROC curve based on the eleven-gene signature also, respectively, showed an average 3, 5, and 10 year AUC values for 0.777, 0.771, and 0.785 (Figure 10B). Strikingly, the concordance index (0.7552) of the combined prognostic model of the eleven-gene signature was more higher than that of the eight-miRNA signature (0.7305) and the nine-TF signature (0.7281), and the Akaike information criterion (1606.0516) of the combined prognostic model of the eleven-gene signature was lower than that of the eight-miRNA signature (1622.2941) and the nine-TF signature (1632.5955) (Table S8), which indicated that the prognostic effect and the credibility of the eleven-gene signature were better than both the eight-miRNA signature and the nine-TF signature did. These findings suggested that the eleven-gene signature could act as a prognostic factor for the overall survival of ccRCC patients.

eleven-gene signature could act as a prognostic factor for the overall survival of ccRCC patients. Here, we also further used Cox proportional hazard regression analysis to characterize the impact of various clinical factors on overall survival of ccRCC patients (Table 2). Age, gender, tumor-pathologic, metastasis pathologic, pathologic stage, neoplasm histologic grade, and the eleven-gene signature were coded as continuous variables. The univariate analysis showed that all factors, except for gender, might serve as prognostic indicators for ccRCC patients (Table 2). Notably, the multivariate analysis demonstrated that only age and the eleven-gene signature could be used as independent prognostic indicators for ccRCC patients (Table 2). Taken together, our present results revealed that the interaction between TFs and miRNAs might have very important effects on the overall survival of ccRCC patients, and the eleven-gene signature might serve as an independent factor to improve the prognosis for ccRCC patients.

**Figure 9. Kaplan Meier survival and ROC curves based on the risk score of the combined five transcription factors and six miRNAs as a signature.** (**A**): Overall survival curves for high-risk and low-risk groups based on the prognostic model of an eleven-gene signature; and, (**B**): ROC curves for **Figure 9. Kaplan Meier survival and ROC curves based on the risk score of the combined five transcription factors and six miRNAs as a signature.** (**A**): Overall survival curves for high-risk and low-risk groups based on the prognostic model of an eleven-gene signature; and, (**B**): ROC curves for high and low risk from the prognostic model of an eleven-gene signature.



factor to improve the prognosis for ccRCC patients. Overall, herein we propose a potential molecular mechanism that the interplay between TFs and miRNAs facilitates the overall survival of ccRCC patients (Figure 10). On the one hand, down-Age, gender, tumor stage, metastasis pathologic, pathologic stage, histologic grade and the eleven-gene signature were coded as continuous variable. Specifically, pathologic stage was coded as I = 1, II = 2, III = 3, IV = 4. Tumor stage was coded as T1 = 1, T2 = 2, T3 = 3, T4 = 4. Histologic grade was coded as G1 = 1, G2 = 2, G3 = 3, G4 = 4.

regulated expressed TFs could down-regulate some miRNA expression to further up-regulate tumor

suppressor gene expression, whilst the down-regulated miRNAs could also up-regulate TF expression to further up-regulate the expression of tumor suppressor miRNA to down-regulate the expression of some oncogenes to improve prognosis for ccRCC. On the other hand, up-regulated TFs could up-regulate miRNA expression to inhibit the expression of some tumor suppressor genes, whilst the up-regulated miRNAs could also down-regulate the expression of some TFs to further down-regulate the expression of some miRNAs to up-regulate the expression of some tumor suppressor genes, which might promote ccRCC development. Taken together, our results seemed to reveal that the interplay between TFs and miRNAs might synergistically regulate the expression of a certain number of oncogenes, driver genes, and tumor suppressor genes to improve prognosis for ccRCC patients in transcriptional and post-transcriptional levels. Overall, herein we propose a potential molecular mechanism that the interplay between TFs and miRNAs facilitates the overall survival of ccRCC patients (Figure 10). On the one hand, down-regulated expressed TFs could down-regulate some miRNA expression to further up-regulate tumor suppressor gene expression, whilst the down-regulated miRNAs could also up-regulate TF expression to further up-regulate the expression of tumor suppressor miRNA to down-regulate the expression of some oncogenes to improve prognosis for ccRCC. On the other hand, up-regulated TFs could up-regulate miRNA expression to inhibit the expression of some tumor suppressor genes, whilst the up-regulated miRNAs could also down-regulate the expression of some TFs to further down-regulate the expression of some miRNAs to up-regulate the expression of some tumor suppressor genes, which might promote ccRCC development. Taken together, our results seemed to reveal that the interplay between TFs and miRNAs might synergistically regulate the expression of a certain number of oncogenes, driver genes, and tumor suppressor genes to improve prognosis for ccRCC patients in transcriptional and post-transcriptional levels.

as G1 = 1, G2 = 2, G3 = 3, G4 = 4.

**Table 2.** Univariate and multivariate Cox regression analysis of overall survival for clinical factors and risk of the combined five prognostic transcription factors and six miRNAs as a signature.

**Variables Univariate Analysis Multivariate Analysis**

Age 1.029 (1.014–1.043) <0.001 1.028 (1.011–1.044) <0.001 Gender 0.920 (0.6595–1.284) 0.624 0.949 (0.673–1.338) 0.765 Tumor\_pathologic\_T 1.855 (1.556–2.211) <0.001 0.840 (0.539–1.311) 0.443 Metastasis\_pathologic\_M 4.671 (3.369–6.475) <0.001 1.707 (0.830–3.513) 0.146 Pathologic\_stage\_Stage 1.884 (1.634–2.172) <0.001 1.410 (0.865–2.300) 0.178 Histologic\_grade\_G 2.238 (1.800–2.783) <0.001 1.420 (1.118–1.803) 0.765 The eleven-gene signature 4.349 (2.943–6.428) <0.001 2.590 (1.706–3.930) <0.001 Age, gender, tumor stage, metastasis pathologic, pathologic stage, histologic grade and the elevengene signature were coded as continuous variable. Specifically, pathologic stage was coded as I = 1, II

**Hazard Ratio (95% CI)** *p***-Value Hazard Ratio (95% CI)** *p***-Value** 

**Figure 10. The molecular mechanism of the interplay between transcription factors and miRNAs improving the prognosis of ccRCC patients.** The red represents an up-regulation of gene expression and the green represents down-regulation of gene expression. The sharp arrow represents activation and the flat arrow represents inhibition. **Figure 10. The molecular mechanism of the interplay between transcription factors and miRNAs improving the prognosis of ccRCC patients.** The red represents an up-regulation of gene expression and the green represents down-regulation of gene expression. The sharp arrow represents activation and the flat arrow represents inhibition.

#### **3. Discussion**

At present, a number of miRNAs have been suggested as diagnosic and prognosic biomarkers for ccRCC patients, but, since a single miRNA mainly fine-tunes gene expression to execute its function, and the functional effect of a single miRNA is relatively weak, which might result in a lack of miRNA application into the clinical diagnosis and prognosis for ccRCC patients. Many studies have shown that miRNAs are more likely to regulate a certain number of gene expressions to control cell fate [52,53]. Remarkably, several studies have revealed that TFs can regulate miRNA expressions, whilst miRNAs may also regulate TF expressions, and TFs and miRNAs interplay can precisely modulate gene expression in transcriptional and post-transcriptional levels [36,45,54–57]. However, the regulatory landscape by miRNAs and TFs interplay is still largely unknown in the tumorigenesis and progression of ccRCC up to now. Especially, the interplay network of miRNA and TF has not yet been systemically studied in ccRCC. Therefore, in this current work, we have integratedly analyzed miRNA, TF, and mRNA profilings, and identified some potential diagnostic and prognostic factors participated in the

survival of ccRCC patients, as well as revealed a possible molecular mechanism that miRNA and TF interplay can serve as an effective prognostic factor to facilitate the survival of ccRCC patients.

In this study, we have identified eight potentially diagnostic and prognostic miRNAs that could significantly distinguish the survival and pathological stratification for ccRCC patients. Among them, four down-regulated miRNAs (hsa-miR-9-5p, hsa-miR-1269a, hsa-miR-183-5p, hsa-miR-335-3p) could serve as good prognostic factors for clinical application of ccRCC patients, conversely the three up-regulated hsa-miR-365b-3p, hsa-miR-223-3p and hsa-miR-3613-5p of them might been acted as poor prognostic factors for ccRCC patients (Figure 2). At present, a few studies have reported that hsa-miR-183-5p can promote canceration by targeting SRSF2 in renal cancer [58], miR-335 could inhibit the proliferation and invasion of clear renal cells through suppressing Bcl-w [59], hsa-miR-365b-3p is poorly associated with ccRCC patient survival [60] and hsa-miR-9-5p is associated with the development and risk of renal cancer recurrence [61]. Although studies on the function roles of the other three miRNAs in ccRCC are still less reported to date, some other studies have revealed that hsa-miR-1269a could function as an onco-miRNA in NSCLC via down-regulating its target SOX6 [62], as well as that miR-1269a could promote colorectal cancer (CRC) metastasis by targeting Smad7 and HOXD10 [63]. The highly expressed hsa-miR-3613 as a oncogene could inhibit apoptosis via the down-regulation of target *APAF1* in human neuroblastoma BE(2)-C cells, as well as serve as a potential prognostic biomarker for pancreatic adenocarcinoma [64]. Some recent studies have indicated that miR-223-3p could down-regulate aNHEJ expression to result in synthetic lethality in human BRCA1-deficient cancers [65], and also act as an oncogenic miRNA in colon cancer through regulating EMT and PRDM1 [66]. These previous studies have revealed that the high expression of the seven miRNAs could down-regulate the respective tumor suppressor genes or driver genes involved in different malignant tumor occurrences. Especially, the same miRNA might regulate different target genes involved in specific tumor formation, which suggests that the different roles of same miRNA might depend on different cell microenvironment and cancer types. Interestingly, our present study has indicated that the seven aforementioned miRNAs mainly regulate tumor suppressor and/or driver genes to involve in ccRCC. Thus, we suggest that the seven miRNAs could served as potentially diagnostic and prognostic factors that may be due to the four down-regulated miRNAs that could restore or elevate the expression level of a certain number of tumor suppressor or driver genes to improve the survival of ccRCC patients, whereas the three up-regulated miRNAs might decrease the expression level of many tumor suppressor genes to result in the worse survival of ccRCC patients. Particularly, our present findings have shown that the highly expressed hsa-miR-144-5p could reduce the expression level of multiple oncogenic genes to promote the survival for ccRCC prognosis. A previous study has also shown that the hsa-miR-144-5p could serve as a tumor-suppressor gene for inhibiting cell growth and arresting cells in the G1 phase in renal cancer [30], which is also agreement with our present result. Taken together, we urge that hsa-miR-223-3p, hsa-miR-365b-3p, hsa-miR-3613-5p, hsa-miR-9-5p, hsa-miR-183-5p, hsa-miR-335-3p, and hsa-miR-1269a could serve as oncogenic miRNAs, as well as the hsa-miR-144-5p could act as a tumor suppressor miRNA for ccRCC patient diagnosis, but it is still further experimental verification.

A single miRNA molecule is known to carry out its function through fine-tuning the expression of target genes [42,43]. Therefore, it is hard to say that the expression of a single gene regulated by a single miRNA could significantly impact the proliferation and migration of cancer cell. Here, our findings seem to imply that single miRNA might need to regulate the expression of a lot of genes involved in the progression of ccRCC, suggesting that a single miRNA might be more suitable for ccRCC diagnosis and prognosis than a single gene. Specially, multiple miRNAs prefer synergistically or antagonistically regulating one or more target genes to control the strength and duration of cell response [42,43]. Thus, the combined multi-miRNAs as diagnosic and prognosic factors might be more suitable for the clinical application of ccRCC patients than a single miRNA. In the present study, just as we wish the prognostic effect and the credibility of the eight-miRNA signature are clearly superior

to a single miRNA (Figure 3, Figure S3 and Table S4), implying that the cooperative regulation of multi-miRNAs might play very important roles in tumorigenesis and progression of ccRCC.

In the present work, we have found that eight prognostic miRNAs could interplay with twenty-two TFs. Many studies have revealed that the interplay between miRNA and TFs play key roles in establishing and maintaining cell phenotype [67–69]. Notably, in the gene regulatory network, TF and miRNA interplay could constitute positive or negative feedback loops to execute similar and opposing functions, which can precisely control the regulation of gene expression to reduce noise and maintain cell homeostasis [40,70,71]. Therefore, while considering the complexity of TFs and miRNAs interplay regulating gene expression, we further constructed an interplay network of TF-miRNA and propose a potential molecular mechanism that the interaction between TFs and miRNAs facilitates the survival of ccRCC patients (Figure 10). As shown in Figure 10, down-regulated transcription factors, such as KLF5 and GATA2, can down-regulate miR-183-3p and miR-1269a to up-regulated FAS and other tumor suppressor genes, whilst down-regulated miR-183-3p and miR-1269a can also up-regulate IKZF1 and IRF1 to down-regulated downstream targets. Additionally, up-regulated MYC and IKZF1 can activate miR-223-3p and miR-365b-3p to down-regulate RPS6KA6, DEPTOR, and other tumor suppressor genes, whilst miR-223-3p and miR-365b-3p can also down-regulate TFAP2A and GATA2 to down-regulate miR-183-5p and miR-1269a to up-regulated FAS and VEGFA. Remarkably, some previous reports showed that the transcription factor GATA2 could activate the expression of miR-194 to promote this distant metastasis of prostate cancer by inhibiting SOCS2 [72], and the transcription factor KLF5 could also promote the expression of miR-145, miR-124, and miR-183 by binding to their promoter involved in the progression of invasive pituitary adenoma [73], as well as the transcription factor TFAP2C could promote lung tumorigenesis and aggressiveness through it activating miR-183 and miR-33a-mediated cell cycle regulation [74]. Especially, the interplay of MYC and hsa-miR-144 has been reported in chronic myelogenous leukemia cell K562 [75], and the transcription factor E2F1 could also up-regulate miR-224/452 expressions to inhibit the expression of TXNIP to drive EMT in malignant melanoma [36], as well as miR-3188, as a tumour suppressor, could control the nasopharyngeal carcinoma proliferation and chemosensitivity through a mechanism where FOXO1 modulated a positive feedback loop of mTOR-p-PI3K/AKT-c-JUN [38]. These above results supported our conclusion that the interplays between transcription factors and miRNAs might play very important roles in the prognosis of ccRCC patients. Thus, based on the importance of TF and miRNA interplay in gene expression regulation, we ultimately screened six miRNAs (hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-1269a, hsa-miR-144-5p, hsa-miR-183-5p, hsa-miR-335-3p) and five TFs (TFAP2A, KLF5, IRF1, MYC, IKZF1) as an integrative prognostic predictor. Interestingly, the prognostic effect and the credibility of the combined six miRNAs and five TFs signature are better than both the eight-miRNA signature and the nine-TF signature. The reason might be that not only TFs can regulate the expression of multiple target genes, including miRNAs, but also miRNAs can fine-tune multiple gene expressions, including TFs, as well as their closely coordinated regulations control cell homeostasis. This also suggests why TFs and miRNAs interplay is effective as a clinical prognostic factor for ccRCC patients. Of course, all of these regulatory pairs predicted by bioinformatics and data integration are still to be further verified experimentally.

#### **4. Materials and Method**

#### *4.1. Data Sources and Pre-Processing*

All KIRC sample RNA-seq data of mRNA and miRNA isoform and corresponding clinical information were downloaded from the The Cancer Genome Atlas (TCGA) database (https://portal. gdc.cancer.gov/, Springer Netherlands, Bethesda, MD, USA). The samples were filtered based on survival days greater than three months and simultaneous possessing mRNA and miRNA expression data. Based on the AnimalTFDB3.0 (http://bioinfo.life.hust.edu.cn/AnimalTFDB/, Oxford University Press, Hubei, China.) and Ensemble (http://asia.ensembl.org/index.html, Oxford University Press, Cambridgeshire, UK.) annotations, we identified 19,780 coding genes, of which 1400 were transcription

factors. The expression level of 2,104 miRNA matures was obtained after miRNA isoform alignment. These lower expressed genes (sum(cpm) < 1) were removed and these genes expressed in at least 50% of the sample were retained.

#### *4.2. Di*ff*erentially Expressed Gene Analysis*

Difference gene expression analysis between all tumor and paracancerous tissues was performed while using the edgeR package with filter parameters |log2FC| > 1 and *p*.adjust < 0.05. Similarly, differentially expressed miRNA was screened by limma package with criteria: |log2FC| > 1, *p*.adjust < 0.05. In addition, the mRNA and miRNA expression profiles were further respectively converted to log2 (normalized value + 1) and log2 (RPKM + 1) to be used for the next operation.

#### *4.3. Survival Analysis and Prognosis Model Establishment*

The samples of ccRCC with OS > 90 days were selected for survival analysis. Firstly, according to the median expression of miRNAs, batch survival analysis was performed to screen out the significantly differentially expressed miRNAs that are associated with survival. Subsequently, the univariate Cox regression was used to further assess the miRNAs related to survival. Only those miRNAs with a *p*-value < 0.001 were selected as candidate biomarker miRNAs. Finally, these candidate miRNAs were subjected to multivariate cox regression to determine independent prognostic marker miRNAs and calculate the risk value constructing prediction model for each miRNA. Based on the above results, the time-dependent receiver operating characteristic (ROC) curve was drawn using the Survival ROC R package, and the classification model was evaluated according to the area under the curve (AUC). In addition, we also analyzed the concordance index (C-index) and the Akaike information criterion (AIC). The C-index represents the consistency of the probability of the actual occurrence of the outcome and the probability of the prediction, and the AIC represents a standard for measuring the goodness of statistical model fitting. The prognostic signature of individual calculated according to the risk values of each marker miRNA and combined miRNAs. Next, the patient were divided into high and low risk groups according to the median risk score, and the survival analysis curve was then performed to check significant difference of patients in two groups over time.

#### *4.4. miRNA Target Prediction and Target Function Analysis*

These target genes of prognostic miRNAs were predicted by integrating Mirwalk3.0 (http://mirwalk. umm.uni-heidelberg.de/, Public Library of Science, Mannheim, Germany.) and a negative correlation between miRNA and mRNA expression. The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed using the clusterProfiler package [76] with the filtration standard: *p*.adjust < 0.05. The PPI network from target genes was derived from the STING database (https://string-db.org/, Oxford University Press, Zurich, Switzerland).

#### *4.5. Prediction of Transcription Factors Regulating miRNAs*

The TransmiR2.0 (http://www.cuilab.cn/transmir, Oxford University Press, Beijing, China.) has collected the human regulatory pair of TFs regulating miRNAs (TFs-miRNAs) based on accurate transcriptional start site (TSS) of miRNA and ChIP-seq sequencing as well as experimental validation. Combining the TFs-miRNAs regulatory pair with the expression positive correlation of miRNAs and TFs, and TFs that may regulate these biomarker miRNAs were screened.

#### *4.6. Data Statistics and Visualization*

All data analysis was performed using R software (version 3.5.1, R Core Team, Vienna, Austria). Cytoscape software [77] (version 3.6.1, Cold Spring Harbor Laboratory Press, Washington, WA, USA.) were used to visualize the network. The survival curve was plotted using Kaplan Meier function, and the difference significance was evaluated by log-rank test.

### **5. Conclusions**

In this work, we have identified eight prognostic miRNAs. Among them, seven miRNAs (hsa-miR-223-3p, hsa-miR-365b-3p, hsa-miR-3613-5p hsa-miR-9-5p, hsa-miR-183-5p, hsa-miR-335-3p, hsa-miR-1269a) can serve as potential oncogenes, whereas hsa-miR-144-5p might act as a tumor suppressor gene for ccRCC diagnosis. In addition, the eleven-gene signature (hsa-miR-365b-3p, hsa-miR-223-3p, hsa-miR-1269a, hsa-miR-144-5p, hsa-miR-183-5p, hsa-miR-335-3p, *TFAP2A*, *KLF5*, *IRF1*, *MYC*, *IKZF1*) can sever as an effective prognostic predictor to significantly improve the overall survival of ccRCC patients. Especially, our study has revealed a possible molecular mechanism that TFs and miRNAs interplay can cooperatively regulate the expression of oncogenes, driver genes, and tumor suppressor genes to facilitate the survival of ccRCC patients. Thus, our findings not only provide a new insight into the mechanism that TFs and miRNAs interplay control the tumorigenesis and progression of ccRCC, but also identify several novel diagnostic and prognostic biomarkers as well as potential therapeutic targets that are very crucial for making individualized therapeutic strategies of ccRCC patients.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-6694/11/11/1668/s1, Figure S1: The volcano map of differentially expressed genes and miRNAs. The volcano map of differentially expressed miRNA (A) and mRNA (B). The red point represents up-regulated expression and the green point represents down-regulated expression with statistical significance (*p*-value < 0.05), Figure S2: Kaplan-Meier survival based on the risk score of eight miRNAs. Overall survival curves of high-risk and low-risk groups for hsa-miR-1269a (A), hsa-miR-183-5p (B), hsa-miR-9-5p (C), hsa-miR-335-3p (D), hsa-miR-3613-5p (E), hsa-miR-223-3p (F), hsa-miR-365b-3p (G), hsa-miR-144-5p (H), Figure S3: ROC curves based on the risk score of eight miRNAs. ROC curves for high and low risk of hsa-miR-1269a (A), hsa-miR-183-5p (B), hsa-miR-9-5p (C), hsa-miR-335-3p (D), hsa-miR-3613-5p (E), hsa-miR-223-3p (F), hsa-miR-365b-3p (G), hsa-miR-144-5p (H), Figure S4: The correlation between the expression level of eight prognostic miRNAs and clinical factors. The correlation heat map between clinical parameters and prognostic miRNA expression levels in ccRCC. The red box represents a *p*-value < 0.05; the green box represents a *p*-value > 0.05, Figure S5: The KEGG signaling pathway analysis of up-regulated and down-regulated targets. Signaling pathways of KEGG enrichment of up-regulated targets(A) and down-regulated targets(B), Figure S6: Transcription factors regulating eight miRNAs. Six down-regulated transcription factors down-regulate three down-regulated miRNAs (A); sixteen up-regulate transcription factors up-regulate four up-regulated miRNAs (B). The red represents up-regulation of gene expression and the blue represents down-regulation of gene expression. Table S1: These 37 miRNAs are related to ccRCC survival, Table S2: These top 21 miRNAs are significant associated with ccRCC survival, Table S3: The correlation between prognostic miRNA expression levels and clinical factors, Table S4: The concordance index (C-index) and Akaike information criterion (AIC) of combined miRNAs and single miRNA model, Table S5: These top 30 hub genes from PPI network. Table S6: Targets of has-miR-144-5p. Table S7: TF-miRNA regulation pairs, Table S8. The C-index and AIC of TFs and miRNAs.

**Author Contributions:** F.M. and P.J. designed this study and edited this manuscript; S.Q. and X.S. acquired and analyzed the data; S.Q. and P.J. wrote the original draft; P.J. supervised this project; C.W. collected literatures. All authors read and approved the manuscript.

**Funding:** This research was funded by grants from the National Natural Science Foundation of China (No. 31970477) and the Natural Science Foundation of Jiangsu Province (No. BK20191368).

**Acknowledgments:** We thank our colleagues for their suggestions and criticisms on the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **MiR-378a-3p Is Critical for Burkitt Lymphoma Cell Growth**

**Fubiao Niu <sup>1</sup> , Agnieszka Dzikiewicz-Krawczyk <sup>2</sup> , Jasper Koerts <sup>1</sup> , Debora de Jong <sup>1</sup> , Laura Wijenberg <sup>1</sup> , Margot Fernandez Hernandez <sup>1</sup> , Izabella Slezak-Prochazka <sup>3</sup> , Melanie Winkle <sup>1</sup> , Wierd Kooistra <sup>1</sup> , Tineke van der Sluis <sup>1</sup> , Bea Rutgers <sup>1</sup> , Miente Martijn Terpstra <sup>4</sup> , Klaas Kok <sup>4</sup> , Joost Kluiver <sup>1</sup> and Anke van den Berg 1,\***


Received: 16 October 2020; Accepted: 25 November 2020; Published: 27 November 2020

**Simple Summary:** MicroRNAs (miRNAs) are small RNAs that regulate expression of specific target genes. We observed elevated levels of miR-378a-3p in Burkitt lymphoma (BL) and studied its role in the pathogenesis of BL. Inhibition of miR-378a-3p reduced growth of BL cells, confirming its significance in BL. Identification of BL specific target genes of miR-378a-3p revealed four candidates. For two of them, MNT and IRAK4, miR-378a-dependent regulation was confirmed at the protein level. Overexpression of MNT and IRAK4 in BL cell lines resulted in a similar effect as observed upon miR-378a-3p inhibition, suggesting their involvement in the growth regulatory role of miR-378a-3p.

**Abstract:** MicroRNAs (miRNAs) are small RNA molecules with important gene regulatory roles in normal and pathophysiological cellular processes. Burkitt lymphoma (BL) is an MYC-driven lymphoma of germinal center B (GC-B) cell origin. To gain further knowledge on the role of miRNAs in the pathogenesis of BL, we performed small RNA sequencing in BL cell lines and normal GC-B cells. This revealed 26 miRNAs with significantly different expression levels. For five miRNAs, the differential expression pattern was confirmed in primary BL tissues compared to GC-B cells. MiR-378a-3p was upregulated in BL, and its inhibition reduced the growth of multiple BL cell lines. RNA immunoprecipitation of Argonaute 2 followed by microarray analysis (Ago2-RIP-Chip) upon inhibition and ectopic overexpression of miR-378a-3p revealed 63 and 20 putative miR-378a-3p targets, respectively. Effective targeting by miR-378a-3p was confirmed by luciferase reporter assays for MAX Network Transcriptional Repressor (MNT), Forkhead Box P1 (FOXP1), Interleukin 1 Receptor Associated Kinase 4 (IRAK4), and lncRNA Just Proximal To XIST (JPX), and by Western blot for IRAK4 and MNT. Overexpression of IRAK4 and MNT phenocopied the effect of miR-378a-3p inhibition. In summary, we identified miR-378a-3p as a miRNA with an oncogenic role in BL and identified IRAK4 and MNT as miR-378a-3p target genes that are involved in its growth regulatory role.

**Keywords:** Burkitt lymphoma; miR-378a-3p; cell growth; microRNA

#### **1. Introduction**

Burkitt lymphoma (BL) is one of the fastest growing human tumors with a cell doubling time of about 24 h. BL mainly affects children and young adults but can also occur at a later age [1]. The tumor cells are derived from germinal center B (GC-B) cells and usually carry the hallmark translocation involving *MYC* and the immunoglobulin heavy or light chain loci which results in high expression of MYC [2,3].

MicroRNAs (miRNAs) are a class of short noncoding RNAs of about 22 nucleotides. They modulate gene expression at the post-transcriptional level by translational inhibition or by inducing mRNA degradation [4,5]. MiRNAs regulate a wide range of cellular processes, including cell cycle, proliferation, and apoptosis, and they are important determinants of B-cell development and maturation [6]. A widespread deregulation of miRNA expression has been observed in all B-cell lymphoma subtypes [7].

We and others identified distinct miRNA expression patterns in BL and demonstrated the central role of MYC in regulating miRNA levels [8–12]. Functional studies showed crucial roles for the miR-17~92 cluster, miR-28, miR-150, and miR-155 as either oncogenic or tumor suppressor miRNAs in the pathogenesis of BL [9,13–19]. Nevertheless, the role of most of the deregulated miRNAs in BL remains to be explored.

In this study, we carried out small RNA-sequencing in BL cell lines and normal GC-B cells and subsequently focused on downstream functional experiments for miR-378a-3p. We show for the first time that this miRNA is upregulated in BL and confirmed its regulation by MYC. Further analysis indicated that miR-378a-3p is essential for the growth of BL cells. Using a combination of genome-wide target gene identification, luciferase reporter assays, and Western blot upon modulating miR-378a-3p levels, we confirmed targeting of Interleukin 1 Receptor Associated Kinase 4 (IRAK4) and MAX Network Transcriptional Repressor (MNT) by miR-378a-3p. Overexpression of these two genes phenocopied the effect of miR-378a-3p inhibition on growth of BL cells.

#### **2. Results**

#### *2.1. MiRNA Expression Profiling in BL and GC-B Cells*

An overview of the total number of reads and percentages of mapped reads per sample is given in Table S1. The top 10 most abundantly expressed miRNAs accounted for 73% of all reads in BL and for 71% in GC-B cells (total GC-B cells sorted on the basis of a CD20+IgD−CD38<sup>+</sup> or IgD−CD138−CD3−CD10<sup>+</sup> phenotype). Seven of the top 10 most abundantly expressed miRNAs were shared between BL and GC-B cells (Figure 1A). Twenty-six miRNAs were significantly differentially expressed between BL and GC-B cells, including eight miRNAs upregulated in BL and 18 downregulated (Figure 1B). qRT-PCR validation on the same set of samples confirmed differential expression for six out of eight selected miRNAs (Figure 1C and Figure S1). Of the six validated miRNAs, miR-378a-3p levels were increased in BL relative to GC-B cells, while expression levels of miR-28-5p, miR-155-5p, miR-363-3p, miR-221-3p, and miR-222-3p were decreased. Further expression analysis in primary BL tissue samples and GC-B cells confirmed the differential expression for five of the six miRNAs, excluding miR-221-3p (Figure 1D).

*Cancers* **2020**, *12* 3 of 18

**Figure 1.** Deregulated expression patterns of microRNAs (miRNAs) in Burkitt lymphoma (BL) compared to germinal center B (GC-B) cells. (**A**) Overview of the top 10 most abundantly expressed miRNAs in Burkitt lymphoma (BL) and normal germinal center B (GC-B) cells as determined by small RNA-sequencing. Asterisks indicate miRNAs present in the top 10 of both BL and GC-B cells. (**B**) Heatmap of miRNAs significantly differentially expressed between BL and GC-B cells. (**C**) qRT-PCR validation results for six of the eight tested miRNAs with significantly differential expression between BL cell lines and GC-B cells. MiRNA expression levels were normalized to Small Nucleolar RNA C/D Box 44 (SNORD44). (**D**) The differential expression pattern was confirmed for five of the six tested miRNAs when BL tissues and GC-B cells were compared. MiRNA expression levels were normalized to Small Nucleolar RNA C/D Box 49 (SNORD49). Significant differences were calculated using an unpaired *t*-test. \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001, and \*\*\*\* *p* < 0.0001 **Figure 1.** Deregulated expression patterns of microRNAs (miRNAs) in Burkitt lymphoma (BL) compared to germinal center B (GC-B) cells. (**A**) Overview of the top 10 most abundantly expressed miRNAs in Burkitt lymphoma (BL) and normal germinal center B (GC-B) cells as determined by small RNA-sequencing. Asterisks indicate miRNAs present in the top 10 of both BL and GC-B cells. (**B**) Heatmap of miRNAs significantly differentially expressed between BL and GC-B cells. (**C**) qRT-PCR validation results for six of the eight tested miRNAs with significantly differential expression between BL cell lines and GC-B cells. MiRNA expression levels were normalized to Small Nucleolar RNA C/D Box 44 (SNORD44). (**D**) The differential expression pattern was confirmed for five of the six tested miRNAs when BL tissues and GC-B cells were compared. MiRNA expression levels were normalized to Small Nucleolar RNA C/D Box 49 (SNORD49). Significant differences were calculated using an unpaired *t*-test. \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001, and \*\*\*\* *p* < 0.0001

#### *2.2. MYC-Induced miR-378a-3p Controls BL Cell Growth 2.2. MYC-Induced miR-378a-3p Controls BL Cell Growth*

We selected miR-378a-3p for further functional analysis, because it was the only significantly upregulated miRNA with a high expression level in BL. Previous studies demonstrated that miR-378a-3p is induced by MYC in human mammary epithelial cells [20]. We assessed the regulatory role of MYC in B cells using the P493-6 B-cell model that has a tetracycline-repressible We selected miR-378a-3p for further functional analysis, because it was the only significantly upregulated miRNA with a high expression level in BL. Previous studies demonstrated that miR-378a-3p is induced by MYC in human mammary epithelial cells [20]. We assessed the regulatory role of MYC in B cells using the P493-6 B-cell model that has a tetracycline-repressible MYC allele [21]. Our results showed that this miRNA is also induced by MYC in B cells (Figure 2A).

MYC allele [21]. Our results showed that this miRNA is also induced by MYC in B cells (Figure 2A). To explore the role of miR-378a-3p in growth of BL cells, we inhibited miR-378a-3p using a lentiviral miRNA inhibition construct (mZip-378a-3p) in four BL cell lines and followed cell growth To explore the role of miR-378a-3p in growth of BL cells, we inhibited miR-378a-3p using a lentiviral miRNA inhibition construct (mZip-378a-3p) in four BL cell lines and followed cell growth in a GFP

competition assay. Compared to the negative control (mZip-SCR), a significant decrease in the number of GFP-positive cells was observed in three (ST486, CA46, and DG75) out of four BL cell lines (Figure 2B). No relationship was observed between the reduction in the percentage of GFP<sup>+</sup> cells and the level of miR-378a-3p expression (Figure S2). Together, our data indicate that inhibition of miR-378a-3p is disadvantageous for BL cells, suggesting miR-378a-3p is indispensable for growth of BL cells.

Overexpression of miR-378a in ST486 using a lentiviral miRNA overexpression construct (pCDH-378a) resulted in a ~47-fold increase in miR-378a-3p level. In a GFP competition assay, miR-378a overexpression had no effects on cell growth, probably due to the already high endogenous levels.

**Figure 2.** MYC-induced miR-378a-3p is essential for BL cell growth. (**A**) Levels of MYC and miR-378a-3p in tetracycline treated ("−", MYC-off) and non-treated ("+", MYC-on) P493-6 B-cells. MYC levels were normalized to RNA Polymerase II Subunit A (POLR2A). MiRNA levels were normalized to SNORD44. (**B**) Green fluorescent protein (GFP) growth competition assay upon miR-378a-3p inhibition in four Burkitt lymphoma (BL) cell lines. The miR-378a inhibitor (mZip-378a-3p) and the scrambled control (mZip-SCR) were stably transduced in BL cells using a lentiviral vector, which co-expresses GFP. The GFP percentage was measured for 18 days, and the GFP percentage at the first day of measurement (day 4) was set to 100%. All assays were performed in triplicate. Significant differences were calculated using a mixed model analysis. \*\*\* *p* < 0.001 and \*\*\*\* *p* < 0.0001.

#### *2.3. Identification of miR-378a-3p Targets*

To identify miR-378a-3p target genes, we performed Ago2-RIP-Chip upon miR-378a-3p inhibition and overexpression in ST486 cells (Figure S3A,B). Efficient pulldown of Ago2-containing RISC and miRNAs was confirmed by qRT-PCR for miR-378a-3p and the unrelated highly expressed miR-181a-5p (Figure S3C,D) and by Western blot for Ago2 protein (Figures S3E and S4). The number of Ago2 immunoprecipitation (IP)-enriched probes was similar in all four conditions, ranging between 6.3% and 9.8% of the probes with consistent expression levels (Table 1).

**Table 1.** Number of genes in the miRNA targetome upon miR-378a-3p overexpression (pCDH) and inhibition (mZip).


EV = pCDH-EV, 378a = pCDH-378a, SCR = mZip-SCR, IP = Ago2 immunoprecipitated fraction, T = total fraction.

A total of 22 probes corresponding to 20 genes showed a ≥2-fold increased IP enrichment upon miR-378a overexpression compared to empty vector control infected cells (Figure 3A and Table 2). Nine of the 20 genes (45%) had at least one putative miR-378a-3p binding site (7mer-A1, 7mer-m8, or/and 8mer). Six of them, i.e., MAX Network Transcriptional Repressor (MNT), Heat Shock Protein Family B (Small) Member 1 (HSPB1), Interleukin 1 Receptor Associated Kinase 4 (IRAK4), Cyclin K (CCNK), Cyclin Dependent Kinase Inhibitor 2A (CDKN2A), and Ring Finger Protein 34 (RNF34), had a Gene Ontology (GO) term related to cell growth, apoptosis, and/or cell cycle. Upon miR-378a-3p inhibition, 74 probes, corresponding to 63 genes, showed a ≥2-fold decreased IP enrichment compared to the negative control infected cells (Figure 3B and Table 3). Nineteen of these 63 genes (30%) contained at least one putative miR-378a-3p binding site, including Cytokine Inducible SH2 Containing Protein (CISH), BCR Activator Of RhoGEF And GTPase (BCR), Tubulin Alpha 1c (TUBA1C), SWI/SNF Related, Matrix Associated, Actin Dependent Regulator Of Chromatin, Subfamily A, Member 4 (SMARCA4), and Forkhead Box P1 (FOXP1) with a GO term related to cell growth, apoptosis, or cell cycle. One of the target genes, i.e., MYC Binding Protein (MYCBP), was identified with both experimental set-ups.



\* IP/T ratios in pCDH-EV (EV) were set to 1.0 in cases where the ratios were <1.0. \*\* The binding site in the noncoding RNA is listed in the 30 -UTR column. 7mA1 = 7mer-A1, 7m8 = 7mer-m8, 8m = 8mer, FC = fold change, 378a = pCDH-378a, IP = Ago2 immunoprecipitated fraction, T = Total fraction, 50UTR = 5 0 untranslated region, CDS = coding sequence, 30UTR = 3 0 untranslated region, GO = gene ontology.

**Table 3.** Identified targets of miR-378a-3p upon inhibition.



**Table 3.** *Cont*.

\* IP/T ratios in mZip-378a-3p were set to 1.0 if <1.0. \*\* Binding sites on noncoding RNAs are listed in 30 -UTR column. 7mA1 = 7mer-A1, 7m8 = 7mer-m8, 8m = 8mer, FC = fold change, SCR = mZip-SCR, N/A = not available, IP = Ago2 immunoprecipitated fraction, T = total fraction, 50UTR = 5 0 untranslated region, CDS = coding sequence, 30UTR = 3 0 untranslated region, GO = gene ontology.

#### *2.4. IRAK4 and MNT Are Involved in the Function of miR-378a-3p*

For further analysis, we selected MYCBP that was identified in both approaches and six candidates that had at least one 7mer-A1, 7mer-m8, or an 8mer and a GO term related to cell growth, apoptosis, or/and cell cycle (CISH, BCR, TUBA1C, FOXP1, MNT, and IRAK4). We also included the lncRNA JPX, as it showed a strong enrichment upon miR-378a-3p overexpression and contained an 8-mer seed binding site (Figure 3C).

*Cancers* **2020**, *12* 7 of 18

**Figure 3.** Identification and validation of miR-378a-3p targets. MiR-378a-3p targets identified by Ago2-RIP-Chip upon (**A**) miR-378a overexpression relative to empty vector (EV) and (**B**) scrambled vector relative to miR-378a-3p inhibition (SCR/mZip-378a). The black bars indicate genes with miR-378a-3p seed binding sites. (**C**) Schematic representations of the eight genes selected for luciferase reporter assay validation. Black boxes indicate positions of the open reading frames (ORFs). Positions and types of miR-378a-3p binding sites are indicated relative to the ORF. The binding site in MNT indicated by an asterisk was not tested. (**D**) Luciferase reporter assay results upon co-transfection of ST486 and DG75 cells with the Psi-check-2 construct containing the wildtype (WT) or mutated (MUT) miR-378a-3p binding sites from the selected genes and either an miR-378a-3p mimic or a negative control mimic. Significant differences were calculated using a paired *t*-test. \* *p* < 0.05, \*\* *p* < 0.01, ns = not significant. **Figure 3.** Identification and validation of miR-378a-3p targets. MiR-378a-3p targets identified by Ago2-RIP-Chip upon (**A**) miR-378a overexpression relative to empty vector (EV) and (**B**) scrambled vector relative to miR-378a-3p inhibition (SCR/mZip-378a). The black bars indicate genes with miR-378a-3p seed binding sites. (**C**) Schematic representations of the eight genes selected for luciferase reporter assay validation. Black boxes indicate positions of the open reading frames (ORFs). Positions and types of miR-378a-3p binding sites are indicated relative to the ORF. The binding site in MNT indicated by an asterisk was not tested. (**D**) Luciferase reporter assay results upon co-transfection of ST486 and DG75 cells with the Psi-check-2 construct containing the wildtype (WT) or mutated (MUT) miR-378a-3p binding sites from the selected genes and either an miR-378a-3p mimic or a negative control mimic. Significant differences were calculated using a paired *t*-test. \* *p* < 0.05, \*\* *p* < 0.01, ns = not significant.

To confirm targeting of the eight selected genes by miR-378a-3p, we carried out luciferase reporter assays for 10 putative miR-378a-3p binding sites in ST486 and DG75. This revealed a strong reduction in the *Renilla*/firefly ratio for four of the miR-378a-3p binding sites in four genes (IRAK4, FOXP1 (site

8-mer seed binding site (Figure 3C).

1), MNT, and JPX) (Figure S5). To further confirm specific binding by miR-378a-3p, we generated constructs with mutations in these four miR-378a-3p binding sites. For IRAK4 (trend), FOXP1, and MNT (both significant), wildtype binding sites showed lower *Renilla*/firefly ratios compared to the mutated binding sites (black bars in Figure 3D), indicating binding of endogenous miR-378a-3p to these sequences. Upon miR-378a-3p overexpression a significantly reduced *Renilla*/firefly ratio was observed for the wildtype, but not for the mutated target sites, confirming efficient and specific targeting (Figure 3D). efficient and specific targeting (Figure 3D). To validate the regulatory role of miR-378a-3p on endogenous IRAK4, FOXP1, and MNT protein levels, we analyzed protein levels upon modulation of miR-378a-3p levels in ST486 and DG75 cells (Figure 4 and Figure S6). For IRAK4, this revealed a significant decrease upon miR-378a-3p overexpression and a significant increase upon miR-378a-3p inhibition in DG75. In ST486, the same pattern was observed albeit not significant. MNT levels were significantly decreased upon miR-378a-3p overexpression in DG75 and ST486. Increased MNT levels upon miR-378a-3p inhibition were observed only in DG75 cells although the increase was not significant.

*Renilla*/firefly ratio was observed for the wildtype, but not for the mutated target sites, confirming

*Cancers* **2020**, *12* 8 of 18

For further analysis, we selected MYCBP that was identified in both approaches and six candidates that had at least one 7mer-A1, 7mer-m8, or an 8mer and a GO term related to cell growth, apoptosis, or/and cell cycle (CISH, BCR, TUBA1C, FOXP1, MNT, and IRAK4). We also included the lncRNA JPX, as it showed a strong enrichment upon miR-378a-3p overexpression and contained an

To confirm targeting of the eight selected genes by miR-378a-3p, we carried out luciferase reporter assays for 10 putative miR-378a-3p binding sites in ST486 and DG75. This revealed a strong reduction in the *Renilla*/firefly ratio for four of the miR-378a-3p binding sites in four genes (IRAK4, FOXP1 (site 1), MNT, and JPX) (Figure S5). To further confirm specific binding by miR-378a-3p, we generated constructs with mutations in these four miR-378a-3p binding sites. For IRAK4 (trend), FOXP1, and MNT (both significant), wildtype binding sites showed lower *Renilla*/firefly ratios compared to the mutated binding sites (black bars in Figure 3D), indicating binding of endogenous

*2.4. IRAK4 and MNT Are Involved in the Function of miR-378a-3p* 

To validate the regulatory role of miR-378a-3p on endogenous IRAK4, FOXP1, and MNT protein levels, we analyzed protein levels upon modulation of miR-378a-3p levels in ST486 and DG75 cells (Figure 4 and Figure S6). For IRAK4, this revealed a significant decrease upon miR-378a-3p overexpression and a significant increase upon miR-378a-3p inhibition in DG75. In ST486, the same pattern was observed albeit not significant. MNT levels were significantly decreased upon miR-378a-3p overexpression in DG75 and ST486. Increased MNT levels upon miR-378a-3p inhibition were observed only in DG75 cells although the increase was not significant. No significant differences in FOXP1 expression were observed in either cell line. No significant differences in FOXP1 expression were observed in either cell line. Next, we performed gain-of-function experiments to determine whether overexpression of the validated targets IRAK4 and MNT could phenocopy the effect of miR-378a-3p inhibition. Overexpression of IRAK4 in DG75 cells resulted in a 74% decrease in GFP+ cells in 11 days (Figure 5). Overexpression of MNT resulted in a >90% decrease in GFP+ cells in 11 days in DG75. Thus, our results indicate that overexpression of IRAK4 and MNT could phenocopy the effect of miR-378a-3p inhibition on growth of BL cells. Altogether, our results suggest that growth of BL cells depends on high miR-378a-3p levels through regulating IRAK4 and MNT.

**Figure 4.** Analysis of the effect of modulating miR-378a-3p levels upon the protein levels of the target genes. Representative examples of the effect of miR-378a overexpression on IRAK4 (**A**), FOXP1 (**C**), **Figure 4.** Analysis of the effect of modulating miR-378a-3p levels upon the protein levels of the target genes. Representative examples of the effect of miR-378a overexpression on IRAK4 (**A**), FOXP1 (**C**), and MNT (**E**). Representative examples of the effect of miR-378a-3p inhibition on IRAK4 (**B**), FOXP1 (**D**), and MNT (**F**). cells. Graphs show the quantification of the protein levels of three to four independent infections. For each protein, the part of the total protein lane corresponding to the position of the protein band is shown as an indication for protein loading. Protein levels were quantified relative to the total protein amount as measured in the complete lane. Uncropped blots can be found in Figure S6. Significant differences were calculated using a paired *t*-test. \* *p* < 0.05, \*\* *p* < 0.01, ns = not significant.

Next, we performed gain-of-function experiments to determine whether overexpression of the validated targets IRAK4 and MNT could phenocopy the effect of miR-378a-3p inhibition. Overexpression of IRAK4 in DG75 cells resulted in a 74% decrease in GFP<sup>+</sup> cells in 11 days (Figure 5). Overexpression of MNT resulted in a >90% decrease in GFP<sup>+</sup> cells in 11 days in DG75. Thus, our results indicate that overexpression of IRAK4 and MNT could phenocopy the effect of miR-378a-3p inhibition on growth of BL cells. Altogether, our results suggest that growth of BL cells depends on high miR-378a-3p levels through regulating IRAK4 and MNT.

< 0.01, ns = not significant.

and MNT (**E**). Representative examples of the effect of miR-378a-3p inhibition on IRAK4 (**B**), FOXP1 (**D**), and MNT (**F**). cells. Graphs show the quantification of the protein levels of three to four independent infections. For each protein, the part of the total protein lane corresponding to the position of the protein band is shown as an indication for protein loading. Protein levels were quantified relative to the total protein amount as measured in the complete lane. Uncropped blots

**Figure 5.** Overexpression of IRAK4 and MNT phenocopies the effect of miR-378a-3p inhibition in DG75 cells. (**A**) Validation of the IRAK4 and MNT overexpression in DG75. An empty vector (EV) was used as a control. For each protein, the part corresponding to the position of the protein band is shown as control of protein loading. Protein levels were quantified relative to the total protein amount in the complete lanes. Uncropped blots can be found in Figure S7. (**B**) GFP competition assays upon overexpression of IRAK4 and MNT resulted in a strong decrease in GFP+ cells over time, while no or only a mild effect was observed for the EV control. The GFP percentages were measured for 11 days, and the GFP percentage as measured on day 4 was set to 100%. Assays were performed in duplicate. Significant differences between IRAK4 and MNT overexpression and EV control were calculated using a mixed model analysis. \*\* *p* < 0.01 and \*\*\*\* *p* < 0.0001. **Figure 5.** Overexpression of IRAK4 and MNT phenocopies the effect of miR-378a-3p inhibition in DG75 cells. (**A**) Validation of the IRAK4 and MNT overexpression in DG75. An empty vector (EV) was used as a control. For each protein, the part corresponding to the position of the protein band is shown as control of protein loading. Protein levels were quantified relative to the total protein amount in the complete lanes. Uncropped blots can be found in Figure S7. (**B**) GFP competition assays upon overexpression of IRAK4 and MNT resulted in a strong decrease in GFP<sup>+</sup> cells over time, while no or only a mild effect was observed for the EV control. The GFP percentages were measured for 11 days, and the GFP percentage as measured on day 4 was set to 100%. Assays were performed in duplicate. Significant differences between IRAK4 and MNT overexpression and EV control were calculated using a mixed model analysis. \*\* *p* < 0.01 and \*\*\*\* *p* < 0.0001.

#### **3. Discussion 3. Discussion**

In this study, we identified 26 miRNAs differentially expressed in BL cell lines compared to GC-B cells. For five of the miRNAs, deregulated expression levels were confirmed in both BL cell lines and primary tissues. Among them, miR-378a-3p is MYC-induced, highly abundant (top 10 within BL), and overexpressed in BL compared to GC-B cells. Inhibition of miR-378a-3p showed a negative effect on BL cell growth. In a genome-wide Ago2-RIP-Chip analysis, 20 and 63 genes were identified as the potential targets of miR-378a-3p upon miR-378a-3p overexpression and inhibition, respectively. MNT and IRAK4 were confirmed as novel targets of miR-378a-3p in BL, and their In this study, we identified 26 miRNAs differentially expressed in BL cell lines compared to GC-B cells. For five of the miRNAs, deregulated expression levels were confirmed in both BL cell lines and primary tissues. Among them, miR-378a-3p is MYC-induced, highly abundant (top 10 within BL), and overexpressed in BL compared to GC-B cells. Inhibition of miR-378a-3p showed a negative effect on BL cell growth. In a genome-wide Ago2-RIP-Chip analysis, 20 and 63 genes were identified as the potential targets of miR-378a-3p upon miR-378a-3p overexpression and inhibition, respectively. MNT and IRAK4 were confirmed as novel targets of miR-378a-3p in BL, and their overexpression phenocopied the effect of miR-378a-3p on BL cell growth.

overexpression phenocopied the effect of miR-378a-3p on BL cell growth. Six of the 26 identified differentially expressed miRNAs were reported to be deregulated in BL by Oduor et al. in a previous study [8]. These six included miR-155-5p, miR-221-3p, miR-222-3p, and miR-28-5p which we validated by qRT-PCR in BL cell lines and tissue samples. Nine additional miRNAs were proven to be differentially expressed in BL compared to other B-cell lymphomas [9– 12,20]. Thus, our small RNA sequencing data confirmed some of the previously identified deregulated miRNAs in BL. Moreover, we showed for the first time that miR-378a-3p is an MYC-induced and significantly upregulated miRNA in BL. Since the role of miR-378a-3p in BL was Six of the 26 identified differentially expressed miRNAs were reported to be deregulated in BL by Oduor et al. in a previous study [8]. These six included miR-155-5p, miR-221-3p, miR-222-3p, and miR-28-5p which we validated by qRT-PCR in BL cell lines and tissue samples. Nine additional miRNAs were proven to be differentially expressed in BL compared to other B-cell lymphomas [9–12,20]. Thus, our small RNA sequencing data confirmed some of the previously identified deregulated miRNAs in BL. Moreover, we showed for the first time that miR-378a-3p is an MYC-induced and significantly upregulated miRNA in BL. Since the role of miR-378a-3p in BL was not studied before, we focused on this miRNA for further functional analysis.

not studied before, we focused on this miRNA for further functional analysis. Inhibition of miR-378a-3p resulted in a strong reduction of BL cell growth, suggesting a possible oncogenic role of miR-378a-3p in BL. The effect on growth upon miR-378a-3p was most pronounced in ST486. DG75 and CA46 showed intermediate phenotypes while no significant effect was observed in Ramos. There was no obvious relationship between endogenous miR-378a-3p levels and the decrease in percentage of GFP+ cells upon miR-378a-3p inhibition. Given the complex interactions between miRNAs and target genes, i.e., multiple targets per miRNA and multiple miRNAs per Inhibition of miR-378a-3p resulted in a strong reduction of BL cell growth, suggesting a possible oncogenic role of miR-378a-3p in BL. The effect on growth upon miR-378a-3p was most pronounced in ST486. DG75 and CA46 showed intermediate phenotypes while no significant effect was observed in Ramos. There was no obvious relationship between endogenous miR-378a-3p levels and the decrease in percentage of GFP<sup>+</sup> cells upon miR-378a-3p inhibition. Given the complex interactions between miRNAs and target genes, i.e., multiple targets per miRNA and multiple miRNAs per target, the differences in the observed phenotypes might be related to endogenous levels of other miRNAs or target genes. Previous studies have shown opposite roles of miR-378a-3p in different cancer types. MiR-378a-3p was reported to inhibit growth or promote apoptosis and, thus, act as a tumor suppressor in colorectal cancer, lung cancer, ovarian cancer, prostate cancer, and rhabdomyosarcoma [22–26]. In contrast, in line with our findings, miR-378a-3p was shown to promote proliferation and reduce apoptosis in gastric cancer, nasopharyngeal carcinoma, colorectal cancer, and acute myeloid leukemia [27–30].

Using an unbiased genome-wide experimental approach, we identified 83 putative miR-378a-3p target genes. Luciferase reporter assays confirmed targeting of four out of eight selected genes, i.e., the protein-coding genes MNT, IRAK4, and FOXP1, and the lncRNA JPX. Targeting of the endogenous *MNT* and *IRAK4* transcripts by miR-378a-3p was confirmed at the protein level, albeit with limited effects. For FOXP1, we could not confirm targeting by miR-378a-3p at the protein level, while we did not follow up potential targeting of endogenous JPX transcripts in BL. Moreover, we showed that overexpression of MNT and IRAK4 strongly inhibited the growth of DG75 cells. Together, these data suggest that the effect of miR-378a-3p might at least in part be dependent on targeting MNT and IRAK4. The apparent discrepancy between the limited effect of miR-378a-3p on the levels of these proteins and the relatively strong phenotype on growth is in line with the general thought that miRNAs most often do not work as on/off switches but rather fine-tune the expression of their targets [5]. A combined moderate effect on MNT, IRAK4, and possible others might explain the strong effect on cell growth. Further experiments in additional cell lines using, e.g., phenotype rescue approaches with more precise control of the overexpression levels are required to confirm the relevance of these targets for the miR-378a-3p-induced phenotype.

Previous studies showed a dual role of MNT in tumorigenesis. On the one hand, MNT was reported as a facilitator of MYC-driven T-cell proliferation and survival [31]. In line with this, a study in Eµ-MYC mice showed that reduced MNT levels reduced tumorigenesis [32], suggesting that MNT is indispensable for MYC-driven oncogenesis. However, in most studies, MNT acted as a tumor suppressor and was a functional antagonist of MYC by repressing its activities related to cell cycle, proliferation, and apoptosis [33,34]. Loss of MNT in mouse embryonic fibroblasts (MEFs) phenocopied the effect of MYC overexpression [35–37]. These findings are in line with our results and suggest that high levels of miR-378a-3p could promote BL tumorigenesis by reducing MNT levels, thereby enabling MYC to execute its oncogenic effects in BL.

IRAK4 plays an essential role in the Toll-like receptor (TLR) pathway [38], which mediates inflammatory signals in B cells and causes activation of nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB). The TLR pathway is hyperactive in mantle cell lymphoma and diffuse large B-cell lymphoma (DLBCL), and activation of NF-κB promotes B-cell survival and proliferation [39–41]. Depletion of IRAK4 showed a negative effect on NF-κB activity and autocrine IL-6/IL-10 engagement of the Janus Kinase (JAK)-Signal Transducer and Activator of Transcription 3 (STAT3) pathway, reducing survival of DLBCL cells [42,43]. Despite the pro-survival role of NF-κB in DLBCL, activation of NF-κB has been reported to be disadvantageous in MYC positive BL consistent with our data, supporting a role of miR-378a-3p-dependent repression of IRAK4 in limiting the activation of NF-κB [44,45].

Although we could not show that FOXP1 protein levels are affected by miR-378a-3p, we cannot exclude that expression changes in FOXP1 are more subtle and not captured by our experimental set-up. FOXP1 is a member of the forkhead box (Fox) transcription factor family and a regulator of early B-cell development [46]. Interestingly FOXP1 expression is downregulated during the normal GC reaction. In BL, FOXP1 levels are comparable to GC-B cells, but lower than the level in other B-cell lymphomas [47,48]. The potential relevance of maintaining low FOXP1 in BL is further supported by the finding that aberrant expression of FOXP1 cooperates with (constitutive) NF-κB activity [49], which might be disadvantageous for the survival of BL. The relevance of the long noncoding RNA JPX in the phenotype induced by miR-378a-3p inhibition remains to be elucidated. JPX is an activator of X Inactive Specific Transcript (XIST) and acts as a molecular switch for X chromosome inactivation [50]. JPX was reported to act as an oncogene in ovarian cancer and non-small-cell lung cancer by promoting cell proliferation, invasion, and migration [51,52]. In hepatocellular carcinoma, JPX-dependent induction of XIST suppresses hepatocellular carcinoma progression by binding to the miR-155-5p oncomiR [53]. The potential role of miR-378a-3p in regulating FOXP1 and JPX levels need to be further established.

#### **4. Materials and Methods**

#### *4.1. BL Cell Lines, Germinal Center (GC) B Cells, and BL Patient Material*

BL cell lines were purchased from American Type Culture Collection (ATCC/LGC Standards, Molsheim Cedex, France) (ST486 and Ramos) and German Collection of Microorganisms and Cell Cultures (DSMZ, Braunschweig, Germany) (CA46 and DG75). BL cells were cultured at 37 ◦C under an atmosphere containing 5% CO<sup>2</sup> in Roswell Park Memorial Institute (RPMI)-1640 medium (Cambrex Biosciences, Walkersville, MD, USA) supplemented with 2 nM ultra-glutamine, 100 U/mL penicillin, 0.1 mg/mL streptomycin, and 10% (CA46, DG75, and Ramos) or 20% (ST486) fetal bovine serum (Sigma-Aldrich, Zwijndrecht, The Netherlands). P493-6 B-cells were cultured as described previously [54]. We routinely confirmed cell line identity using the PowerPlex® 16HS System (Promega, Leiden, The Netherlands) and absence of mycoplasma contamination.

GC-B cells and frozen BL tissue sections were obtained previously as described in [9,54,55]. GC-B cells (defined as CD20+IgD−CD38+, *n* = 6 or IgD−CD138−CD3−CD10+, *n* = 1) were sorted from routinely removed tonsil specimens of children. Specifically, for small RNA seq experiments, fluorescence-activated cell sorter (FACS)-sorted (CD20+IgD−CD38+, *n* = 2) and magnetic-activated cell sorter MACS-sorted (IgD−CD138−CD3−CD10+, *n* = 1) GC-B cells were used as controls, while FACS-sorted (CD20+IgD−CD38+, *n* = 4) GC-B cells were used for qRT-PCR validation experiments. Written permission for the use of the tonsil tissues to isolate GC-B cells was obtained from the parents of the children. The study protocol was consistent with international ethical guidelines (the Declaration of Helsinki and the International Conference on Harmonization Guidelines for Good Clinical Practice). According to the Medical ethics review board of the University Medical Center Groningen our studies fulfilled requirements for patient anonymity and were in accordance with their regulations. The Medical ethics review board waives the need for approval if rest material is used under law in the Netherlands, and waives the need for informed consent when patient anonymity is assured (BL tissue samples).

#### *4.2. RNA Isolation*

RNA was isolated using miRNeasy Mini or Micro kit (Qiagen, Hiden, Germany) according to the manufacturer's instructions. RNA concentration was measured by a NanoDropTM 1000 Spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA) and integrity was evaluated on a 1% agarose gel.

#### *4.3. Small RNA Library Preparation and Sequencing*

Firstly, 1–2 µg of total RNA from four BL cell lines and three samples of GC-B cells was used to generate small RNA libraries using Truseq Small RNA Sample Preparation Kit and TruSeq small RNA indices (Illumina, San Diego, CA, USA). Sequencing was performed on an Illumina 2000 HiSeq platform. After removal of 30 - and 50 -adaptor sequences from the raw reads using the CLC Genomics Workbench (CLC Bio, Cambridge, MA, USA), sequencing data were analyzed using miRDeep version 2.0 (Max Delbrück Center for Molecular Medicine in the Helmholtz Association, https://www.mdcberlin.de/8551903/en) [56] and annotated against miRbase version 21 (http://www.mirbase.org) [57] allowing one nucleotide mismatch. Read counts of miRNAs with the same mature miRNA sequence were merged. Total read counts per sample were normalized to 1,000,000. For statistical analysis, we included all unique miRNAs with at least 50 read counts in all seven samples. Genesis software v1.7.6 (Institute for Genomics and Bioinformatic Graz, Graz, Austria) was used to generate the heat map. The small RNA sequencing data were deposited in the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo; accession number GSE92616).

#### *4.4. qRT-PCR*

For validation of small RNA sequencing results, we selected differentially expressed miRNAs with expression log<sup>2</sup> reads per million (RPM) >8 in BL or GC-B cells. We selected the most optimal Taqman assay on the basis of isoform abundance as observed in the small RNA sequencing data (Table S2). The miRNA expression levels were analyzed using Taqman miRNA quantitative PCR assays (Thermo Fisher Scientific Inc.) in a multiplexed fashion as described previously [58]. Cycle crossing point (Cp) values were determined with Light Cycler 480 software version 1.5.0 (Roche, Basel, Switzerland). Relative expression levels of miRNAs to the housekeeping gene (SNORD44 or SNORD49) were determined by calculating 2−∆Cp (∆Cp <sup>=</sup> CpmiRNA <sup>−</sup> Cphouse-keeping gene). MYC transcript levels were analyzed as indicated previously [54].

#### *4.5. Lentiviral Constructs, Transduction, and GFP Competition Assay*

Lentiviral constructs to inhibit (mZip-378a-3p) or overexpress (pCDH-miR-378a) miR-378a-3p were purchased from System Biosciences (Palo Alto, CA, USA). A nontargeting mZip-scrambled (SCR) and an empty vector pCDH (EV) construct were used as negative controls. Lentiviral pLV-EGFP:T2A:Puro-EF1A IRAK4, MNT, and empty vector control constructs were purchased from VectorBuilder (Chicago, IL, USA). Lentiviral particles were produced in HEK-293T cells by calcium phosphate precipitation transfection using a third-generation packaging system as described previously [55]. Briefly, HEK-293T cells were seeded in six-well plates and grown until ~80% confluence. A plasmid mix consisting of 15 µL of CaCl<sup>2</sup> (2.5 M), 1 µg of pMSCV-VSV-G, 1 µg of pRSV.REV, 1 µg of pMDL-gPRRE, 2 µg of lentiviral vector, and 150 µL of 2× HEPES buffered saline (HBS) was prepared to transfect the HEK-293T cells. The virus was harvested and filtered using a 0.45 µm filter 48 h after transfection. The virus was either used directly or stored at −80 ◦C.

For GFP competition assays, BL cell lines were infected with the mZip-378a-3p and the negative control mZip-SCR in three biological replicates per construct, aiming at an infection efficiency of 20% to 50% GFP<sup>+</sup> cells on day 4. The percentage of GFP<sup>+</sup> cells was monitored by flow cytometry (BD Biosciences, San Jose, CA, USA) three times per week for a total period of 18 days. For the IRAK4 and MNT GFP competition assays, two biological replicates per construct were performed, and the percentage of GFP<sup>+</sup> cells was monitored for 11 days.

#### *4.6. AGO2-IP Procedure*

Immunoprecipitation (IP) of the Ago2-containing RNA-induced silencing complex (RISC) (Ago2-IP) procedure was done as described previously [59]. To identify miR-378a-3p target genes, we applied the Ago2-RIP-Chip on BL cells infected with lentiviral miR-378a-3p inhibition or overexpression constructs. For both constructs, a parallel infection with appropriate control constructs (nontargeting or empty vector) was performed. We aimed at a high infection percentage and harvested the cells at day 5, either directly or after sorting to reach a GFP<sup>+</sup> percentage >95% for inhibition and >85% for overexpression. For each AGO2-IP experiment, we started with ~30 million cells. RNA was isolated from the Ago2-IP and total (T) fractions. Efficiency of the Ago2-IP procedure was confirmed by qRT-PCR for miR-378a-3p and miR-181a and by Western blot for the Ago2 protein.

#### *4.7. Western Blotting*

Infected cells were harvested and lysed in lysis buffer (#9803, Cell Signaling Technology, Danvers, MA, USA) supplemented with protease inhibitor. After centrifugation at 14,000 rpm for 10 min (4 ◦C), supernatants were collected, and protein concentrations were measured using the bicinchoninic acid (BCA) Protein Assay Kit (Thermo Fisher Scientific Inc., Waltham, MA, USA) according to the manufacturer's protocol. Then, 20 µg of protein was separated on a polyacrylamide gel and transferred to a nitrocellulose membrane followed by incubation overnight at 4 ◦C with primary antibodies, anti-Ago2 (2E12-1C9, Abnova, Taipei, Taiwan), anti-MNT (#A303-626A, Bethyl Lab, Montgomery, TX, USA), anti-IRAK4 (ab32511, Cambridge, UK), and anti-FOXP1 (#2005, Cell Signaling Technology, Danvers, MA, USA)). All antibodies were diluted 1000× in 5% milk in Tris-buffered saline with Tween-20 (TBST). After incubation with the secondary and tertiary (for MNT, IRAK4, and FOXP1) antibodies and with the enhanced chemiluminescence (ECL) substrate (Thermo Fisher Scientific

Inc.), chemiluminescence was detected. Ago2 was quantified with Image Lab 4.0.1 software (BioRad, Hercules, CA, USA), and MNT, IRAK4, and FOXP1 were quantified using ImageJ (NIH, Bethesda, MD, USA). MNT, IRAK4, and FOXP1 protein levels were normalized relative to the total protein amount in the complete lane.

#### *4.8. Microarray Analysis*

About 50 ng RNA of both the total (T) and the IP fractions were labeled and hybridized on an Agilent gene expression microarray (AMADID no.: 072363, SurePrint G3 Human Gene Exp v3 array kit, Agilent Technologies, Santa Clara, CA, USA). The microarray contained 58,341 probes against coding and noncoding transcripts. The procedure and data analysis were performed as previously described [54]. Briefly, after complementary RNA (cRNA) synthesis and amplification, labeling was done with cyanine 3-CTP (Cy3) or cyanine 5-CTP (Cy5) using the LowInput QuickAmp Labeling kit (catalog no.: 0006322867). Equal amounts of Cy3- or Cy5-labeled cRNA samples were mixed and hybridized on the microarray slide overnight. Raw data were quantile normalized without baseline transformation using GeneSpring GX 12.5 software (Agilent Technologies). Probes were selected for further analysis if they were flagged present in all samples, expressed in the 25th to 100th percentile in at least half of the total (T) fractions, and showed consistent expression in the duplicate measurements (<2-fold change). The average signals of replicates were used to calculate the IP/T ratio, and probes with a ≥ 2-fold enrichment in the IP fraction as compared to total (T) fraction were considered as potential miRNA targets.

For pCDH-378a transduced cells, we next assessed miR-378a-3p targets by identifying probes that were enriched at least ≥2-fold higher in miR-378a-overexpressing cells as compared to empty vector (pCDH-EV). For mZip-378a-3p transduced cells, we assessed miR-378a-3p targets by identifying probes with at least ≥2-fold higher enrichment in mZip-SCR transduced cells as compared to mZip-378a-3p cells. Gene expression microarray data are deposited in the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo; accession number: GSE141691).

#### *4.9. Identification of miR-378a-3p Seed Sites and Gene Ontology (GO) Terms*

For all Ago2-RIP-Chip identified targets of miR-378a-3p, we used a Perl script to search for 7mer-A1, 7mer-m8, and 8mer [60] miR-378a-3p seed sites in 50 -untranslated region (UTR), coding sequence (CDS), and 30 -UTR). Ensemble transcript isoforms were selected on the basis of a refseq identifier (ID) conversion using Biomart (https://www.ensembl.org/biomart). For lncRNA transcripts without an Ensembl ID, we used the LNCipedia or a XLOC/TCONS (BROAD Institute, Cambridge, MA, USA) transcript ID. The growth-related Gene Ontology (GO) terms (proliferation, cell cycle, and apoptosis) for selected genes were retrieved from the Ensembl database (https://www.ensembl.org/biomart).

#### *4.10. Cloning of miRNA Binding Sites and Luciferase Reporter Assay*

We adapted the psi-Check-2 vector (Promega, Madison, WI, USA) to remove the predicted miR-378a-3p target site (7mer-A1) in the open reading frame of the *Renilla* luciferase gene. The binding site was mutated by changing two nucleotides in the seed region without affecting the amino-acid sequence. This was accomplished by substituting the 460 nt long fragment between the EcoRV and XhoI sites of the psi-Check-2 vector (Figure S8; Integrated DNA Technologies, Leuven, Belgium). Effective *Renilla* luciferase production independent of miR-378a-3p levels was confirmed before cloning putative binding sites of target genes.

Ten potential miRNA binding sites of eight miR-378a-3p target genes were ordered as 58-mer oligo duplexes (Integrated DNA technologies) and cloned into the *Xho*I and *Not*I restriction sites of the modified luciferase reporter vector (Table S3). For the binding sites with a positive result in the first luciferase reporter assay, mutant controls were generated by cloning oligo duplexes with mutations in three nucleotides in the seed region. The reporter vectors with miR-378a-3p wildtype or mutated binding sites were co-transfected with 10 µM of either miR-378a pre-miRNA (Cat. NO.: AM17100, Ambion, Austin, TX, USA) or control oligos (Cat. NO.: AM17111, Ambion) to ST486 and DG75 cells using an Amaxa nucleofector device (program A23) and the Amaxa Cell Line Nucleofector Kit V (Cat NO.: VACA-1003) (Amaxa, Gaithersburg, MD, USA). Cells were harvested 24 h after transfection. *Renilla* and firefly luciferase activity in cell lysates was measured using a Dual-Luciferase Reporter Assay System (Promega). Each experiment was measured in duplicate and results were averaged per experiment. For each construct, the luciferase assay was performed in three independent biological replicates.

#### *4.11. Statistical Analysis*

MiRNAs significantly differentially expressed in the small RNA-seq profiling were identified with a moderated *t*-test and Benjamini–Hochberg correction for multiple testing using the GeneSpring GX software (version 12.5, Agilent Technologies, Santa Clara, CA, USA). For confirmation of differentially expressed miRNAs by qRT-PCR, we used the nonparametric Mann–Whitney U-test (GraphPad Software Inc., San Diego, CA, USA). Statistical analysis of GFP competition assays was performed as described previously [55]. Briefly, the percentage of mZip-378a-3p infected cells at day 4 was set to 100% and compared to percentages in the control over time using a mixed model, with time and the interaction of time and miRNA construct type as a fixed effect and the measurement repeat within miRNA construct type as a random effect in SPSS (22.0.0.0 version, IBM, Armonk, NY, USA). For the luciferase reporter assay, significance was calculated on the basis of the *Renilla*-to-firefly (RL/FL) luciferase ratios between experimental samples and negative controls using a paired *t*-test (GraphPad Software Inc.). MNT, IRAK4, and FOXP1 protein levels were normalized to total protein loading as visualized with Ponceau S staining. The average change in protein levels of three to four experiments was calculated relative to pCDH-EV and miRZip-SCR controls, which were set to 1. Significance was calculated using a one-tailed paired *t*-test (GraphPad Software Inc.).

#### **5. Conclusions**

In conclusion, we identified 26 miRNAs differentially expressed between BL cells and GC-B cells and confirmed deregulated expression of five out of eight miRNAs in both BL cell lines and tissue samples. For one of the differentially and highly abundant (top 10 in BL) miRNAs, miR-378a-3p, we showed a negative effect on BL cell growth upon inhibition. Overexpression of two experimentally proven miR-378a-3p targets, i.e., IRAK4 and MNT, phenocopied the effect observed upon miR-378a-3p inhibition. Together, our data show a critical role for miR-378a-3p in promoting BL cell growth and suggest that this involves controlling IRAK4 and MNT levels.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-6694/12/12/3546/s1, Figure S1: Validation of the differential expression pattern observed in the BL cell lines relative to GC-B cells by qRT-PCR, Figure S2: miR-378a-3p levels based on small RNA-seq data derived from individual BL cell lines and germinal center B-cells, Figure S3: Validation of Ago2-IP procedure on ST486 cells upon inhibition and overexpression of miR-378a-3p, Figure S4: Uncropped western blot results for Ago2 protein, Figure S5: Results of the luciferase reporter assays for 10 putative miR-378a-3p binding sites identified in eight Ago-2 IP-enriched genes in (A) ST486 and (B), Figure S6: Uncropped Western blot results for IRAK4, FOXP1, and MNT upon miR-378a-3p inhibition or overexpression, Figure S7: Uncropped Western blot results for IRAK4 and MNT upon lentiviral overexpression, Figure S8: Sequence of the minigene used to generate a luciferase reporter vector without the putative miR-378a-3p binding site in the *Renilla* gene, Table S1: Small RNA sequencing read summary, Table S2: Taqman miRNA assays used for qRT-PCR validation, Table S3: Oligos for cloning miR-378a-3p binding sites from selected genes, wildtype, and with mutations in miR-378a-3p seed. Green letters indicate miR-378a-3p binding sites and red letters indicate the mismatches.

**Author Contributions:** Conceptualization, A.D.-K., J.K. (Joost Kluiver), K.K., and A.v.d.B.; methodology, J.K. (Jasper Koerts), D.d.J., L.W., M.F.H., W.K., T.v.d.S., I.S.-P., M.W., M.M.T., and B.R.; software, F.N., A.D.-K., M.M.T., and J.K. (Joost Kluiver); validation, F.N., J.K. (Jasper Koerts), and D.d.J.; formal analysis, F.N., A.D.-K., and J.K. (Joost Kluiver); investigation, F.N., I.S.-P., A.D.-K., J.K. (Joost Kluiver), and A.v.d.B.; resources, J.K. (Joost Kluiver) and A.v.d.B.; data curation, A.D.-K. and J.K. (Joost Kluiver); writing—original draft preparation, F.N., A.D.-K., J.K. (Joost Kluiver), I.S.-P., K.K., and A.v.d.B.; writing—review and editing, F.N., A.D.-K., J.K. (Joost Kluiver), and A.v.d.B.; visualization, F.N.; supervision, A.D.-K., J.K. (Joost Kluiver), and A.v.d.B.; project administration, F.N.; funding acquisition, F.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the China Scholarship Council (CSC), grant number 601317 (awarded to F.N.) and by a grant from the National Science Center, Poland (2015/19/D/NZ1/03443 to I.S.P.).

**Acknowledgments:** We thank the UMCG Genomics Coordination center, the UG Center for Information Technology, and their sponsors BBMRI-NL & TarGet for storage and computing infrastructure.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Communication*
