Long noncoding RNAs associated with nonalcoholic fatty liver disease in a high cholesterol diet adult zebrafish model

The mechanism of nonalcoholic fatty liver disease (NAFLD) has not been completely revealed. In this study, we investigated the association of liver histological changes and long noncoding RNAs (lncRNAs) in the NAFLD zebrafish model. Forty zebrafish were fed a high-cholesterol diet (1.5 g per day) for 8 weeks. We measured fatty liver changes in the zebrafish liver using oil red O staining and divided them into two groups based on high and low scores. We pooled each group of zebrafish livers and identified lncRNAs, miRNAs, and mRNAs using Next-generation sequencing. Human homologs of lncRNAs were identified using ZFLNC, Ensembl, and NONCODE. We found several significant genes, including 32 lncRNAs, 5 miRNA genes, and 8 protein-coding genes, that were associated with liver metabolism and NAFLD-related functions in zebrafish. In particular, eight conserved human homologs of lncRNAs were found. We discovered the human homologs of eight lncRNA candidates from fatty liver zebrafish for the first time. The spectrum of biological mechanisms by which lncRNAs mediate their functional roles in NAFLD in a high cholesterol diet adult zebrafish model remains to be uncovered.


Scientific Reports
| (2021) 11:23005 | https://doi.org/10.1038/s41598-021-02455-0 www.nature.com/scientificreports/ lncRNAs are highly conserved due to their short lengths to preserve functional or secondary structures 13 . Although very few lncRNAs exhibit sequence conservation across species 14,15 , since lncRNAs contribute more to cell-line specificity than protein-coding genes proven by Djebali et al. 16 , meaningful research on lncRNAs associated with NAFLD is possible if one can find lncRNAs in certain animals that are homologous in humans. In this work, to identify lncRNAs in the NAFLD zebrafish model, we performed high-cholesterol diet feeding, measured and grouped NAFLD in the liver of zebrafish, and performed next-generation sequencing in the zebrafish liver ( Fig. 1). Using bioinformatics approaches, we found eight human homologs of lncRNA candidates that might be associated with liver metabolism and NAFLD.

