{'fraction_non_alphanumeric': 0.036194051739308966, 'fraction_numerical': 0.023200563148941162, 'mean_word_length': 4.938860825640132, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Whole transcriptome analysis of RNA expression profiles reveals the potential regulating action of long noncoding RNA in lactating cows fed a high concentrate diet

Subacute ruminal acidosis (SARA) is a common metabolic disease in the dairy farming industry which is usually caused by an excessive amount of high concentrate diet. SARA not only threatens animal welfare but also leads to economic losses in the farming industry. The liver plays an important role in the distribution of nutritional substances and metabolism; however, a high concentrate diet can cause hepatic metabolic disorders and liver injury. Recently, noncoding RNA has been considered as a critical regulator of hepatic disease, however, its role in the bovine liver is limited. In this study, 12 mid-lactating dairy cows were randomly assigned to a control (CON) group (40% concentrate of dry matter, n = 6) and a SARA group (60% concentrate of dry matter, n = 6). After 21 d of treatment, all cows were sacrificed, and liver tissue samples were collected. Three dairy cows were randomly selected from the CON and SARA groups respectively to perform whole transcriptome analysis. More than 20,000 messenger RNA (mRNA), 10,000 long noncoding RNA (lncRNA), 3,500 circular RNA (circRNA) and 1,000 micro RNA (miRNA) were identified. Furthermore, 43 mRNA, 121 lncRNA and 3 miRNA were differentially expressed, whereas no obvious differentially expressed circRNA were detected between the 2 groups. Gene Ontology (GO) annotation revealed that the differentially expressed genes were mainly enriched in oxidoreductase activity, stress, metabolism, the immune response, cell apoptosis, and cell proliferation. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that the deferentially expressed genes were highly enriched in the phosphatidylinositol 3 kinase (PI3K)-serine/threonine kinase (AKT) signaling pathway (P < 0.05). According to KEGG pathway analysis, the differentially expressed lncRNA (DElncRNA) target genes were mainly related to proteasomes, peroxisomes, and the hypoxia-inducible factor-1 signaling pathway (P < 0.005). Further bioinformatics and integrative analyses revealed that the lncRNA were strongly correlated with mRNA; therefore, it is reasonable to speculate that lncRNA potentially play important roles in the liver dysfunction induced by SARA. Our study provides a valuable resource for future investigations on the mechanisms of SARA to facilitate an understanding of the importance of lncRNA, and offer functional RNA information.


Introduction
The dairy industry is gradually developing towards a largescale, intensive farming model in China with the increasing demand for dairy products. Although the production performance of dairy cows has improved gradually, it has become a great challenge for researchers and breeders to more effectively meet the nutritional needs of cows (Sundrum, 2015). In order to meet the nutritional requirement of cows, farmers have increased the usage of high concentrate diets (Patel, 2012). A high concentrate diet improves the production performance of dairy cows, however, it can also induce a variety of metabolic diseases (Kleen et al., 2003). Subacute ruminal acidosis (SARA) is one of the most representative metabolic diseases. SARA decreases rumen pH, which leads to microbiota disorders and chronic metabolic diseases (Ametaj et al., 2010). More importantly, SARA can cause a series of secondary diseases, increasing the elimination rate of cows (McCann et al., 2016). The incidence of SARA is increasing annually. Therefore, SARA not only poses a threat to the health of highyielding cows, but also results in great economic losses to the cattle industry (Kitkas et al., 2013).
RNA-seq-based transcriptome data provide a helpful platform for investigating the molecular aspects of metabolic diseases in animals. Therefore, we concentrated on the transcriptomic data to analyze the molecular aspects of SARA in the present study. Previous studies have suggested that eukaryotic genomes encode a large number of functional noncoding RNA (ncRNA) transcripts, including housekeeping and regulatory RNA (Huttenhofer et al., 2005;Brosnan and Voinnet, 2009). miRNA are the most widely known small/short ncRNA, which can regulate the expression of targeted genes (Gindin et al., 2015;Mogilyansky et al., 2016). circRNA is a new type of ncRNA that eliminates miRNA (Hansen et al., 2013;Memczak et al., 2013). Research in lncRNA has rapidly emerged and an increasing number of long noncoding RNA (lncRNA) have been reported in cows (Sun et al., 2016). Moreover, previous studies have revealed that the regulatory function of lncRNA at the transcriptional and post-transcriptional levels is related to disease (Gutschner and Diederichs, 2012;He et al., 2015). However, limited research has been conducted on ncRNA in dairy cows. Studies have indicated that the hepatic metabolic disorders caused by a high concentrate diet occur at the transcriptional level (Dong et al., 2013;Xu et al., 2015). Therefore, we speculated that ncRNA are responsible for hepatic metabolic abnormalities in dairy cows. Understanding the pathogenic mechanisms of a high concentrate diet by studying ncRNA would be helpful.
The liver is an important metabolic organ responsible for glucose and lipid distribution (Knegsel et al., 2005;Xu et al., 2015). In ruminants, the liver is the major site for gluconeogenesis which is the main way to maintain an adequate glucose supply in mammary tissue; and it is also the major site for fat synthesis which is important for the physical properties and quality of milk (Knegsel et al., 2005;Zhao and Keating, 2007;Dong et al., 2013). Therefore, it is of great significance to study the molecular mechanism of SARA by exploring the hepatic transcriptome in dairy cows. Our study used RNA-seq technology for the first time to conduct a joint analysis of mRNA and ncRNA to provide useful information for further research in liver metabolic disorders of dairy cows with SARA.

Ethics statement
In this study, all animal experiments were approved by the Institutional Animal Care and Use Committee of Nanjing Agricultural University. All procedures of the Guidelines on Ethical Treatment of Experimental Animals (2006) No. 398 set by the Ministry of Science and Technology, China and the Regulation regarding the Management and Treatment of Experimental Animals' ' (2008) No. 45 set by the Jiangsu Provincial People's Government, was strictly followed during the slaughter and sampling procedures.

Animals, diets and experimental design
Twelve healthy Chinese Holstein cows (mean body weight, 651 ± 54 kg) during mid-lactation (mean milk yeild, 17.43 ± 4.04 kg/ d) were raised in dairy farm (Taizhou, Jiangsu Province, China) for the study. All the cows were in their second parity or third parity when the liver tissues were collected. These cows were kept in freestall housing in the same cubicle partition of the barn during the experimental period. After a week of prefeeding, they were randomly assigned to the CON or SARA group. The CON group received a low concentrate diet (40% of dry matter) (n ¼ 6), and the SARA group received a high concentrate diet (60% of dry matter) (n ¼ 6) for 3 weeks. The experimental diets were formulated according to the lactation nutritional requests of a 30 kg milk/d producing cow. The ingredients in the diets and the nutritional composition are shown in Appendix Table 1. Crude protein, starch, ash and other extracts were measured according to the AOAC methods (Helrich and Helrich, 1990). Acid detergent fiber (ADF) and neutral detergent fiber (NDF) were measured by following the methods of Van Soest et al. (Van Soest et al., 1991). The amounts of NE L were calculated using NRC (2001). Ca and P were assayed by inductively coupled plasma spectrometry. The NFC content of the grain was determined by calculation. Cows were given free access to both food and water throughout the experimental time. Dry matter intake (DMI) (kg/d) was recorded daily (CON ¼ 23.72, SARA ¼ 20.82, SEM ¼ 0.68, P ¼ 0.639).

Sample collection and illumina deep sequencing
After fasting, all cows were sacrificed using electric stunning slaughtering. Samples of about 40 mL of the ruminal fluid were immediately collected through spontaneous efflux after aspirating with a modified 60-mL syringe, and strained through 4 layers of cheesecloth. The ruminal fluid pH was measured directly using a calibrated pH meter (Starter 300; Ohaus Instruments Co. Ltd., Nanjing, China). The liver samples were collected from the same sites and stored under À80 C for further RNA extraction.
Three dairy cows were randomly selected from the CON and SARA groups respectively. Total RNA was extracted from 50 mg of liver tissue using total RNA extraction reagent following the manufacturer's protocol (Sangon, Shanghai, China). RNA degradation and contamination was monitored on 1.5% agarose gels. RNA concentration and purity was measured using advanced molecular biology equipment to ensure the high quality of samples for transcriptome sequencing.
1.5 mg RNA per sample was used for rRNA removal using the Ribo-Zero rRNA removal kit (Epicentre, USA). Sequencing libraries were generated by an ultra directional RNA library preparation kit for Illumina RNA library preparation (NEB, USA), following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Divalent cations were used to carry out fragmentation under elevated temperature in first strand synthesis reaction buffer (NEB, USA). Random hexamer primer and reverse transcriptase were used to synthesize first strand cDNA. Then second strand cDNA was synthesized by DNA polymerase I and RNase H. Through exonuclease/polymerase activities, the remaining overhangs were converted into blunt ends. After adenylation of 3 0 ends of DNA fragments, NEBNext adaptors with hairpin loop structures were ligated to prepare for hybridization. For preferentially selecting insert fragments of 150 to 200 bp in length, AMPure XP beads (Beckman, USA) were used to purify the library fragments.
Then 3 mL USER enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 C for 15 min before PCR. Next, highfidelity DNA polymerase, universal PCR primers and index primer were used to perform PCR (NEB, USA). Lastly, PCR products were purified and library quality was assessed on the Agilent 2100 bioanalyzer (Agilent, USA).
The clustering of the index-coded samples was performed on a cBot cluster generation system (Illumia, USA) according to the manufacturer's instructions. Then the library preparations were sequenced on an Illumina Hiseq platform and paired-end reads were generated.

Quality control
First of all, raw data of fastq format were cleaned by selfdeveloped Perl scripts. Through removing ploy-N-containing, low-quality reads, and adapter-containing reads from the raw data, the clean data were acquired. The Q20, Q30 and GC contents of the clean data were calculated. All the subsequent analyses were based on the high-quality clean data.

ncRNA identity
For lncRNA identity, the transcriptome was assembled using the StringTie based on the reads mapped to the reference genome. The assembled transcripts were annotated using the gffcompare program. CPC2/CNCI/Pfam/CPAT were combined to sort nonprotein coding RNA candidates. Putative protein-coding RNA were filtered out using a minimum length and exon number threshold. Transcripts with lengths more than 200 nt and 2 exons were selected as lncRNA candidates. Then, the protein-coding genes could be distinguished from the noncoding genes by CPC2/CNCI/Pfam/CPAT. Furthermore, cuffcompare was used to select the different types of lncRNA including lincRNA, intronic lncRNA, anti-sense lncRNA, and sense lncRNA; For miRNA identity, we compared the known reads to the miRBase (v22) database, allowing at most one mismatch, so that the identified reads were considered to be the known miRNA. For the sequences not identified to known miRNA, we used miR-Deep2 software to predict novel miRNA; for circRNA identity, we used circRNA identifier (CIRI) tools. During the first scanning of SAM alignment, CIRI detected junction reads with PCC signals that reflect a circRNA candidate. Preliminary filtering was implemented using paired-end mapping (PEM) and GT-AG splicing signals for the junctions. After clustering junction reads and recording each circRNA candidate, CIRI scaned the SAM alignment again to detect additional junction reads. At the same time, CIRI performed further filtering to eliminate false positive candidates produced from incorrectly mapped reads of homologous genes or repetitive sequences.

Quantification of gene expression levels
FPKM of both mRNA and ncRNA were calculated using StringTie (1.3.1) in each sample. FPKM were calculated through summing the FPKM of transcripts in each group. FPKM means fragments per kilobase of exon per million fragments mapped, calculated based on the length of the fragments and reads count mapped to this fragment.

Differential expression analysis
Analysis of 2 deferentially expressed conditions/groups was implemented using the DESeq R package (1.10.1). The resulting Pvalues were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). RNA absolute value of FDR < 0.05 and |(fold change)| ! 2 found by DESeq were assigned as deferentially expressed.

