The Effects of Naringenin on miRNA-mRNA Profiles in HepaRG Cells

Naringenin, a natural flavonoid widely found in citrus fruits, has been reported to possess anti-oxidant, anti-inflammatory, and hepatoprotective properties as a natural dietary supplement. However, the regulatory mechanism of naringenin in human liver remains unclear. In the present study, messenger RNA sequencing (mRNA-seq), microRNA sequencing (miRNA-seq), and real-time qPCR were used to distinguish the expression differences between control and naringenin-treated HepaRG cells. We obtained 1037 differentially expressed mRNAs and 234 miRNAs. According to the target prediction and integration analysis in silico, we found 20 potential miRNA-mRNA pairs involved in liver metabolism. This study is the first to provide a perspective of miRNA–mRNA interactions in the regulation of naringenin via an integrated analysis of mRNA-seq and miRNA-seq in HepaRG cells, which further characterizes the nutraceutical value of naringenin as a food additive.


Introduction
Naringenin, a natural flavonoid, is found abundantly in citrus fruits and other edible fruits, like grapefruit, oranges, bergamot, tomatoes, and figs [1]. After oral administration, naringenin is widely distributed in the gastrointestinal tract and liver [2,3] and exhibits a direct antioxidant property by free radical scavenging activity on account of its molecular structure, and induces the endogenous antioxidant system in the liver [4]. Growing evidence from both in vitro and in vivo studies has identified various protective capacities of naringenin, such as anti-inflammatory [5], antioxidant [6], anti-fibrosis [7], and hepatoprotective [4,8] activities. These capacities suggest that dietary naringenin could be applied to prevent metabolic syndrome and malignant diseases, including fulminant hepatitis, fatty liver disease, fibrosis, etc. [9,10]. However, there are few detailed studies about the regulatory effects of naringenin on the overall genes of liver [4].
MicroRNAs (miRNAs) are a class of non-coding, single-stranded RNA molecules with a length of 18-26 nucleotides encoded by endogenous genes, which exhibit a broad range of biological regulatory functions in phylogeny, differentiation, proliferation, and apoptosis [11]. miRNA, binding messenger RNAs (mRNAs) by sequence-specific recognition, negatively regulates gene expression at the post-transcriptional level through degradation of target mRNAs [12]. Under exogenous stimulation, miRNA expression is altered, and then target mRNA expression is regulated. Eventually, extensive physiological functions are changed to cope with the challenges caused by exogenous stimulation [13]. A more reliable method for predicting miRNA-mRNA target relations is to simultaneously integrate mRNA-seq analysis with miRNA-seq analysis using a particular processing context in silico [14]. According to the Gene Ontology (GO) classification system, 1037 DEGs were classified into three major functional categories (biological process, cellular component, and molecular function) and 61 subcategories ( Figure 2). Among the biological process, cellular process (731) was the most commonly represented, followed by single-organism process (672) and biological regulation (595). In the category of cellular component, a significant proportion of clusters was assigned to cell (727), cell part (721), and organelle (580). Genes involved in binding (666) and catalytic activity (238) groups were notably represented in the molecular function category. According to the Gene Ontology (GO) classification system, 1037 DEGs were classified into three major functional categories (biological process, cellular component, and molecular function) and 61 subcategories ( Figure 2). Among the biological process, cellular process (731) was the most commonly represented, followed by single-organism process (672) and biological regulation (595). In the category of cellular component, a significant proportion of clusters was assigned to cell (727), cell part (721), and organelle (580). Genes involved in binding (666) and catalytic activity (238) groups were notably represented in the molecular function category.
The Kyoto Encyclopedia of Genes and Genomes (KEGG) classification was found for 1037 DEGs that were further classified into the top twenty biochemical pathways according to the smallest q-value and the largest GeneNumber in pathway annotation ( Figure 3).The GeneNumber and ratio of annotated genes of the top five pathways are systemic lupus erythematosus (33,8.4%), alcoholism (38, 9.67%), transcriptional misregulation in cancers (22,5.6%), PI3K-Akt signaling pathway (34,8.65%) and complement and coagulation cascades (12, 3.05%). There were more than 10 enriched genes among the five pathways. Overall, undergoing naringenin treatment had a significant impact on the global gene expression profile of HepaRG cells. These results implied that the genes involved in these pathways may play crucial roles in naringenin regulation.  The Kyoto Encyclopedia of Genes and Genomes (KEGG) classification was found for 1037 DEGs that were further classified into the top twenty biochemical pathways according to the smallest q-value and the largest GeneNumber in pathway annotation ( Figure 3).The GeneNumber and ratio of annotated genes of the top five pathways are systemic lupus erythematosus (33,8.4%), alcoholism (38, 9.67%), transcriptional misregulation in cancers (22,5.6%), PI3K-Akt signaling pathway (34,8.65%) and complement and coagulation cascades (12, 3.05%). There were more than 10 enriched genes among the five pathways. Overall, undergoing naringenin treatment had a significant impact on the global gene expression profile of HepaRG cells. These results implied that the genes involved in these pathways may play crucial roles in naringenin regulation.