Results
Nonalcoholic fatty liver disease in the adult zebrafish model fed a high-cholesterol diet. The fatty liver change scores of 36 adult zebrafish after 8 weeks of a high-cholesterol diet are summarized in Table 1. For example, 6 frozen sections (Fig. 2a) of a liver were obtained from zebrafish No. 5 and stained with oil red O. If the 10% fatty change in oil red O was the highest among the six slides, it became the representative score of zebrafish No. 5 (Fig. 2b,c). A representative section of oil red O stain from the zebrafish before feeding high cholesterol diet showed negative expression (Fig. 2d) (see the workflow diagram of Fig. 1). All zebrafish sections after the high cholesterol diet were stained with oil red O, divided into two groups: the low and the high fatty liver groups. The high group contained zebrafish that have proportion scores of no less than 50 (Fig. 2e). At higher magnification, in the representative higher score group, both macrovesicular and microvesicular fatty changes were severely detected compared to the negative internal control (Fig. 2f).
Long noncoding RNAs associated with nonalcoholic fatty liver disease in zebrafish. After the library preparation and validation quality check, the high fatty change group and low fatty change group were included in the final RNA sequencing. Through a next-generation sequencing approach, differentially expressed gene (DEG) analyses between high and low fatty change group using edgeR was conducted and 5065 selected genes satisfying the |fold change| ≥ 2 and p-value < 0.05 condition was extracted. The heatmap using the Z-score of the log 2 based Top 8 differentially expressed protein-coding genes with normalized value by Morpheus Soft-   (Fig. 3a,b). The smear plot and volcano plot between the high-and low-fatty change groups are featured in Fig. 3c,d. We selected 548 lncRNAs that had a greater than twofold change, either downregulated or upregulated, than the control group with statistical significance with a p-value of less than 0.05. Then, we chose 28 downregulated and 8 upregulated lncRNAs that had greater than tenfold changes than the control group with p-values of less than 0.05 (Fig. 3b). The identified 28 lncRNAs are listed in Table 2.
Identification of human homologs of eight long noncoding RNAs of zebrafish. To obtain the human homologs of each lncRNA of the NAFLD zebrafish, we performed ZFLNC BLAST and NCBI FASTA. First, in the NCBI section, we copied the sequence of the lncRNA and put it in the BLAST section of ZFLNC. The genes with the highest score were designated human homologs of the specific lncRNAs. Eight (ZFLNCT18125,    Table 2). The human homologs of eight lncRNAs according to ZFLNC are NONHSAT149958, ENST00000475761, NONHSAT180755, NONHSAT174445, NONHSAT205239, ENST00000415679, ENST00000426250, and NONHSAT030573 (Table 2). We searched these human loci in the Ensembl database (http:// www. ensem bl. org) and the NONCODE database (http:// www. nonco de. org) to obtain general information, sequences, expression profiles in human tissues, and exosome expression profiles. We put the FASTA link for each conserved lncRNAs in supplementary data 3. Since these genes are specific lncRNA candidates that were identified for the first time, there is limited information, especially regarding their structures. However, among the eight conserved lncR-NAs, six lncRNAs were downregulated and two lncRNAs were upregulated in the zebrafish liver with high fatty changes compared to that of low fatty changes, as demonstrated by Next-generation sequencing and ZFLNC. The human tissue expression profiles represented by FPKM/TPM (Fig. 4a) and exosome expression profiles represented by FPKM ( Fig. 4b) were evaluated for each conserved lncRNA.
Identification of miRNA genes and protein-coding genes in zebrafish fed a high-cholesterol diet using Next-generation sequencing. Through a next-generation sequencing approach, among the 5065 genes that were associated with NAFLD of zebrafish, we detected 5 miRNA genes and 8 protein coding genes that had greater than a 40-fold change than the lower score group, either downregulated or upregulated. Among them, two miRNA genes and 8 protein-coding genes showed a statistically significant difference from the lower cholesterol group (p-values less than 0.05). However, we could not find the human homologs of each miRNA and protein-coding gene. The identified miRNA and protein-coding gene are listed in Table 3. In contrast, we found that there was a lncRNA-mRNA (protein-coding gene) interaction network among 5 out of 8 conserved lncRNAs in the zebrafish using ZFLNC (http:// www. zflnc. org/) (Fig. 5).

Discussion
Given its clinical importance and enormous socioeconomic impacts, the underlying mechanism, as well as various clinical interventions, have been extensively studied for NAFLD. However, since the pathophysiological mechanism is very complex, drugs that directly act on NAFLD have yet to be developed. Traditional methods of developing NAFLD drugs include in vitro cell line experiments and in vivo mammalian models, such as mice 17 . In the case of mouse models, it takes a long time to observe pathological lesions in the liver through the standard high-fat diet in which mice start to have fatty liver changes from the 8th week to the 12th week and a nonalcoholic steatohepatitis feature at the 12th week 18 . In addition, it is difficult to proceed due to the high cost and low capacity of newly developed drugs. Moreover, in vitro cell line experiments are limited to multiplestage NAFLD 17 . Therefore, an in vivo model is needed that can proceed with a high-fat diet in a short period of www.nature.com/scientificreports/ time and can be used to screen the effectiveness of a drug with a small amount of NAFLD medicine. Recently, a high-cholesterol diet-induced larval zebrafish model (Danio rerio) was highlighted as a screening drug for the progression of NAFLD 19 . In addition, the model was used to perform mechanistic research. Therefore, in this work, we wanted to explore the association of fatty liver changes and lncRNAs in the NAFLD zebrafish model. We fed a high-cholesterol diet to wild-type adult zebrafish, measured the occurrence of NAFLD in the liver of zebrafish, and checked the lncRNAs associated with NAFLD outbreaks using Next-generation sequencing in the liver of zebrafish. Then, bioinformatics approaches, including ZFLNC (ZFLNC BLASTN 2.7.1 +), Ensembl, and NONCODE were undertaken as part of this study. We found eight conserved lncRNA candidates that might be associated with liver metabolism and NAFLD-related functions. ZFLNC is a database of vast amounts of information and knowledge about zebrafish lncRNAs. It is a well- From an evolutionary point of view, conserved sequences are similar or identical in RNA, DNA, or proteins across a species, within a genome, or in gene transfer between two organisms. Conservation, in the evolutionary sense, means that a sequence has been genetically preserved due to natural selection. From generation www.nature.com/scientificreports/ to generation, nucleic acid sequences can be gradually changed due to deletions and random mutations 21,22 .
In addition, when chromosomal rearrangement occurs, sequences may be recombined or deleted. Conserved sequences, unlike other sequences, are rarely mutated and preserved within the genome 23 . Either coding or noncoding sequences are conserved. Noncoding genes, which are essential for gene regulation, are conserved in a genome. In lncRNAs, however, sequence conservation generally does not occur as well as in protein-coding genes, and instead, sequences related to secondary structure or function are primarily conserved 24,25 . Nevertheless, the sequence conservation of lncRNAs is meaningful. While the protein-coding gene or amino acid sequence requires a strong functional restriction in the length or number of genes, lncRNAs, due to their relatively short structures of length, are highly conserved, preserving their functional or secondary structure 13 . Although very few lncRNAs exhibit sequence conservation across species, meaningful research on lncRNAs associated with certain diseases is possible if we can identify lncRNAs in certain animals that are homologous in humans. This method has the advantage of allowing experiments to be conducted in animals with a conserved gene that cannot be tested in humans for economic or ethical reasons. Herein, we conducted our study using bioinformatics approaches, including ZFLNC, Ensembl, and NONCODE. We found several significant genes that were associated with liver metabolism and NAFLD-related functions in zebrafish. Specifically, eight conserved lncRNAs were found to have homologs in humans ( Table 2).  www.nature.com/scientificreports/ Eight conserved lncRNAs are either downregulated (6 lncRNAs) or upregulated (2 lncRNAs) in the liver with high fatty changes compared to that in the liver with low fatty changes, as demonstrated by Next-generation sequencing. The 6 downregulated lncRNAs are expected to protect or prevent NAFLD, and the 2 upregulated lncRNAs are expected to cause fatty changes in the liver. The human tissue expression profiles represented by FPKM/TPM based on data from the human BodyMap (Fig. 4a) and exosome expression profiles represented by FPKM based on data from NCBI GEO (Fig. 4b) were evaluated for each conserved lncRNA. A few lncRNAs have human tissue expression profiles in adipose tissue (NONHSAT149958, and NONHSAT205239). It is not clear whether they contribute to the mechanisms affecting fatty changes in the liver with the information described above. In addition, there were two lncRNAs (NONHSAT205239, and NONHASTAT030573) with exosome expression profiles in HepG2 cells (hepatocellular carcinoma cell line exosomes). Since NONHSAT205239 has both adipose tissue expression and an exosome expression profile from the hepatocellular carcinoma cell line, it will be the most powerful candidate for revealing the mechanistic potential for NAFLD. In addition, we found that there was a lncRNA-mRNA (protein-coding gene) interaction network among 5 out of 8 conserved lncRNAs in the zebrafish (Fig. 5). Among the eight candidates, the information provided by ZFLNC (http:// www. zflnc. org/) and zfin (https:// zfin. org/) was identified, five (ZFLNCG11695, ZFLNCG07435, ZFLNCG06796, ZFLNCG10246, and ZFLNCG03169) of which were identified as having related protein-coding genes of zebrafish, and three (ZFLNCG07435, ZFLNCG10246, and ZFLNCG03169) of which were identified as having genes (prkcha, ankha, si:dkey-226m8.9, kcnjla.2, trpv6, kcnjla.1, and znf148) related to human disorder. Evaluating the regulatory interactions between those genes may help to identify biological mechanisms and treat NAFLD.
To date, the effects and mechanisms of lncRNA-mediated regulation in the pathogenesis of a particular disease have not been properly identified, especially the effects of lncRNAs on NAFLD. Recently, a regulatory interaction between a conserved pair of miRNAs and lncRNAs in zebrafish was discovered 13 . In addition, coding RNA-miRNA-lncRNA network and ceRNA network have been analyzed to discover the effect of tumor microenvironments in 1p/19q codeletion in oligodendrogliomas 26 . In this study, we discovered human homologs of eight lncRNA candidates from the fatty liver of zebrafish for the first time. Although we could not evaluate all of the lncRNAs, up to 547 genes identified by Next-generation sequencing, we identified human homologs of eight out of 32 lncRNAs, between high and low-fatty zebrafish liver groups. Since information on zebrafish lncRNA has not yet been properly organized, there is a limit to building a ceRNA network or finding the regulatory interaction of miRNA in a specific conserved lncRNA. In this study, the zebrafish gene was used to identify conserved lncRNAs for NAFLD mechanical research, and the detailed and practical methods for NAFLD animal model experiments were elucidated. Using a high cholesterol diet adult zebrafish model, we expect to uncover the biological mechanism and treatment by which lncRNA candidates mediate their functional roles in NAFLD. Among the eight candidates, five (ZFLNCG11695, ZFLNCG07435, ZFLNCG06796, ZFLNCG10246, and ZFLNCG03169) of which were identified as having related protein-coding genes of zebrafish, and three (ZFLNCG07435, ZFLNCG10246, and ZFLNCG03169) of which were identified as having genes (prkcha, ankha, si:dkey-226m8.9, kcnjla.2, trpv6, kcnjla.1, and znf148) related to human disorder. Information provided by ZFLNC (http:// www. zflnc. org/) and zfin (https:// zfin. org/). Measurement of fatty changes in the liver using oil red O staining. We used frozen sections of each zebrafish (Fig. 6a) and cut them into longitudinal sections for fat staining. When the liver that surrounds the stomach, the largest organ in the zebrafish, began to appear in the cut surface (Fig. 6b) of the frozen section, it was cut into serial sections with a thickness of 4 µm (Fig. 6c), and oil red O staining was performed on the best section. We obtained six or more frozen liver sections from each zebrafish and performed oil red O staining according to the modified method of our previous study 15 . Sections were cut at a thickness of 4 µm. In addition, 80% propylene glycol was added for 2 min, and oil red O staining was performed for 35 min. After washing the dyed slides twice with distilled water, further washing was carried out microscopically to prevent overstaining. Otherwise, the method was the same as in our previous study 15 . The proportion of fatty changes on a slide with oil red O staining was evaluated by two pathologists through a blinded method. The highest proportional value among the six sections of each zebrafish was the representative score.
RNA extraction and Next-generation sequencing analysis. Each piece of liver tissue was obtained from the remaining frozen specimen using a scalpel. Then, those tissues were immediately crushed in 10 µl of quiazole and stored at -70 °C. Liver specimens of zebrafish numbers 4, 15, 22, 25, 32, 34, and 36 with high fatty liver scores based on oil red O staining were pooled into a "high" group. Zebrafish numbers 1, 5, 11, 13, 14, 20, and 29 were pooled into a "low" group. For RNA-sequencing, we extracted RNA from two pooled groups. The method of RNA extraction from the tissue is described in Supplementary data 4. We performed library preparation and validation quality check before the Next-generation sequencing. Rebo-depletion was performed by Ribo-zero H/M/R Gold, during the library preparation. Total RNA integrity was checked using an Agilent Technologies 2100 Bioanalyzer with an RNA Integrity Number (RIN) value greater than or equal to 7. To verify the size of PCR enriched fragments, template size distribution was checked running on an Agilent Technologies 2100 Bioanalyzer using a DNA 1000 chip. Transcriptome sequencing was performed using an Illumina platform. The raw data of sequencing were extracted as fragments per kilobase of exon per million fragments mapped (FPKM) across each sample. The statistical significance in the fold change of transcript expression profile was determined by paired t-tests (Macrogen, https:// www. macro gen. co. kr). Data of the transcriptome sequencing and gene ontology analyses are summarized in Supplementary data 5 and 2, respectively.