Function annotations of RNA in bovine liver
Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KEGG (Kyoto Encyclopedia of Genes and Genomes); GO (Gene Ontology). GO enrichment analysis of the differentially expressed genes (DEG) was implemented by the topGO R packages. We used KOBAS (Mao et al., 2005) software to test the statistical enrichment of differential expression genes in KEGG pathways. The sequences of the DEG was blast (blastx) to the genome of a related species (the protein interaction of which exists in the STRING database:http:// stringdb.org/) to get the predicted proteineprotein interaction (PPI) of these DEG. Then the PPI of these DEG were visualized in Cytoscape (Shannon et al., 2003).

Validation of RNA-seq results analysed by quantitative PCR (qPCR)
In order to validate the gene expression data obtained from RNA sequencing, 11 DEG and 12 differentially expressed lncRNA were radomly selected to perform qPCR. Two micrograms of total RNA was treated with RNase-free DNase (M6101; Promega, Madison, WI, USA) and reverse transcribed according to manufacturer's instructions. 2 mL of diluted cDNA (1:40, vol/vol) was used for qPCR, which was performed in an Mx3000 P (Stratagene, USA). SYBR green qPCR master mix was used for qPCR analysis (Vazyme, Nanjing, China). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH), which is not affected by the experimental factors, was chosen as the reference gene. All the primers, listed in Appendix Table 2, were synthesized by Tsingke Company (Nanjing, China). The relative gene expression level was calculated using the 2 ÀDDCt method. The data of gene expression are shown as the mean ± standard error. Statistical significance of the results was tested using the paired t-test.

Plasma triglyceride and glucose detection
During slaughter, the jugular vein was cut and blood samples were collected using 5 mL Vacutainer tubes (Jiangsu Kangjian Medical Apparatus Co., Ltd., Nanjing, China) with sodium heparin as an anticoagulant. Plasma was isolated from the blood samples by centrifugation at 3,000 Â g at 4 C for 15 min and stored at À 20 C. Plasma glucose and triglyceride were measured using an automatic biochemical analyzer (7020, HITACHI, Tokyo, Japan).
In this study, a total of 22,103 mRNA were identified. To identify ncRNA in bovine livers, transcripts were filtered following a rigorous set of criteria. Finally 10,376 lncRNA, 1,290 miRNA and 3,508 circRNA were obtained (Fig. 1A). To assess the distribution of the putative RNA, we analyzed the RNA expression ratio across the chromosomes. Description of Fig. 1B, the putative RNA were not equally distributed across the bovine chromosomes. Most mRNA and lncRNs were found on chromosome 10, whereas most miRNA were found on chromosome 4. In addition, most circRNA were found on chromosome 10 and X.
The length distribution of mRNA and ncRNA showed that mRNA and lncRNA have similar length distribution within 3,000 nt. As described in Fig. 2A and B, the length of most mRNA and lncRNA ranged from 200 to 400 nt. With the increase of the length, the number of mRNA and lncRNA gradually reduced. Different from mRNA and lncRNA, the length of most circRNA ranged from 400 to 600 nt, followed by 200 to 400 nt. As the length exceeded 600 nt, the number of circRNA decreased gradually (Fig. 2D). The final generated mature miRNA length mainly ranged from 20 to 24 nt, especially those enriched in 22 nt (Fig. 2C).

Analysis of differentially expressed RNA
As shown in Fig. 3A, under the requirement of FDR < 0.05 and |(fold change)| ! 2, we acquired 167 deferentially expressed RNA, including 43 mRNA, 121 lncRNA, and 3 miRNA. However, no deferentially expressed circRNA was found. The 121 DElncRNA consisted of 70 down-regulated and 51 up-regulated lncRNA transcripts. Among the 43 DEG, 12 mRNA up-regulated whereas 31 mRNA down-regulated. Moreover, 3 DEmiRNA were up-regulated compared to the CON group. The distribution of DEG and DElncRNA across the chromosomes showed that most DEG were located on chromosome 15 and 18, and most DElncRNA were located on chromosome 6 and 19 (Fig. 3B).

Analysis of ncRNA target genes
In order to find the relationship between DEG and ncRNA, 2 methods, position (cis-acting) and correlation between lncRNA and mRNA (trans-acting), were used to predict the lncRNA target genes. The results are shown in Fig. 4A and B, respectivelly. A total of 7,451 lncRNA and 2,826 miRNA target genes were found. All 42 DEG were the lncRNA target genes and they were trans-acted by lncRNA. Six DEG were the miRNA target genes and 5 DEG were both lncRNA and miRNA target genes. These results indicated that DEG tend to be directly trans-acted by lncRNA.

Functional annotation and classification of DEG
A total of 39 DEG were annotated on the GO database. In the biological process ontology, cellular process, single-organism, and biological regulation were the most abundant terms. In the cellular component ontology, cell part, cell, and organelle were the most abundant terms. In the molecular function ontology, binding, catalytic activity, and transporter activity were the most abundant terms. Furthermore, the biological adhesion, extracellular matrix part, collagen trimer and channel regulator activity were significantly different in proportion between the DEG and all genes (Fig. 5A). A total of 24 DEG were annotated to KEGG pathway database. We found that most annotated genes were involved in phosphatidylinositol 3 kinase (PI3K)-serine/threonine kinase (AKT) signaling pathway. Moreover, KEGG enrichment analysis showed that DEG were highly enriched in PI3K-AKT signaling pathway (P < 0.05) (Fig. 5B).

Functional annotation and classification of DElncRNA target genes
To infer the biological function of the identified DElncRNA, genes trans-acted by DElncRNA were annotated to the GO and KEGG databases. GO enrichment analysis showed that RNAdependent DNA biosynthetic process, transmembrane transport and polyamine catabolic process were the most enriched terms. RNA-dependent DNA biosynthetic process and oxidation-reduction  process were the most abundant terms (P < 0.001) (Fig. 6A). KEGG enrichment analysis showed that proteasome, peroxisome, and hypoxia inducible factor-1 (HIF-1) signaling pathway were the most enriched pathways (P < 0.005) (Fig. 6B).

Detailed information of selected mRNA and lncRNA
In order to find out potential mRNA and lncRNA involved in SARA, we screened 20 DEG and 60 DElncRNA with the minimum FDR for further investigation. Except for unannotated new genes, 20 DEG sorted by FDR are shown in Table 1. Among these DEG, only 4 DEG up-regulated and the others were down-regulated. GO annotation showed that AKR1B10 and DDH3 were related to the oxidation-reduction process; 6 DEG (AKR1B10, ACSM3, SGK1, CHI3L2, KRTAP5-AS1 and CCDC80) were related to substance metabolism including lipid, carbohydrate and protein metabolism; 5 DEG (C1orf198, SGK1, C21orf5, NEDD9 and PLEKHF1) were related to response to stimulus; 4 DEG (AKR1B10, SGK1, ZBTB16 and PLEKHF1) were related to cell proliferation or apoptotic process; 3 DEG (FERMT1, SGK1 and CCDC80) were related to cell adhesion; and 5 DEG (A2M, SGK1, ICAM3, CHI3L2 and HAMP) were related to inflammation or immune defense. Furthermore, the DElncRNA which have at least 3 target DEG were sorted by FDR, and 20 DElncRNA with the minimum FDR are shown in Table 2. From Table 2, we found LYPD1, DDH3, and ACSM3 were the most annotated genes. In addition, we analyzed the distribution of these DElncRNA on chromosomes. There were 3 DElncRNA on chromosome 5, 9, and 19, 2 DElncRNA on chromosome 1 and 15.

Validation of RNA-seq results by qPCR
Eleven DEG and 12 DElncRNA were subjected to qPCR analysis to validate the RNA-seq result. As is shown in Fig. 7A and C, the expression profiling of these DEG and lncRNA using qPCR had similar trends with the RNA-seq samples. The Pearson's correlation coefficients of the differential expression ratios between RNA-seq and qPCR were 0.9738 (P < 0.001) and 0.8192 (P ¼ 0.001), respectivelly ( Fig. 7B and D). It indicated that the expression differences of mRNA and lncRNA observed in RNA-seq between the CON and SARA groups were highly reliable.

Detection of plasma glucose, triglyceride and ruminal pH
Plasma glucose showed no significant difference between the CON and SARA groups (P > 0.05). However, the SARA group had a lower plasma triglyceride concentration than the CON group (P < 0.05). Compared with the CON group, the ruminal pH was lower in the SARA group (P < 0.01) ( Table 3).

Oxidative stress increased in SARA group
The antioxidant and lipid peroxide profiles in the livers of the dairy cows are presented in Fig. 8. The SOD levels were not changed (P > 0.05) by a high concentrate diet (Fig. 8B). However, CAT and TAOC dramatically decreased in the SARA group (P < 0.01) (Fig. 8A and D). In addition, there was an increase in MDA content in the SARA group (P < 0.05) (Fig. 8C).

Discussion
A high concentrate diet improves the performance of dairy cows over the short-term. However, long-term usage of a high concentrate diet could result in a series of nutritional and metabolic diseases (Jouany, 2006). SARA is one of these representative diseases, which leads to reduced dry matter intake and milk production, and increased elimination rate of cows (McCann et al., 2016). SARA has received increasing attention in recent years because of its damaging effects on animal welfare and profits. Therefore, investigating the occurrence and mechanisms of SARA is of great importance. Most studies attribute SARA to abnormal fermentation, volatile fatty acid accumulation, and reduced chewing frequency induced by a high concentrate diet (Yun and Han, 1989;Shaani et al., 2017). The liver plays important roles in nutrition metabolism and distribution and is key for maintaining milk fat, lactose, and lactoprotein levels. Although a high concentrate diet can   induce liver dysfunction, the mechanism remains unclear. Thus, we investigated whole-transcriptome changes in the liver via highthroughput RNA sequencing.
The quality control and alignment efficiency results indicated that a successful library was constructed with high sequencing quality. In this study, 22,103 mRNA, 10,376 lncRNA, 1,290 miRNA, and 3,508 circRNA were obtained from the clean data. Furthermore, we screened 43 DEG, 121 lncRNA, and 3 miRNA. These results indicate that lncRNA play an important role in SARA at the transcriptional level. Analysis of ncRNA target genes also confirmed that 42 of 43 DEG were DElncRNA target genes. In contrast, only 6 DEG were DEmiRNA target genes. We performed a qPCR analysis to validate the RNA-seq results. The qPCR results were highly consistent with the RNA-seq results, suggesting that the RNA-seq data were reliable. In addition, we counted the number of exons and the length of the open reading frames (ORF) of the mRNA and lncRNA. An overwhelming majority of the mRNA had only one exon, whereas most of the lncRNA had 2 exons (Appendix Fig. 1A and B). The ORF length ranged from 100 to 400 nt and was dominated by mRNA. The ORF length of the lncRNA was shorter and most were <150 nt (Appendix Fig. 1C and D). In general, the lncRNA in the liver were expressed at lower levels, with fewer exons and shorter ORF compared to the mRNA. This result was similar to previous studies on bovine rumen (Zheng et al., 2018).
We conducted further research on the DEG and DElncRNA. The KEGG analysis revealed that most of the DEG were annotated in the PI3K-AKT signaling pathway. AKT is an important regulatory factor in cell signal transduction pathways. It is involved in biological signal transduction, such as nutrient metabolism, cell apoptosis, and inflammation. Previous studies have reported that the PI3K- Table 1 Detailed information of the 20 DEG with the lowest FDR.

CCDC80
Fibronectin binding, protein binding, positive regulation of cell-substrate adhesion, cellular metabolic process 0.00295 À1.2423 LOC514257 Transporter activity, transmembrane transport 0.00295 À1.09758 PLEKHF1 Apoptotic process, positive regulation of autophagy, cellular response to stimulus 0.00295 1.144043 HAMP Hormone activity, cellular iron ion homeostasis, defense response to bacterium 0.003941 À1.14331 DEG ¼ differentially expressed genes; FDR ¼ false discovery rate; GO ¼ Gene Ontology; FC ¼ fold change. Q. Chen, C. Wu, Z. Yao et al. Animal Nutrition 7 (2021 AKT pathway is essential for the synthesis of milk components, such as lipids and lactose (Oliver and Watson, 2013). Furthermore, a high concentrate diet can increase lipopolysaccharide (LPS) from the gastrointestinal tract into the peripheral circulation (Emmanuel et al., 2008). Studies on dairy cow mammary epithelial cells have confirmed that LPS significantly disturbs cell proliferation and lipid metabolism via the PI3K-AKT signaling pathway (Liu et al., 2015).
Little information is available about the roles of lncRNA in the liver during SARA. Hence, we annotated the DElncRNA target genes to the GO and KEGG databases. The GO enrichment analysis showed that oxidationereduction processes were highly enriched and contained the second most targeted mRNA. Similar to the DEG results, the KEGG classification for the DElncRNA target genes indicated that most were classified to the PI3K-AKT pathway. The KEGG enrichment analysis also showed that the DElncRNA target genes were highly enriched in the proteasome, peroxisome, and the HIF-1 signaling pathways. These pathways are strongly correlated with inflammation and energy metabolism; therefore, they could be potentially involved in the regulation of SARA. For example, Xu et al. (2015) reported that LPS derived from the rumen upregulates mRNA expression of peroxisome proliferator-activated receptor-a, carnitine palmitoyl transferase 1, and acyl-CoA oxidase 1 in the liver, possibly due to the increased energy demand to resist inflammation during SARA. Kim et al. (2007) showed that activated HIF-1 induced by LPS binds directly to the tumor necrosis factor-a (TNF-a) promoter to increase TNF-a expression in hepatocytes. However, the specific roles of such pathways during SARA require further study. Subsequently, we selected 20 candidate DEG and DElncRNA with regulatory potential for further study. Table 1 shows that these DEG were mainly involved in oxidoreductase activity, stress, metabolism, the immune response, cell adhesion, apoptosis, and Fold change indicates the ratio of average expression of SARA-affected lncRNA relative to controls. Fold difference for each gene was calculated via the FPKM values of CON and SARA. GAPDH was used as the reference gene for genes expression. The data were analyzed by an independent-samples t-test using the Compare Means of SPASS 11.0 for Windows (StaSoft Inc, Tulsa, OK, USA). Values are mean ± SEM. *P < 0.05, **P < 0.01 (CON vs SARA). qPCR ¼ quantitative PCR; FPKM ¼ fragments per kilo-base of exon per million fragments; CON ¼ control; SARA ¼ subacute ruminal acidosis; GAPDH ¼ glyceraldehyde-3-phosphate dehydrogenase. proliferation. This result was highly consistent with our previous studies. Our previous studies on dairy goats reported that a high concentrate diet increases plasma acute phase protein and cortisol concentrations, and induces activation of the hypothalamicepituitaryeadrenal axis. This result suggests that goats with SARA experience stress (Dong et al., 2013;Jia et al., 2014). Furthermore, many studies have found that a high concentrate diet causes liver metabolic abnormalities, including glucose metabolism, lipid metabolism, and protein metabolism (Dong et al., 2013;Xie et al., 2015). In addition, a high concentrate diet induces oxidative stress in the liver. The enzyme activities of SOD, glutathione peroxidase, and catalase change during SARA. TAOC significantly decreases, whereas malondialdehyde notably increases (Abaker et al., 2017). Our comparative proteomic results suggested that a high concentrate diet activates the inflammatory response and decreases antioxidant capacity (Duanmu et al., 2016). Some of these DEG are worthy of note, as they were strongly correlated with the disease process. For example, AKR1B10, which had the greatest fold-change, not only participates in cell proliferation but also participates in lipid metabolism (Cao et al., 2014). AKR1B10 promotes cell growth and proliferation but it regulates lipid synthesis via stabilizing acetyl-coenzyme A carboxylase a (Ma et al., 2008). However, its role in SARA has not been reported.
Most of the target genes of the 20 candidate DElncRNA involved in SARA overlapped with the 20 DEG in Table 1, suggesting that such DEG may be regulated by lncRNA. It would be important to explore and validate the functions of the lncRNA to alleviate the negative effects of SARA. For example, DDH3, which is a DElncRNA target gene, has a strong ability to resist reactive oxygen species production in vitro, and inhibiting its expression can lead to cell death (Chen et al., 2008). Therefore, increased DDH3 expression in the liver via lncRNA may be an effective strategy to reduce oxidative stress that has been induced by a high concentrate diet. Although lncRNA research in ruminants is limited compared to humans and rodents, the findings of the present study suggest it is of value for further studies on SARA.
The result of whole transcriptome indicated that SARA induced metabolic disorders and oxidative stress. Therefore, we detected plasma triglyceride and glucose concentrations, as well as hepatic antioxidant and lipid peroxide profiles. Although plasma glucose did not show an obvious difference, there was a significant decrease in plasma triglyceride concentration. In accordance with the result of whole transcriptome, decreased TAOC, CAT and increased MDA suggested that the SARA group suffered from hepatic oxidative stress. In addition, there was no significant difference in the dry matter intake, daily yields of milk, milk lactose and protein. However, milk fat dramatically reduced in the SARA group (data are not shown). Such results were highly consistent with our transcriptome data. For example, AKR1B10, DDH3 and ACSM3, which are related to fatty acid synthesis or antioxidation significantly decreased in the SARA group.
In this study, transcriptome analysis was used to reveal the DEG in the liver of SARA dairy cows compared to normal counterparts. The DElncRNA target genes were mainly related to some cellular process or signaling pathways by KEGG pathway analysis. However, some limitations of this study should be mentioned. Considering the animal welfare and the 3Rs principles, the number of dairy cows used in this study was limited. Moreover, the underlying mechanism of altered lncRNA expression in the liver of SARA animals is still unclear. The target molecules as well as the regulatory functions of these DElncRNA still need further investigation.

Conclusion
This study identified 10,376 lncRNA, 1,290 miRNA, 3,508 circRNA, and 22,103 mRNA in Holstein cow livers. Among these RNA, 121 lncRNA, 3 miRNA, and 43 mRNA were differentially expressed. Our results indicate that mRNA and lncRNA rather than miRNA or circRNA play major roles in SARA at the transcriptional level. Moreover, a bioinformatics analysis showed that the mRNA were strongly correlated with the lncRNA. The GO and KEGG analysis revealed that the DEG were mainly related to oxidoreductase activity, stress, metabolism, the immune response, cell apoptosis and proliferation. Similarly, DEG and DElncRNA target genes were enriched in pathways, such as the PI3K-AKT, proteasome, and HIF-1 signaling pathways. Such pathways are related to inflammation, metabolism, and cell apoptosis and proliferation. Taken together, our results suggest that DEG, mainly regulated by lncRNA at the transcriptional level, and lncRNA play an important regulatory role in SARA. Above all, the findings from the current study provide valuable resources for further studies on SARA.

Conflict of interest
We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, and there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the content of this paper.