Comprehensive Analysis of the Expression Profiles of Hepatic lncRNAs in the Mouse Model of Alcoholic Liver Disease

Background and Aim: The worldwide prevalence of alcoholic liver disease (ALD) due to escalating alcohol consumption has presented an unprecedented pressure on human health. A few studies have determined long non-coding RNAs (lncRNAs) involved in the pathogenesis of liver diseases. However, the roles of lncRNAs in ALD development is still poorly understood. Methods: An ALD mouse model was established and confirmed. Expression profiles of lncRNAs were obtained by whole transcriptome sequencing. The altered lncRNAs in ALD mice were further verified by qRT-PCR. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used to enrich the functions of these lncRNAs. In combination with miRNA and mRNA profiles, we constructed concise endogenous RNA (ceRNA) networks. The function of the most up/downregulated lnRNA was further verified and investigated in both ALD model and AML-12 cells. Results: Totally, five downregulated lncRNAs were obtained and verified in ALD mice. The GO term and KEGG pathway analyses revealed that the identified lncRNAs were associated with alcohol-induced hepatic oxidative damage, cellular inflammation, and lipid metabolism. Combination the differentially modulated miRNAs and mRNAs with ceRNA network analysis, we constructed five ceRNA networks and obtained 30 miRNAs and 25 mRNAs that may participate in ALD. Further, we verified and investigate the function of the most downregulated lnc_1700023H06Rik. Depletion lnc_1700023H06Rik reduced genes encoding for lipid metabolism, especially mRNA Acat2 (ENSMUST00000159697) and Pgrmc2 (ENSMUST00000058578) both in vivo and in vitro. Knocking down lnc_1700023H06Rik induced triglyceride accumulation and lactate dehydrogenase leakage in AML12 cells, consisting with that in alcohol-treated cells. Conclusion: The five remarkably downregulated lncRNAs in ALD mouse model were identified as novel biomarkers, highlighting the key role of lncRNAs in the development of ALD. The effect of lnc_1700023H06Rik plays a pivotal role in lipid deposition and its pathological pathway in ALD needs further investigation.


INTRODUCTION
Continued and excessive alcohol consumption is exacerbating the worldwide incidence of alcoholic liver disease (ALD) (Rehm et al., 2013), According to the Global Status Report on Alcohol and Health (2018), alcohol consumption caused approximately 3 million deaths and nearly half can be attributed to ALD (Lv et al., 2020;Shield et al., 2020). The pathological progression of ALD comprises a spectrum of diseases from liver steatosis and hepatitis to severe fibrosis or even cirrhosis (Liangpunsakul et al., 2016). Although ALD is wellcharacterized, it cannot be halted or reversed because the underlying mechanisms for ALD initiation or progression are still unclear.
In the past few decades, the proposed mechanisms for ALD have mainly focused on protein-coding genes of alcohol metabolism, reactive oxygen species, formation and inflammation responses (Beier and Arteel, 2010). Technologies such as RNA sequencing and microarray have revealed that non-coding RNAs (ncRNAs) compring up to 80% of the "transcriptional noise," and thus, also participated in ALD evelopment. However, the role of lncRNAs in ALD is still inconclusive as compared to that of the well-studied endogenous miRNA.
Long non-coding RNAs (lncRNAs, with lengths exceeding 200 nucleotides) are transcribed from various genomic regions, including introns and exons. The competitive endogenous RNA (ceRNA) regulates target gene expression through chromatin remodeling, transcription, and post-transcriptional processing Zhao et al., 2018). They play a critical role in various biological functions and disease processes (Peng et al., 2017). Recently, lncRNAs has been reported to involve in alcohol-related diseases. For instance, serum levels of AK054921 and AK128652 were downregulated in patients with alcohol-cirrhosis, and they were inversely correlated with the survival of these patients, indicating the pathological role in alcoholic cirrhosis (Yang et al., 2017). The lncRNAs NEAT1, Gm5091, and MEG3 participated in alcoholinduced hepatic steatosis and fibrosis in animal models Zhou et al., 2018;Ye et al., 2020). Despite the definite role of lncRNAs in alcoholic diseases, the profiles of lncRNAs in ALD models have not been comprehensively investigated.
This study objective to screen the potential lncRNAs that may involved in the development of ALD in mice. We performed RNA sequencing to identify the differentially expressed lncRNAs in liver. The functions of lncRNAs were analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Several differentially expressed lncRNAs were further verified. A ceRNA network analysis was constructed to understand the lncRNA-miRNA-mRNA crosstalk. In addition to the miRNA sequencing assay, we proposed several key lncRNA pathways that may participate in the development of ALD. Our results may provide new insights into the underlying mechanisms of ALD.