Analysis of miRNA Transcript Levels in Response to Naringenin
In this study, we aimed to determine whether naringenin exposure alters the expression levels of miRNAs in HepaRG cells. After exposure, we collected small RNAs and measured their relative abundance using Illumina HiSeq2500 (Genedenovo Biotechnology Co., Ltd, Guangzhou, China). As shown in Table 2, clean reads of eight samples were generated, respectively, after removing contaminant reads. An overview of reads for small RNA sequencing from raw data to high quality and with quality filtering is provided in Table 2. The length distributions of small RNAs were similar among libraries in that 21-23 nt RNAs were the most abundant ( Figure 4). In this study, we aimed to determine whether naringenin exposure alters the expre sion levels of miRNAs in HepaRG cells. After exposure, we collected small RNAs an measured their relative abundance using Illumina HiSeq2500 (Genedenovo Biotechno ogy Co., Ltd, Guangzhou, China). As shown in Table 2, clean reads of eight samples wer generated, respectively, after removing contaminant reads. An overview of reads fo small RNA sequencing from raw data to high quality and with quality filtering is pro vided in Table 2. The length distributions of small RNAs were similar among libraries i that 21-23 nt RNAs were the most abundant ( Figure 4).    All the small RNAs were aligned in the GeneBank database (Release 209.0) and the Rfam database (11.0) to identify and remove ribosomal RNA (rRNA), small conditional RNA (scRNA), small nucleolar (snoRNA), small nuclear (snRNA), and transfer RNA (tRNA). In accordance with reference genome, these small RNAs, mapped to exons, introns, and repeat sequences, were also removed. The filtering small RNAs were searched against miRBase database (Release 21) to identify miRNAs. The heat map of 3373 miRNAs shows the cluster analysis of the control group and the naringenin group in Figure 5a. Illumina HiSeq2500 profiling of the 3373 miRNAs analyzed in naringenin exposure vs. control samples showed that a total of 234 differentially expressed miRNAs (DEMs, 174 up-and 60 down-regulated; Table S2, Supplementary Materials) were detectable in HepaRG cells (Figure 5b).
Rfam database (11.0) to identify and remove ribosomal RNA (rRNA), small conditional RNA (scRNA), small nucleolar (snoRNA), small nuclear (snRNA), and transfer RNA (tRNA). In accordance with reference genome, these small RNAs, mapped to exons, introns, and repeat sequences, were also removed. The filtering small RNAs were searched against miRBase database (Release 21) to identify miRNAs. The heat map of 3373 miRNAs shows the cluster analysis of the control group and the naringenin group in Figure 5a. Illumina HiSeq2500 profiling of the 3373 miRNAs analyzed in naringenin exposure vs. control samples showed that a total of 234 differentially expressed miRNAs (DEMs, 174 up-and 60 down-regulated; Table S2

Target Prediction and Integration Analysis of mRNA and miRNA Expression Profiles in Response to Naringenin
Acting at the post-transcriptional level, miRNAs silence and/or down-regulate cellular mRNA gene expression by target RNA cleavage. To predict the target genes of 234 DEMs, we performed computational analyses using the RNAhybrid (v2.1.2) + svmlight (v6.01), Miranda (v3.3a), and TargetScan (Version: 7.0). The simultaneous profiling of 234 DEMs and 1037 DEGs levels in silico can identify the presumptive target mRNAs of miR-NAs. We selected the intersection of DEGs and target genes of DEMs, and then performed bioinformatics analysis on these intersection genes. A total of 5607 negative miRNA-mRNA pairs for naringenin treatment were obtained, with the involvement of 216 DEMs and 681 DEGs (Table S3, Supplementary Materials). In line with the GO classification sys-

Target Prediction and Integration Analysis of mRNA and miRNA Expression Profiles in Response to Naringenin
Acting at the post-transcriptional level, miRNAs silence and/or down-regulate cellular mRNA gene expression by target RNA cleavage. To predict the target genes of 234 DEMs, we performed computational analyses using the RNAhybrid (v2.1.2) + svmlight (v6.01), Miranda (v3.3a), and TargetScan (Version: 7.0). The simultaneous profiling of 234 DEMs and 1037 DEGs levels in silico can identify the presumptive target mRNAs of miRNAs. We selected the intersection of DEGs and target genes of DEMs, and then performed bioinformatics analysis on these intersection genes. A total of 5607 negative miRNA-mRNA pairs for naringenin treatment were obtained, with the involvement of 216 DEMs and 681 DEGs (Table S3, Supplementary Materials). In line with the GO classification system, 681 DEGs were classified into 58 subcategories ( Figure 6). The first three subcategories had not changed in three major functional categories compared with transcriptome sequencing analysis (Figures 2 and 6).
Pathway enrichment analysis for 681 DEGs of 5607 negative miRNA-mRNA pairs identified the top twenty pathways according to the smallest q-value and the largest GeneNumber in pathway annotation after naringenin exposure ( Figure 7): with PI3K-Akt signaling pathway (25,8.96%) being the eighth pathway and the second-most abundant. tem, 681 DEGs were classified into 58 subcategories ( Figure 6). The first three subcategories had not changed in three major functional categories compared with transcriptome sequencing analysis (Figures 2 and 6).

Real-Time qPCR Validation of Naringenin Regulation in Liver Metabolism and Potential Regulatory miRNA-mRNA Pairs
According to global gene function annotations, literature review, and their potential

Discussion and Conclusions
The findings discussed here reveal the first detailed information regarding parallel mRNA and miRNA expression changes in HepaRG cells in response to naringenin. We performed an integrative analysis of these data including 234 DEMs and 1037 DEGs induced by naringenin, which provide global insight into the miRNA-mRNA interactions of naringenin in the regulatory mechanism. According to the gene function annotations and literature review, 19 DEGs related to metabolism were screened out. In particular, the PI3K-Akt signaling pathway was significantly enriched both in analysis of transcriptome sequencing (the third-most abundant, and ranked fourth) and integration analysis of miRNA-mRNA expression profiles (the second-most abundant, and ranked eighth) in responses to naringenin. In addition, 11 DEGs in the PI3K-Akt signaling pathway were further validated using real-time qPCR analysis. In this work, we constructed a miRNA-mRNA regulatory network according to the DEMs and DEGs datasets and miRNA-targeting information. Some studies have demonstrated that the miRNA-mRNA regulatory network responds to liver damage, including hepatocellular carcinoma and oxidative stress [15]. Although several miRNA-induced RNA activation phenomena were identified [16], under most circumstances, the negative correlation between miRNAs and their target mRNAs is often considered support for miRNA targeting [17]. Ultimately, 20 miRNA-mRNA negative correlation pairs were identified with the involvement of liver metabolism and the PI3K-Akt signaling pathway.
With regard to global genes, we addressed our particular research question using pathway analysis to highlight 19 DEGs related to the functional clusters: metabolism, including Lipid metabolism (ALOX15, ACSL5, PLA2G4C, and B4GALT6), Energy metabolism (CA9 and NDUFA4L2), Metabolism of cofactors and vitamins (TH, LIPT2, and

Discussion and Conclusions
The findings discussed here reveal the first detailed information regarding parallel mRNA and miRNA expression changes in HepaRG cells in response to naringenin. We performed an integrative analysis of these data including 234 DEMs and 1037 DEGs induced by naringenin, which provide global insight into the miRNA-mRNA interactions of naringenin in the regulatory mechanism. According to the gene function annotations and literature review, 19 DEGs related to metabolism were screened out. In particular, the PI3K-Akt signaling pathway was significantly enriched both in analysis of transcriptome sequencing (the third-most abundant, and ranked fourth) and integration analysis of miRNA-mRNA expression profiles (the second-most abundant, and ranked eighth) in responses to naringenin. In addition, 11 DEGs in the PI3K-Akt signaling pathway were further validated using real-time qPCR analysis. In this work, we constructed a miRNA-mRNA regulatory network according to the DEMs and DEGs datasets and miRNA-targeting information. Some studies have demonstrated that the miRNA-mRNA regulatory network responds to liver damage, including hepatocellular carcinoma and oxidative stress [15]. Although several miRNA-induced RNA activation phenomena were identified [16], under most circumstances, the negative correlation between miRNAs and their target mRNAs is often considered support for miRNA targeting [17]. Ultimately, 20 miRNA-mRNA negative correlation pairs were identified with the involvement of liver metabolism and the PI3K-Akt signaling pathway.
The PI3K-Akt signaling pathway may offer clues for the molecular mechanism involved in metabolism, inflammation, and oxidative stress [37], which plays a pivotal role in the response to naringenin. The PI3K-Akt signaling pathway has diverse downstream effects on cellular metabolism through either direct regulation of nutrient transporters and metabolic enzymes or the control of transcription factors that regulate the expression of key components of metabolic pathways [38,39], including glucose metabolism, biosynthesis of macromolecules, and maintenance of redox balance. It was reported that PDGFRB silencing inhibited the activation and proliferation of hepatic stellate cells and ameliorated liver fibrosis [40]. The collaborative inhibition of CSF1R and FGFR2 is expected to enhance the antitumor effects by targeting immune evasion and angiogenesis in the tumor microenvironment [41]. Knocking down ITGB4 suppressed glycolysis in cancer-associated fibroblasts [42]. Genetic knockdown of PCK1 prevented fatty-acid-induced rise in oxidative flux, oxidative stress, and inflammation [43], which was also correlated with signaling pathways governed by insulin [44]. It was shown that overexpression of nuclear CREB3L3 induced systemic lipolysis, hepatic ketogenesis, and insulin sensitivity with increased energy expenditure [45]. Inhibition of CREB3L1 reportedly blocked cancer invasion and metastasis [46]. NFκB is a key regulator of immune development, immune responses, inflammation, and cancer [47]. It has been well-established that suppressing NFκB transduces anti-inflammatory signals and reduces inflammation [48]. In the PI3K-Akt signaling pathway, our results showed that naringenin significantly down-regulated the mRNA expressions of PDGFRB, PCK1, CREB3L1, and NFκB1 with related miRNAs (has-miR-1306-5p and hsa-miR-627-3p) being significantly down-regulated. Therefore, our data suggested that naringenin may play a salutary role in anti-inflammatory, anti-oxidative stress, and ameliorative metabolism via the inhibition of the PI3K-Akt signaling pathway.
This study was the first to integrated the analysis of mRNA-seq and miRNA-seq in the liver in response to naringenin, and provide a perspective of metabolism in naringenin regulation. In terms of metabolism and the PI3K-Akt signaling pathway, the 11 DEMs, 30 DEGs, and 20 miRNA-mRNA pairs need more research for the activities of naringenin. There are some limitations of this research; for example, the dose-dependent and timedependent effects of naringenin in miRNA-mRNA interactions were still unclear. Although given miRNAs analyzed in silico suggested regulatory capacities, their functions need to be further certified in a specific context within a living system. In summary, we provided preliminary research analyzing mRNA and miRNA expression and profiling of metabolism. The regulatory mechanism of miRNA-mRNA pairs could be additional possible evidence for annotating the nutraceutical value of naringenin.

Chemicals and Reagents
The HepaRG cell line was originally purchased from Biopredic International (Rennes, France). RPMI-1640 medium and penicillin-streptomycin-glutamine solution were obtained from Gibco (Gaithersburg, MD, USA). Fetal bovine serum (FBS) was purchased from Corning (Auckland, New Zealand). Naringenin and dimethyl sulfoxide (DMSO) were from Sigma-Aldrich (St. Louis, MO, USA). TRIzol™ reagent was supplied by Thermo Fisher (Carlsbad, CA, USA). Ultrapure water was purified by a Milli-Q academic water purification system (Millipore, Bedford, MA, USA). All other reagents were commercialized products of the highest analytical grade available.

Total RNA Extraction
HepaRG cells were washed twice with ice-cold phosphate-buffered saline (PBS) and harvested with TRIzol™ reagent as recommended by the manufacturer. Thermo Scientific NanoDrop™ 2000 c Spectrophotometers (Wilmington, DE, USA) were used to measure the RNA quality and quantity of each sample according to the manufacturers' protocol.

RNA Sequencing
The mRNA was enriched by Oligo(dT) beads, then the enriched mRNA was fragmented and reverse-transcripted into cDNA with random primers by QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands). RNA molecules in a size range of 18-30 nt were enriched by polyacrylamide gel electrophoresis. The 3 and 5 adapters were added, then enriched RNAs were reverse-transcripted by the QiaQuick PCR extraction kit, according to the manufacturer's instructions (Qiagen, Venlo, The Netherlands). The ligation products were size-selected by agarose gel electrophoresis. There were four samples in the naringenin group and four samples in the control group. Each sample generated two cDNA libraries: one for mRNA-seq and the other for miRNA-seq. PCR amplified products were enriched to respectively generate 16 cDNA libraries and sequenced using Illumina HiSeq2500 by Genedenovo Biotechnology Co. (Guangzhou, China). The RNA and small RNA sequencing data were deposited in the NCBI Sequence Read Archive (accession numbers from SRR13675952 to SRR13675963).

Real-Time qPCR
Transcription of mRNA into cDNA was conducted with the GoScript™ Reverse Transcription System (Promega, Madison, WI, USA) from 3 µg of total RNA, according to the manufacturer's instructions. For miRNA analysis, cDNA-synthesis was performed with the miRNA First Strand cDNA Synthesis (Tailing Reaction, Sangon Biotech, Shanghai, China) from 2 µg of total RNA. The real-time qPCR was carried out with GoTaq ® qPCR Master Mix (Promega, Madison, WI, USA) on a LightCycler 480 (Roche, Mannheim, Germany), as recommended by the manufacturer. The thermal cycling procedure started with an initial denaturation at 95 • C for 10 min. This was followed by 45 cycles of denaturation for 10 s at 95 • C, primer binding for 20 s at 60 • C, and elongation for 20 s at 72 • C. The procedure ended with a final amplification at 95 • C for 5 s, 65 • C for 1 min, the addition of a dissociation curve step, and a cooling step. Primers were purchased from Sangon Biotech (Shanghai, China). The primer pairs' sequences used for the validation of the signature are described in Tables 3 and 4.

Bioinformatic Analysis and SStatistics
DEGs and DEMs were identified using an R-based software package. The threshold value for selection of DEGs and DEMs was q-value (adjusted p-value) ≤ 0.05 and fold change (FC) ≥2 or ≤0.5. KEGG and GO classification including molecular functions, biological processes, and cellular components were used to analyze DEGs and DEMs.
The results of biological assay are presented in the form of mean ± SD based on four independent experiments in GraphPad Prism 8.