Abstract
Background/Aims: To analyze the long noncoding (lncRNA)-mRNA expression network and potential roles in rat hepatic stellate cells (HSCs) during activation. Methods: LncRNA expression was analyzed in quiescent and culture-activated HSCs by RNA sequencing, and differentially expressed lncRNAs verified by quantitative reverse transcription polymerase chain reaction (qRT-PCR) were subjected to bioinformatics analysis. In vivo analyses of differential lncRNA-mRNA expression were performed on a rat model of liver fibrosis. Results: We identified upregulation of 12 lncRNAs and 155 mRNAs and downregulation of 12 lncRNAs and 374 mRNAs in activated HSCs. Additionally, we identified the differential expression of upregulated lncRNAs (NONRATT012636.2, NONRATT016788.2, and NONRATT021402.2) and downregulated lncRNAs (NONRATT007863.2, NONRATT019720.2, and NONRATT024061.2) in activated HSCs relative to levels observed in quiescent HSCs, and Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analyses showed that changes in lncRNAs associated with HSC activation revealed 11 significantly enriched pathways according to their predicted targets. Moreover, based on the predicted co-expression network, the relative dynamic levels of NONRATT013819.2 and lysyl oxidase (Lox) were compared during HSC activation both in vitro and in vivo. Our results confirmed the upregulation of lncRNA NONRATT013819.2 and Lox mRNA associated with the extracellular matrix (ECM)-related signaling pathway in HSCs and fibrotic livers. Conclusion: Our results detailing a dysregulated lncRNA-mRNA network might provide new treatment strategies for hepatic fibrosis based on findings indicating potentially critical roles for NONRATT013819.2 and Lox in ECM remodeling during HSC activation.
Introduction
Fibrosis is characterized by an excessive deposition of collagen due to an imbalance between collagen synthesis and degradation [1]. Hepatic fibrogenesis results from the accumulation of extracellular matrix (ECM), including collagen, proteoglycan, and adhesive glycoproteins, which is mainly produced by hepatic stellate cells (HSCs). HSCs are types of mesenchymal cells located in the perisinusoidal space (space of Disse) beneath the endothelial barrier [2]. Upon activation by fibrogenic stimuli, HSC morphology and physiology change significantly during myofibroblastic transdifferentiation, which is regulated by multiple cellular signal-transduction pathways and cytokines [3]. Among these, Lox and the Lox-mediated ECM-related pathway are probably the most important [4]. As the key enzyme in liver fibrosis, Lox controls collagen and elastin maturation, promoting ECM cross-linking and participating in the initiation and/or progression of hepatic fibrogenesis [5].
Previous studies attempted to reveal the mechanisms underlying HSC activation mainly at the mRNA level [6]; however, a comprehensive understanding of fibrogenesis remains beyond our reach due to its complexity, especially involving the intricate regulation of gene expression. Fortunately, recent progress in RNA-Seq technology focusing on long noncoding RNAs (lncRNAs) and bioinformatics analysis added new depth to our knowledge of gene regulation during the process of HSC transformation [7].
To better understand the intricate regulatory network associated with HSC activation, we performed RNA sequencing and comparative bioinformatics analysis of quiescent and activated HSCs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to identify the functions of target mRNAs exhibiting differential expression levels correlating with changes in lncRNA expression. Additionally, we established co-expression networks and potential mRNA-targeting relationships and also verified changes in NONRATT013819.2 and lysyl oxidase (Lox) expression during HSC activation both in vitro and in vivo. Moreover, we discussed the possible lncRNA-mRNA interactions associated with fibrogenesis and enhanced our knowledge regarding the mechanisms associated with liver disease.
Materials and Methods
Isolation, culture, and identification of rat HSCs
HSCs were isolated from three normal male Sprague-Dawley rats (weight: 300–400 g; Shanghai Laboratory Animal Center of Chinese Academy of Sciences, Shanghai, China) by in situ perfusion and density gradient centrifugation [8]. All materials for HSC isolation were obtained from commercial sources as previously described [9]. The rats received humane care according to the Guide for the Care and Use of Laboratory Animals of the Chinese Academy of Sciences.
Isolated rat HSCs were cultured in Dulbecco’s modified Eagle medium supplemented with 10% fetal bovine serum, 100 U/mL penicillin, 100 mg/mL streptomycin, and 2 mmol/L glutamine. Thereafter, quiescent and activated HSCs were harvested at days 2 and 14 for experimentation.
Cell purity was verified by immunohistochemical examination using monoclonal antibodies specific for desmin (1:100; Sigma-Aldrich, St. Louis, MO, USA). Immunohistochemical staining for desmin was performed as previously described [10]. Nuclei were labeled with Hoechst 33258 (Roche, Basel, Switzerland).
Establishment of an animal model of liver fibrosis and histological examination
Thirty-six Sprague-Dawley rats were divided into three groups (normal, control, and fibrosis model). Fibrosis-model rats were injected subcutaneously with 40% CCl4 (3 mL/kg CCl4:olive oil = 2:3) every 3 days for 10 weeks [11]. Control rats received only olive oil. Model rats were sacrificed at 4, 6, 8, or 10 weeks, respectively, and the degree of liver fibrosis was evaluated by microscopy.
Intrahepatic HSCs were isolated from CCl4-induced fibrotic rat livers as previously described. Liver tissues were fixed with 4% formaldehyde in phosphate-buffered saline and embedded in paraffin, and 5-μm-thick sections were prepared. All sections were stained with hematoxylin and eosin and underwent standard Van Gieson staining. Total RNA from rat HSCs was extracted, and fibrosis was graded according to the Ishak modified-staging system [12]. Histopathology was interpreted by two independent board-certified pathologists who were blind to the study.
RNA isolation, library preparation, and sequencing analysis
Total RNA was extracted from HSCs using TRIzol (Invitrogen, Carlsbad, CA, USA) and purified according to manufacturer instructions. Total RNA was quantified using a NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity was determined using an Agilent 2100 system and an RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, USA). Ribosomal RNA was removed using an EpicentreRibo-zero rRNA removal kit (Thermo Fisher Scientific) to enable accurate lncRNA analysis. RNA-Seq was conducted on an Illumina genome analyzer IIx from Shanghai Biotechnology Corporation (Shanghai, China) using next-generation sequencing analysis. Genome mapping was performed using TopHat software (v2.0.9; https://ccb.jhu.edu/software/tophat/index.shtml).
We determined differentially expressed genes and lncRNAs using edgeR analysis [13] for gene-expression profiling (http://bioconductor.org/packages/release/bioc/html/edgeR.html) and the Cufflinks algorithm [14] for identification of transcripts from RNA-Seq data. False-discovery rate (FDR) was used to determine the significance threshold of the p-value for multiple tests. A fold change ≥2 and an FDR ≤0.05 represented the thresholds signifying genes that were considered as being differentially expressed.
LncRNA expression, identification, and classification
Raw RNA-Seq samples were obtained from sequencing Poly-A+ fractions using paired-end Illumina reads (Illumina, San Diego, CA, USA). The mapped reads of each sample were assembled using the Cufflinks. Application function of the Cuffdiff algorithm (http://cufflinks.cbcb.umd.edu/manual.html#cuffdiff) to determine differences between samples for lncRNA analysis. Differences in lncRNA-screening criteria were determined by a q-value ≤0.05 and a fold-change in expression ≥2.
Cuffcompare [15] was used to compare candidate sequences with known lncRNAs, and we used the following process to select new lncRNAs: 1) transcription length ≥ 200 bp and number of exons ≥2; 2) a predicted open reading frame length <300 bp; 3) intersecting prediction results from the protein family database [16] (Pfam; http://www.uniprot.org/database/DB-0073), the coding potential calculator [17] (CPC; http://cpc.cbi.pku.edu.cn/), and the coding-noncoding index [18] (CNCI; https://omictools.com/coding-non-coding-index-tool) selected based on a CPC score <0, a CNCI score <0, and Pfam for determination of non-significant transcripts associated with potential lncRNAs; and 4) sequences matching those of known lncRNAs were removed.
According to their locations relative to the nearest protein-coding genes, the annotated lncRNAs were subdivided into intergenic lncRNAs (lincRNAs), intronic-antisense lncRNAs, exonic-sense lncRNAs, bidirectional lncRNAs, or exonic-antisense and intronic-sense lncRNAs[19].
Feature comparisons of transcripts in HSCs
To analyze the differences between IncRNAs and protein-coding genes, the distributions of lncRNA- and mRNA-transcript length and exon number were analyzed in HSCs. To compare transcript-expression levels and verify the consistency of predicted lncRNA with general characteristics, lncRNA- and mRNA-transcript levels were averaged, and a box diagram was created using log10 [fragments per kilobase of transcript per million (FPKM) + 1] values.
Target-gene prediction
To determine lncRNA function, we performed transcriptome analysis by lncRNA sequencing (lncRNA-Seq) to predict potential lncRNA targets. First, we predicted the cis- and trans-target genes of lncRNAs, with the trans-target genes predicted based on the rat mRNA database (http://rgd.mcw.edu/). Complementary or similar sequences were selected using BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi), and similar free energy values between sequences were calculated using RNAplex [20] (www.bioinf.uni-leipzig. de/Software/RNAplex/), with similarities above a specified threshold selected. For the cis-target genes, lncRNAs >10 kb in length were selected.
GO analysis and KEGG pathway annotation
The functions of predicted lncRNA targets were studied via bioinformatics analysis. The DAVID program (http://david.abcc.ncifcrf.gov/) was used to perform GO analysis to identify the molecular functions of potential target genes. Additionally, KEGG (http://www.genome.ad.jp/kegg/) and BioCarta (http://www. biocarta.com) were used to analyze the potential functions of these genes in given pathways. Significant correlations between target genes and their associated functions and/or pathways were assessed based on a threshold of p < 0.05.
Analysis of lncRNA-mRNA regulatory networks
Co-expression-network analysis was performed to determine the correlated expression of genes in order to construct weighted co-expression networks [21]. To determine relationships between lncRNA expression and regulation of mRNA co-expression, Pearson’s correlation coefficients were calculated [22], and the R value was used to compare the Pearson’s correlation coefficients between lncRNAs and differentially expressed mRNAs, with relationships exhibiting coefficients ≥0.98 used to construct the network via Cytoscape (http://www.cytoscape.org/). Co-expression networks of possible lncRNA-mRNA interactions representing critical lncRNAs and target mRNAs were established accordingly.
Quantitative reverse transcription polymerase chain reaction (qRT-PCR) validation
RNAs were extracted according to manufacturer instructions, and RT reactions were performed according instructions associated with the ReverTra Ace kit (Toyobo, Osaka, Japan). The expression of selected upregulated lncRNAs (NONRATT012636.2, NONRATT016788.2, and NONRATT021402.2), downregulated lncRNAs (NONRATT007863.2, NONRATT019720.2, and NONRATT024061.2), and a lncRNA-mRNA pair (NONRATT013819.2 and Lox) were analyzed by qRT-PCR using a SYBR Green PCR kit (TaKaRa, Shiga, Japan). Glyceraldehyde 3-phosphate dehydrogenase (Gapdh) mRNA was used as an internal control. To determine relationships between lncRNA-mRNA pairs (NONRATT013819.2 and Lox), we also calculated correlation coefficients (using Pearson’s correlation). For quantitative results, expression of each lncRNA was represented as a fold change using the 2-ΔΔCt method, followed by statistical analysis. Primer sequences are listed in Table 1.
Statistical analyses
All data were expressed as the mean ± standard deviation. Statistical analysis was performed using Student’s t test for comparison of two groups. A p < 0.05 was considered statistically significant, and the FDR was calculated to correct the p-value.
Results
In vitro activation of rat HSCs
HSCs (∼5 × 108/mL) were harvested from each rat, with the isolated cells exhibiting a round morphology and a quiescent phenotype at 48 h. Gradually, activation was achieved after a 7- to 14-day culture, with only cultures >95% viable (as determined by Trypan blue exclusion) used for experiments. The purity of the quiescent (>95%) and activated (>95%) HSCs was demonstrated by immunohistochemical examination to determine desmin levels (Fig. 1).
Differentially expressed lncRNAs and mRNAs during HSC activation
After discarding adaptor and low-quality sequences, we obtained raw data consisting of >84,970,970 clean reads (93.3–94.1%) from the Illumina HiSeq platform. TopHat was subsequently used for genome mapping of the clean reads, and Cufflinks was used to reconstruct the transcripts. After basic filtering, we obtained 16,816 tentative lncRNAs containing 3559 lincRNAs, 281 intronic-antisense lncRNAs, 9235 exonic-sense lncRNAs, 1503 bidirectional lncRNA, 281 exonic-antisense lncRNAs, and 2055 intronic-sense lncRNAs (Fig. 2). Following sequence verification, 1454 novel lncRNAs were identified for further analysis.
We then determined differentially expressed lncRNAs between quiescent and activated rat HSCs, identifying 24 lncRNAs meeting the selection threshold (fold change ≥2 and p ≤ 0.05). Among these, 12 lncRNAs were identified as upregulated by >2-fold in activated HSCs as compared with their levels in quiescent HSCs (Fig. 3). Additionally, we found that among 24 different lncRNAs, only two (ENSRNOT00000081905 and TCONS_00042359) did not predict the cis-target gene, with no trans-target genes predicted (Table 2).
To investigate possible lncRNA-mRNA interactions, we performed RNA-Seq analysis of the mRNA data, revealing 529 differentially expressed mRNAs between quiescent and activated rat HSCs. Among these, 155 mRNAs were upregulated in activated rat HSCs, whereas 374 mRNAs were downregulated (p < 0.01) (Fig. 3).
Distribution of transcript length and exon number between HSC lncRNAs and protein-coding genes
Previous studies showed that both plant and animal lncRNAs are shorter and harbor fewer exons than protein-coding genes. To determine whether HSC lncRNAs share these features, we compared the transcript length, exon number, and expression level of 16,816 lncRNAs with those of protein-coding genes present in HSCs. Figure 4 shows that >80% of lncRNAs ranged in size from 700 to 900 nucleotides, whereas the length of all mRNAs exceeded 1000 bp. Interestingly, most (90%) of the genes encoding HSC lncRNAs only contained one to three exons, whereas the number of exons in protein-coding genes ranged from three to 20. Compared to mRNA-expression levels, lncRNA expression was lower than that of protein coding genes (Fig. 4). These results indicated that the majority of HSC lncRNAs were shorter than HSC protein-coding genes, contained fewer exons, and exhibited lower expression levels.
GO analysis and KEGG pathway annotations
Bioinformatics analyses were used to analyze lncRNA function based on their target genes. GO analysis focusing on biological processes, molecular functions, and cellular components resulted in classification of the top 30 of 783 GO terms (Fig. 5). The most significant GO functions associated with potential lncRNA targets were related to the extracellular exosome, extracellular vesicle, and extracellular region. KEGG pathway analysis (Fig. 6) identified associations with 11 pathways (p < 0.05) ranked by target-gene count, with some highly enriched pathways related to ECM remodeling involving mucin-type O-glycan biosynthesis and glycosaminoglycan biosynthesis shown to participate in HSC activation. Moreover, Lox, a critical member of ECM-related pathways, is a significant target of the lncRNA NONRATT013819.2, both of which were revealed as differentially expressed in this study.
Co-expression relationships between lncRNAs and mRNAs
Networks describing the co-expression of genes and lncRNAs can provide a deeper understanding of the relationship between lncRNAs and regulation of mRNA levels. According to the similarities in gene-/lncRNA-expression profiles, we calculated Pearson’s correlation coefficients between lncRNA and gene expression and used expression relationships ≥0.98 to construct a co-expression network. Our network consisted of 632 matched lncRNA-mRNA pairs representing 24 differentially expressed lncRNAs and 122 differentially expressed mRNAs (Fig. 7), with 352 and 280 pairs representing positive and negative correlations, respectively. Our results showed that lncRNA NONRATT013819.2 exhibited the highest correlation with co-expressed mRNAs, indicating a potentially important role in regulating differential gene expression between activated and quiescent HSCs.
Validation of genes present in the co-expression network by qRT-PCR
Twenty four lncRNAs were randomly chosen for qRT-PCR validation of significantly differential expression patterns exhibited between activated and quiescent HSCs (p < 0.05). Additionally, genes exhibiting higher FPKM values were also selected for qRT-PCR analysis. As shown in Figure 8, qRT-PCR analysis confirmed that NONRATT007863.2, NONRATT019720.2, and NONRATT024061.2 expression was downregulated, and NONRATT012636.2, NONRATT016788.2, and NONRATT021402.2 expression was upregulated in activated HSCs as compared with their levels in quiescent HSCs, which agreed with RNA-Seq results.
Correlation analysis of NONRATT013819.2-Lox co-expression during HSC activation in vitro and in vivo
Because of similarities in NONRATT013819.2-Lox co-expression in the network, we analyzed their chromosomal location and promoter-driven expression. NONRATT013819.2 is an antisense lncRNA located on rat chromosome 18:47500329–47634127 and adjacent to the Lox gene located on rat chromosome 18:47500330–47577819 (http://www.bioinfo.org/noncode/) (Fig. 9). We observed upregulated expression of NONRATT013819.2 and Lox according to qRT-PCR analysis during HSC activation (on days 2, 7, and 14 and in the HSC-T6 rat HSC line) and relative to levels observed in quiescent HSCs, as well as significantly positive correlations in their co-expression patterns (Figs. 1 and 9).
Histological examination revealed that the degree of liver fibrosis progressed in rats that received CCL4 injections as compared with that observed in controls. We also observed that fatty tissue degeneration, necrosis, and fibrogenesis became more evident over time following CCl4 administration. Furthermore, intrahepatic HSCs isolated from rat livers and subjected to qRT-PCR analysis confirmed the significantly upregulated co-expression of NONRATT013819.2 and Lox associated with HSC activation status and fibrogenesis progression as compared with that observed in HSCs isolated from control rats (Fig. 10).
Discussion
HSCs represent the primary mesenchymal cells of the liver and play a major role in ECM maintenance and as inducers of hepatic fibrogenesis upon exposure of the liver to stress [23]. In such cases, HSCs exhibit altered morphology and physiology through the process of activation, which is accompanied by a series of pathological changes based on large-scale gene overexpression and/or repression [24-27]. However, limited understanding of gene regulation associated with these processes rendered these HSC-specific mechanisms vague and controversial.
LncRNA are small noncoding RNAs >200 nucleotides in length and intimately involved in various critical cell processes, such as development, proliferation, differentiation, and pluripotency [28]. Recent studies reported low expression levels of hepatocellular carcinoma lncRNA maternally expressed gene 3 in CCl4-induced mouse liver fibrosis and significantly lower expression in the human-shaped cell line LX-2 stimulated by transforming growth factor-β1 [29]. However, the expression and biological function of lncRNAs in HSCs remain largely unknown. Therefore, we performed RNA-Seq analysis of rat quiescent and activated HSCs, with results revealing significant differential expression of 24 lncRNAs and 529 mRNAs. Tests of NONRATT012636.2, NONRATT016788.2, NONRATT021402.2, NONRATT007863.2, NONRATT019720.2, NONRATT024061.2, NONRATT013819.2, and Lox further validated the reliability of RNA-Seq results.
To gain insight into the functions affected by lncRNA activity, GO and KEGG analyses were performed to elucidate the possible biological functions and mechanisms of lncRNAs in HSCs, resulting in the identification of 783 GO terms and 11 pathways. GO terms involving cellular components (extracellular vesicles, extracellular exosomes, and intracellular membrane-bound organelles) dominated the significant GO categories. Moreover, we noted that the significantly enriched terms derived from KEGG analysis agreed with the GO results. Among these pathways and predicted target genes, ECM-related pathways (mucin-type O-Glycan biosynthesis and glycosaminoglycan biosynthesis) and ECM protein cross-linking by Lox appear to play important roles in HSC activation. Integrated bioinformatics analysis of differentially expressed lncRNAs verified that NONRATT013819.2 was the only lncRNA targeting Lox, with its location adjacent to the Lox gene on the same rat chromosome. Furthermore, construction of a lncRNA-mRNA co-expression network displayed critical lncRNAs and mRNAs potentially involved in HSC activation. Additionally, co-expression analyses identified a subset of lncRNAs tightly linked to collagen genes and multiple proteins that regulate ECM remodeling during hepatic fibrogenesis. In this network, positively correlated co-expression of the lncRNA NONRATT013819.2 and Lox mRNA exhibited the highest level of enrichment in ECM-related pathways. Lox, also known as protein-lysine 6-oxidase, is rarely detected in normal or non-proliferative senescent HSCs, but is expressed in culture-activated HSCs and intra-hepatic HSCs following CCl4 injury [30, 31]. Recent studies attributed ECM remodeling, spontaneous reversion of fibrogenesis, and induced phenotypic reversion of HSCs to Lox overexpression [7, 32]. These findings suggested that the dysregulated lncRNA-mRNA network might be closely associated with ECM remodeling and HSC activation.
Our experimental findings were consistent with these reports. Subsequent examination of NONRATT013819.2 and Lox co-expression both in vitro and in vivo confirmed their upregulation during HSC activation relative to levels observed in quiescent HSCs and suggested their association with HSC activation and progression during fibrogenesis. Moreover, the most significant finding from data mining was the high degree of association between lncRNA NONRATT013819.2/Lox mRNA and the ECM, the production of which by HSC myofibroblast-like cells is the main cause of fibrogenesis and progression to liver failure in chronic liver diseases [33-35]. These results indicated that NONRATT013819.2-Lox co-expression and possible interaction might promote HSC transformation and proliferation by activating signaling pathways associated with ECM remodeling. However, the function of lncRNA in liver fibrosis and related mechanisms remains to be confirmed, given that a single lncRNA might potentially regulate the expression of multiple protein-coding genes at diverse levels [36]. Therefore, additional evidence regarding the cross-talk between NONRATT013819.2 and Lox remains to be experimentally validated in our future work.
Conclusion
In conclusion, we investigated the differential expression of lncRNAs and mRNAs during HSC transformation, with results of RNA-Seq analysis indicating that lncRNAs were involved in pathological processes associated with the progression of fibrogenesis. Furthermore, bioinformatics analysis and in vitro and in vivo experiments identified correlations between lncRNA NONRATT013819.2 and Lox mRNA expression and potential interaction related to ECM-remodeling pathways. These findings provide valuable insight into HSC activation and also facilitate access to novel therapeutic strategies for the treatment of liver fibrogenesis.
Acknowledgements
This work is supported by the National Natural Science Foundation of China (81570544 and 81100296, 81670546).
Disclosure Statement
All authors declare no conflicts of interest.
References
C.J. Guo and X. Xiao X. contributed equally to this work