Analysis of miRNA‐mRNA regulatory network revealed key genes induced by aflatoxin B1 exposure in primary human hepatocytes

Abstract Background Aflatoxin B1 (AFB1) exposure is a crucial factor to initiate hepatocellular carcinoma (HCC). However, comprehensive microRNA (miRNA)‐message RNA (mRNA) regulatory network regarding AFB1‐associated HCC is still lacking. This work was aimed to identify miRNA‐mRNA network in primary human hepatocytes after AFB1 exposure. Methods A miRNA expression dataset GSE71540 obtained from the gene expression omnibus (GEO) was used to identify differentially expressed miRNAs (DEMs) after AFB1 exposure using GEO2R. Target genes of these DEMs were identified using TargetScan V_7.2, miRDB, PITA, miRanda, and miRTarBase. Gene ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed at Database for Annotation, Visualization and Integrated Discovery (DAVID). miRNA‐mRNA regulatory network was established by analyzing three enriched KEGG pathways significantly correlated with HCC onset and then visualized at CytoScape. Results In this work, nine upregulated and nine downregulated DEMs were identified. Functional enrichment analyses showed that these predicted target genes were significantly associated with cancer development. Analysis of three enriched pathways related to the onset of HCC identified 13 and nine target genes for upregulated DEMs and downregulated DEMs, respectively. Subsequently, the miRNA‐mRNA regulatory networks were constructed. Conclusions In conclusion, miRNA‐mRNA regulatory network was established, which will help to understand the mechanism underlying the AFB1‐induced onset of HCC.


| INTRODUCTION
Liver cancer ranks the fourth place of all newly identified cancer cases but the third place of cancer-related deaths in China (Chen et al., 2016). Hepatocellular carcinoma (HCC) alone is reported to represent about 75%-85% of all liver cancer cases (Bray et al., 2018). Exposure to aflatoxin B1 (AFB1) is considered as a trigger for the initiation of HCC (Llovet et al., 2016). However, genes alteration in primary human hepatocytes after AFB1 exposure remains to be elucidated.
Aflatoxin B1 is believed to be transformed to AFB1-8,9epoxide by a specific P450 enzyme and then interact with cellular proteins or DNA to increase the risk of developing into HCC (Groopman, Kensler, & Wild, 2008). It has also been recognized that AFB1 exposure will cause persistent epigenetic changes (Ferreira et al., 2019;Rieswijk et al., 2016). For example, the mutation of P53 (191170), a well-characterized tumor suppressor gene, is often found in tissues or cells encountered with AFB1, which will lead to the loss function of P53 (Shen & Ong, 1996). Moreover, it was found that exposure of AFB1 will cause aberrant methylation status in genes including RUNT RELATED TRANSCRIPTION FACTOR 3 (600210), LONG INTERSPERSED NUCLEAR ELEMENT 1 (151626), and so on to initiate carcinogenesis . In addition, specifically altered microRNA (miRNA) F I G U R E 1 Volcano Plot of differentially expressed miRNAs.
MicroRNAs are endogenous noncoding RNAs with the length of 18-25 nucleotides, which can directly regulate gene expression through 3'-untranslated region binding (Bartel, 2009). It has been reported that miRNAs are crucial modulators for various cellular behaviors including cell growth, invasion, apoptosis, and so on (Gandellini, Doldi, & Zaffaroni, 2017;Jafri, Al-Qahtani & Shay, 2017). Studies focusing on the investigations of abnormally expressed miRNAs in human cancers revealed that miRNA has crucial roles in tumorigenesis (Gandellini et al., 2017;Jafri et al., 2017).
Here, we are interested to investigate abnormally expressed miRNAs after AFB1 exposure with the aim to understand key miRNAs that can trigger HCC carcinogenesis. Therefore, in this study, we investigated dysregulated miR-NAs in human hepatocytes treated with AFB1. Finally, novel miRNA-message RNA (mRNA) regulatory networks associated with HCC onset were constructed.