Animal Experiments
Male C57BL/6J mice weighing 25 ± 0.5 g (mean ± SD) were maintained at the Animal Center of Zhejiang Chinese Medical University. The ALD mouse model was established as previously described (Ding et al., 2017). Mice were divided into two groups of six animals each. The control groups, pair-fed (PF), were fed with an isocaloric control liquid diet (Bioserv, Frenchtown, NJ) and the treated groups, alcohol-fed (AF), were treated with Lieber-DeCarli diet (Bioserv, Frenchtown, NJ) for 5 weeks. In the first 3 days, mice were fed with Lieber-DeCarli diet without alcohol. At the following 2 days, the caloric content of alcohol in the diet was set at 5.5%, and increased to 11% on the 6th and 7th day, 22% at the 2nd week, 27% at the 3rd week, and maintained at 32% at the last 2 weeks. Animals were allowed ad libitum access to diet and water. Mice were euthanized using an intravenous injection of pentobarbital sodium (Merck, Darmstadt, Germany) without any pain after 12 h fasting.

RNA Interference
Based on the results of qRT-PCR, we selected the lncRNA with the most significant differential expression to investigate its function in AML12 cells. Cultured cells were transfected with small interfering RNA (siRNA) against mouse lnc_1700023H06Rik (Ribobio, Guangzhou, China) using Lipofectamine 2000 (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. In the control group (NC), cells were transfected with scrambled siRNA (Ribobio). Cells were collected after 24-48 h for qRT-PCR. Sequences of siRNA against mouse lnc_1700023H06Rik were listed in Supplementary Table S1, the sequences of scrambled siRNA were not disclosed.
To investigate the function of lnc_1700023H06Rik in the setting of alcohol treating. Cells were transfected with siRNAs as described aforementioned. Then, 200 μmol/L of ethanol was added into each well. After 12 h, the culture supernatants were for detection of Lactate dehydrogenase (LDH) and the cells were collected for triglyceride (TG) detection.

Histological and Biochemical Assays
Small pieces of liver were isolated and fixed immediately with 4% paraformaldehyde (Biosharp, Guangzhou, China). Samples were stained with hematoxylin (Sigma) and eosin (Sigma) to determine the histological change. Oil red O (Sigma) staining was performed on fresh liver samples to detect the lipid depostion according to our previous study (Ma et al., 2019). The amounts of liver cholesterol (TC), TG, plasma alanine aminotransferase (ALT), and aspartate aminotransferase (AST) were determined according to the manufacturer's instructions (Jiancheng, Nanjing, China).

Whole Transcriptome Sequencing
Total RNA was extracted, and quantified using Bioanalyzer 2100 (Agilent, CA, United States) with RNA Integrity Number (RIN) > 7.0. Subsequently, a small RNA library and a chain-specific library without ribosomal RNA were established and sequenced. The ribosome-deficient chain-specific library obtain the sequence information of lncRNA but also the sequence information of mRNA and circRNA. The small RNA library provided a miRNA sequence.
Transcripts that overlapped with known mRNAs and those that were shorter than 200 bp were discarded. The coding potential of transcripts was determined using Coding Potential Calculator (CPC) (Kong et al., 2007), Coding-Non-Coding Index (CNCI) (Sun et al., 2013), and Pfam (Punta et al., 2012). Transcripts with CPC scores <−1 and CNCI scores <0 was removed. StringTie and fragments per kilobase of transcript per million fragments (FPKM)  were used to determine the levels of mRNAs and lncRNAs, respectively. LncRNAs were analyzed using the Ballgown R package. Those with log 2 fold change (FC) > 1 or log 2 FC <−1 and p value <0.05 were considered differentially expressed (Frazee et al., 2015).

Validation of Key lncRNAs and mRNA Using qRT-PCR
Total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA). The content and purity was determinated by, NanoDrop ND-1000 (NanoDrop, Wilmington, DE). The reaction systems containing cDNA templates, SYBR Green (Bimake, Houston, TX), ddH 2 O (Bimake), and primers were mixed for qRT-PCR on ABI 7300 PCR instrument (Thermo Fisher Scientific, Waltha, MA). The gradient temperatures were set at 95°C for 5 min for Hot-Start DNA Polymer; 95°C for 15 s, 60°C for 1 min (40 cycles) for PCR; 95°C for 15 s, 60°C for 1 min, and 95°C for 15 s for melt curve. The relative gene expression was calculated using the 2 −ΔΔCt method, and it was normalized to the control group. All primer sequences were listed in Supplementary Table S2-S3.