| Microarray dataset
GSE71540, obtained from Gene Expression Omnibus (https :// www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE71540), is a miRNA expression profile dataset submitted by Linda Rieswijk. This dataset contained data from primary human hepatocytes subjected to 1 µM of AFB1 for 5 days followed by a 3 days wash-out period or not.

| Identification of differentially expressed miRNAs
GEO2R, an algorithm based on R language limma package, was employed to analyze differentially expressed miRNAs (DEMs) in this dataset. Cut-off criteria were set as p < .05 and |log2 fold-change| ≥ 1.5. Volcano plot was established using the plug-in tool at SangerBox (http://sange rbox.com) to help us identify the overlapping genes. Heatmap was also built by plug-in tool at SangerBox (http://sange rbox.com) to give a direct presentation of the expression status of genes in different samples.

| Identification of target genes for DEMs
Targets of these identified DEMs were predicted at five miRNA target prediction websites including TargetScan V_7.2 (http://  Chou et al., 2018). Target predicted by at least three algorithms were selected for the following studies.

| Functional analyses
Database for Annotation, Visualization and Integrated Discovery (DAVID; https ://david.ncifc rf.gov, Jiao et al., 2012) was employed to perform gene ontology (GO) including cellular component (CC), molecular function (MF), and biological process (BP) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis with the purpose of understanding the biological functions of these genes. Parameter p < .01 was used as threshold.

| Gene expression level validation
The expression level of DEMs and their target was validated at StarBase (Li, Liu, Zhou, Qu, & Yang, 2014).

| miRNA-mRNA regulatory network
construction miRNA-mRNA regulatory network related to the onset of HCC was visualized using CytoScape V_3.6.0 software (Shannon et al., 2003).

| Identification of DEMs
As presented in Figure 1, nine upregulated and nine downregulated DEMs were identified in primary human hepatocytes with AFB1 exposure compared to those without. Moreover, heatmap of identified DEMs was presented in Figure 2. We showed that the samples with or without AFB1 treatment were clearly divided into two groups, suggesting the reliability to analyze DEMs in the following experiments.

| Target genes of DEMs
Five miRNA target prediction algorithms were employed to analyze potential targets of these identified DEMs.
Targets predicted by at least 3 algorithms were selected for followingly analyses. Tables 1 and 2 displayed the targets that contained the most and second-most binding sites for upregulated and downregulated DEMs, respectively.

| Functional analyses of the upregulated or downregulated DEMs
Gene ontology analysis including BP, CC and MF terms for the targets of upregulated or downregulated DEMs were performed to understand the roles of these DEMs. As displayed in Tables 3 and 4, our results revealed that these DEMs have crucial roles in regulating gene expression and cell behaviors.

F I G U R E 4 Validation of the expression pattern of downregulated DEMs in StarBase. DEMs, differentially expressed microRNAs
KEGG pathway enrichment analysis in Figure 3a,b revealed that the most enriched pathways were related to cancer development.

| Validation of DEMs expression in HCC using StarBase
Next, we validated the expression of the DEMs identified in HCC. As presented in Figure 4, the expression pattern of hsa- MIR-505-3p, hsa-MIR-1228-3p, hsa-MIR-4286, and hsa-MIR-5100 of the nine downregulated DEMs was similar to the GEO dataset. For the nine upregulated DEMs, we found that all of them exhibited a similar expression trend in HCC as compared with the GEO dataset ( Figure 5).

| Regulatory network analysis
We were interested to identify miRNA-mRNA networks after AFB1 exposure. Hence, three KEGG pathways related to cancer initiation were selected after pathway enrichment analyses. Three KEGG pathways including miRNAs in cancer, Pathways in cancer, and PI3K-Akt signaling pathway for upregulated DEMs were selected and 13 overlapping genes were identified for following analyses (Figure 6a). The expression of these 13 genes in HCC F I G U R E 5 Validation of the expression pattern of upregulated DEMs in StarBase. DEMs, differentially expressed microRNAs was validated in StarBase, and 10 of them were found to be abnormally expressed (Figure 6b). Similarly, the pathways including Pathways in cancer, Ras signaling pathway, and VEGF signaling pathway for downregulated DEMs were analyzed with nine genes identified as shown in Figure 7a. Also, we confirmed the expression of these 9 genes in HCC using Starbase. We found 6 (MAPK1

| DISCUSSION
Aflatoxin B1 is a highly toxic reagent produced by Aspergillus flavus and the exposure to AFB1 is reported to account for the initiation of approximately 4%-28% of all HCC cases (Liu & Wu, 2010). Therefore, detection of early fingerprints after AFB1 exposure may help us to prevent the cell development into HCC (Fedeles, Chawanthayatham, Croy, Wogan, & Essigmann, 2017). It was well-recognized that miRNAs can regulate almose all cellular processes. Therefore, we hypothesized that miRNAs expression may also be altered in AFB1 treated cells. Importantly, studies comprehensively investigating the connection between miRNAs and AFB1 treatment are still lacking. In this study, a miRNA expression profile dataset that contains data from primary human hepatocytes treated with AFB1 for 5 days and then subjected to 3 days wash-out period was used to investigate the changed miRNAs after AFB1 withdraw. A total of 18 DEMs including nine upregulated DEMs and nine downregulated DEMs were identified. After targets prediction for these identified DEMs, we analyzed the biological functions using GO annotation and KEGG pathway enrichment methods. We found that the pathways of these genes involved were significantly associated with human cancer progression, which indicated that AFB1 exposure will indeed activate the development of cancers. On the basis of these analyses results, we constructed the miRNA-mRNA regulatory network to help us understand the mechanisms related to HCC onset after AFB1 exposure. hsa-miR-96-5p in the upregulated DEMs associated miRNA-mRNA network has the largest target numbers. hsa-miR-34a was previously reported to correlate with liver tumorigenesis after AFB1 exposure (Zhu et al., 2015). Furthermore, hsa-miR-96-5p, hsa-MIR-30a-5p (612326), hsa-MIR-199a-3p (610719), and hsa-  were previously reported to be associated with HCC progression (Callegari et al., 2018;de Conti et al., 2017;Iwai et al., 2018;Zhang, Liu, Zhang, & Liu, 2017). However, the role of hsa-MIR-34b-5p (611374) in HCC has not been reported to date and therefore our study implied that hsa-MIR-34b-5p may be related to HCC onset. For the regulatory network constructed by downregulated DEMs, we have noted that the role of hsa-MIR-4286 was not reported to be associated with HCC previously. Therefore, it deserves to be investigated in the future. This work established miRNA-mRNA regulatory networks related to the HCC onset by comprehensive bioinformatic analyses. The limitation in this work was that the predicted miRNA-mRNA link was not experimentally validated, which need further investigations. Hence, future investigations on in vitro cell lines or in vivo animal models are required to validate the prediction results.
In summary, our in-silico analyses results indicated miRNA-mRNA connections contributing to the HCC onset after AFB1 exposure. We hope that our results can accelerate the process of tackling the AFB1-induced onset of HCC.