Gene Ontology Enrichment Analysis
Gene Ontology functional enrichment analysis maps and calculates the numbers of all the differentially expressed genes to each term in the GO database, Then, the enrichment fractions and the number of differentially expressed genes (S gene number) enriched in each term were determined according to p values. Finally, the top 10 terms with p < 0.05 and largest numbers of S gene number were analyzed.

Kyoto Encyclopedia of Genes and Genomes Enrichment Analysis
Biological functions of differentially expressed lncRNAs were determined by KEGG (www.genome.jp/kegg), a public database for genome deciphering. The basic threshold for significant enrichment was set at p < 0.05. Subsequently, each pathway was sorted in a descending order according to S gene numbers, and the top ten pathways were displayed.

ceRNA Network Construction
After verification, lncRNAs with absolute values of log 2 FC > 2 and p < 0.05 were selected. Combination with the differentially expressed miRNA and mRNA in sequence analysis we constructed the ceRNA networks. The RNAs which didn't obeyed the rule of negative regulation between lncRNA and miRNA, miRNA and mRNA, were removed. Finally, ceRNA network diagrams were prepared using Cytoscape (v3.8.2, http:// www.cytoscape.org/) (Liu et al., 2019).

Detection of Triglyceride and Lactate Dehydrogenase Leakage
The supernatant of AML12 was collected and the activity of LDH leakage was detected according to the instructions of LDH assay kit (Beyotime, Shanghai, China). Cells were then washed with PBS and the triglyceride content in each group was measured according to the instructions of TG assay kit (Jiancheng, Nanjing, China).

Statistical Analysis
Statistical analyses were carried out using the Student's t-test and ANOVA-test accordingly via the SPSS 25.0 software (Chicago, IL, United States). Significant level was set at p < 0.05. The bar charts were prepared using GraphPad Prism 7.0 (GraphPad Software, California, United States).

Hepatic Lipid Accumulation and Hepatic Injury Under Alcohol Feeding
The successful establishment of ALD mouse model was verified by a sequence of well-characterized endpoints. As described, the amounts of TC ( Figure 1A), TG ( Figure 1B), ALT ( Figure 1C), and AST ( Figure 1D) remarkably increased in the AF group as compared to those in the PF group, indicating liver injury and hyperlipidemia. Additionally, a moderate amount of lipid droplets and obvious lipid vacuoles were observed in the liver ( Figure 1E). The pathological morphology of the liver pointed toward hepatic steatosis and hepatitis, implicating a well-established ALD mouse model.

Differentially Expressed lncRNAs Identified in Alcoholic Liver Disease Mice
Totally, twenty-nine differentially expressed lncRNAs were identified by using RNA sequencing. We divided them into two groups: twelve upregulated and Seventeen downregulated, as separated and clustered by volcano graph and heat map (Figure 2A,B). In summary, eighteen lncRNAs with eleven upregulated and seven downregulated genes changed nearly by 2 to 2.8 fold. Three lncRNAs (one upregulated and two downregulated) altered by 2.8 to 4-fold. Six lncRNAs, three upregulated and three downregulated, changed by 4 to 8-fold. The remaining two lncRNAs changed nearly by 8 to 16-fold ( Figure 2C). Finally, eight lncRNAs with log 2 FC change >2.0 (4fold) and p < 0.05 were sent for PCR verification ( Table 1).

Verification 5 Differentially Expressed lncRNAs via qRT-PCR
The differentially expressed lncRNAs were further verified by RT-PCR. Five of the selected lncRNAs in the AF group were maintained as the PF group. Among which, five (mmu_lnc_AW495222(39,807), mmu_lnc_1700023H06Rik, mmu_lnc_0610005C13Rik, mmu_lnc_Gm12265, and  Frontiers in Pharmacology | www.frontiersin.org July 2021 | Volume 12 | Article 709287 mmu_lnc_Gm45724) were consistent with the sequencing results, and two were reversed (mmu_lnc_Gm38357 and mmu_lnc_AW495222(39,805)) ( Figure  3). Notably, mou_lnc_AW495222 (39,807) and mou_lnc_Gm45724 demonstrated maximum inhibition both in RNA sequencing and PCR verification assays. The expression of mmu_lnc_Rian was not different in PF and AF groups. Altogether, about 5 lncRNAs were verified in RT-PCR.

Functional Enrichment of Differentially Expressed Genes
Functions of the differentially expressed lncRNAs were separately determined using GO and KEGG enrichment analyses. The first ten pathways enriched from BP, CC, and MF terms were selected using GO analysis. The enriched BP term ( Figure 4A) included metabolic processes of lipids, fatty acids, steroid hormones, and oxidation-related pathways with rich factors higher than 0.1. The epoxygenase P450 pathway was ranked the first rich factor in the BP term. All the enriched CC terms had a high rich factor (>0.6) with the most significantly enriched pathway as intermediate density lipoprotein particle ( Figure 4B). The MF terms were enriched in oxidative metabolism enzymes, aromatase, and lyase activity as well as iron and heme binding ( Figure 4C). Specifically, the rich factor of aromatase activity peaked near 0.6.
The potential pathogenesis of ALD was determined using KEGG enrichment analysis with a rich factor higher than 0.1. In total, ten pathways were enriched, seven of which had rich factors higher than 0.2. They were ranked as retinol metabolism, arachidonic acid metabolism, steroid hormone biosynthesis, metabolism of xenobiotics by cytochrome P450 (CYP450), drug metabolism by CYP450, peroxisome proliferators activated receptor (PPAR) signaling pathway, and glutathione metabolism ( Figure 4D).
The ceRNA interaction network of 5 differentially expressed lncRNAs.
The ceRNA network represents a novel regulatory mechanism among lncRNA, miRNA, and mRNA. We chose the top five significantly downregulated lncRNAs (mou_lnc_0610005C13Rik, mou_lnc_1700023H06Rik, mou _ lnc _ Gm12265, mou _ lnc _ AW495222(39,807), and mou_lnc_Gm45724) to construct the ceRNA network. The complex networks were further modified  Frontiers in Pharmacology | www.frontiersin.org July 2021 | Volume 12 | Article 709287 according to our whole transcriptome sequencing analysis data comprising thirty upregulated miRNAs and twenty-five downregulated mRNAs (Supplementary Table S4). Finally, the five aforementioned lncRNAs targeted eight miRNAs and eleven mRNAs ( Figure 5A), ten miRNAs and eleven mRNAs ( Figure 5B), five miRNAs and four mRNAs ( Figure 5C), four miRNAs and seven mRNAs ( Figure 5D), and five miRNAs and eight mRNAs ( Figure 5E), respectively. Verification on lnc_1700023H06Rik ceRNA network in both in vivo and in vitro model.

Interference lnc_1700023H06Rik Resulted in Triglyceride Accumulation and Cell Damage
The biological function of lnc_170023H06Rik was further investigated in AML-12 cells. In line with the PCR and RNAseq analysis, lnc_170023H06Rik was down-regulated in AML-12 cells after exposure to alcohol ( Figure 7A). Depletion lnc_170023H06Rik was accompanied with apparent TG accumulation and LDH leakage in AML12 cells, which was consistent with that of ethanol treatment ( Figure 7B,C).

DISCUSSION
In recent decades, many studies have been conducted to determine the pathology of ALD. Recently, the cardinal role of lncRNAs in ALD has been recognized. A few lncRNAs have been identified as biomarkers for predicting survival in patients with alcoholic cirrhosis (Yang et al., 2017). For the first time, we FIGURE 4 | GO and KEGG enrichment analysis of differentially expressed lncRNAs. BP (A), CC (B), and MF (C) terms for enrichment (rich factor) of differentially expressed lncRNAs in GO analysis ("*a, b, c" refers to "oxidoreductase activity, acting on paired donors, and molecular oxygen, respectively.") (D) KEGG pathways enriched by differentially expressed lncRNAs. The size of the circle represents S gene numbers. The color depth represents the statistical significance.
Frontiers in Pharmacology | www.frontiersin.org July 2021 | Volume 12 | Article 709287 prioritized and confirmed eight new lncRNAs that may be involved in the pathological of ALD. Furthermore, we established the top five lncRNA-modulated networks based on the results of whole transcriptome sequencing. Based on the co-expressed mRNAs functions of differentially expressed lncRNAs could be mainly ascribed to the following categories: oxidative stress, lipid metabolism and inflammatory reactions. According to the KEGG enrichment analysis, we determined the metabolism of xenobiotics by CYP450, drug metabolism by CYP450, retinol metabolism, and arachidonic acid metabolism as the most enriched pathways. Consistent with our present and previous findings, alcohol intake triggers RNAs involved in xenobiotic and drug biodegradation pathways (Dou et al., 2020). CYP450 is crucial for chemical metabolism. Alcohol-related hepatic oxidative stress is caused by ethanol oxidation by its superfamily CYP2E1, which ultimately promotes the development of ALD (Lu et al., 2018). Normally, the liver would exert anti-oxidative and anti-inflammatory processes by converting arachidonic acid into biologically active cyclic trienoic acid through arachidonic acid epoxygenase (Wells et al., 2016). However, CYP450 activation induces arachidonic acid hydroxylation and epoxidation, which partially lead to oxidative stress and liver peroxidation after alcohol intake (Seitz and Mueller, 2019). Additionally, alcohol can cause the decrease of retinol content in liver, and activation CYP2E1 is the main reason for this effect (Clugston et al., 2015). Results from the KEGG analysis were consistent with those of the GO analysis, as the most enriched CC, BF, and MF terms were occupied by lipid and alcoholic metabolism.
Although the functions of lncRNAs have not been fully recognized, the interacting with miRNAs and mRNAs from the ceRNA network enable the understanding of the pathological roles of lncRNAs in ALD. As mentioned above, the initial pathological feature of ALD is hepatic steatosis, that is, FIGURE 5 | ceRNA networks of mou_lnc_0610005C13Rik, mou_lnc_1700023H06Rik, mou _ lnc _ Gm12265, mou _ lnc _ AW495222(39,807), and mou_lnc_Gm45724. Blue hexagon represents the five lncRNAs. The green triangles and orange quadrilaterals represent differentially expressed miRNA and mRNA, respectively.
Frontiers in Pharmacology | www.frontiersin.org July 2021 | Volume 12 | Article 709287 ALD is closely related to lipid metabolism. This study found that lnc_0610005c13Rik-miR540, lnc_Gm45724-miR200c-3p, lnc_1700023H06Rik-miR294-3p, miR200c-3p, miR200b-3p, lnc_12,265-miR294-3p, and -miR 330-5p regulated the mRNA of Acox1, the first enzyme in fatty acid oxidation. Moreover, induction of miR540 inhibited the proteins involved in fatty acid oxidation, such as Acox1, PPAR⍺, and cpt1⍺, in ALR-Lknockout mice (Kumar et al., 2019). And the expression of Acox1 was down-regulated in alcohol-treated mice with increased lipid accumulation in the mouse liver . Thus, mitigating these lncRNAs remarkably blocked the activities of enzymes involved in fatty acid oxidation, leading to  Frontiers in Pharmacology | www.frontiersin.org July 2021 | Volume 12 | Article 709287 8 lipid accumulation. These predictions were further verified in lnc_1700023H06Rik. Both lnc_1700023H06Rik and its related mRNAs, especially Acat2 and Pgrmc2 were downregulated in ALD mice and alcohol-treated cells. Genetic depletion of Acat2 and Pgrmc2 leads to excessive deposition of free fatty acid in vivo, and lipid metabolic dysfunction Galmozzi et al., 2019). Ectopic Lipid accumulation and the consequence of lipotoxicity is the hallmark for liver injury in alcohol feeding mice (Alpini, 2019). Our study in vitro confirmed the blunted lnc_1700023H06Rik tightly contributed to lipid accumulation and cell injury when response to alcohol treatment.
Besides lipid deposition, oxidative stress and inflammation are two pathological factors in ALD. The mRNAs including Prdx2, slc38a6, and sult2a8, engaging in the network of lnc_1700023H06Rik-miR223-3p and lnc_AW49522-miR494-3p are responsible for the anti-oxidative and glutamate transportation pathways (Chan et al., 2016;Shimohira et al., 2018;Wang et al., 2020). Diminishing Prdx2 in hepatocyte of alcohol-treated mice endows lipid deposition (Han et al., 2020). Thus, attenuation of both lnc_1700023H06Rik and lnc_AW49522 in ALD mice indicated the deterioration of anti-oxidative system and lipid metabolism. Additionally, downregulation of lnc_Gm12265, lnc_Gm45724 and lnc_1700023H06Rik may involve in inflammation signaling. Activation of miR326 by inhibition the expression of lnc_Gm12265 may promoted hepatic inflammation as miR326 participates in the secretion of inflammation cytokines via TLR4/myD88/ NFκB signaling (Liao et al., 2019). In contrast, miR223-3p, an emerging negative regulator of NLRP3 inflammation and inactivator of hematopoietic stem cells (Jimenez Calvente et al., 2020), was also triggered in the ALD mouse model through the downregulation of lnc_Gm45724 or lnc_1700023H06Rik. Thus, we posited that miR223-3p may participated in ALD given that NLRP3 is responsible for innate immunity during acute liver injury (Huang et al., 2013). The mRNA of Akt2, engaging in the network of lnc_0610005C13Rik, may also contribute to the inflammation status in ALD mice since inhibition on Akt2 attenuated the transcription of inflammation cytokine under alcohol challenge (Reyes-Gordillo et al., 2019).
The limitation of this study lies in the low sequence alignment between mice lncRNAs and the whole human genome via the sequence homological alignment. Most of lncRNAs that identified by our studies cannot be successfully paired in the sequences of humans. Only mmu _ lnc _ 170023H06Rik has a high homology between mice and human (ARRDC3-AS1, 89.58%, 600 bases, Supplementary  Table S5). About 20 to 30 bases of mmu_lnc_Rian in mice has a full homology with human lncRNAs (LINC02315, LINC02307). Unexpectedly, none of the compared three human lncRNAs (ARRDC3-AS1, LINC02315, LINC02307) has been reported to be related to ALD. Despite this limitation, our study still sets the stage for further investigation on ALD development because the pathological mechanism of ALD in human is still in its infant stage.
In conclusion, the present study identified five lncRNAs functioning as anti-inflammation, anti-oxidation, and lipid metabolism in ALD mice model. Combination ceRNA network with miRNAs and mRNAs sequencing data, we got five ceRNA networks. Taken the most downregulated lncRNA as an example, we, for the first time, verified lnc_1700023H06Rik is tightly related to lipotoxicity in ALD mice model. However, the exact activities of the identified lncRNAs and the proposed pathways still need verification. Despite this limitation, our study highlights the identified lncRNAs as key factors responsible for ALD pathogenesis.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/geo/, GSE179648; https://www.ncbi.nlm.nih.gov/ geo/, GSE175804.

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Ethical and Welfare Committee of Zhejiang Chinese Medical University.