Whole Transcriptome Profiling of the Effects of Cadmium on the Liver of the Xiangxi Yellow Heifer

Cadmium (Cd) is a major heavy metal toxicant found in industrial zones. Humans and animals are exposed to it through their diet, which results in various physiological problems. In the current study, the toxic effects of Cd on the liver were investigated by whole-transcriptome sequencing (RNA-seq) of the livers of Xiangxi heifers fed a diet with excess Cd. We randomly divided six healthy heifers into two groups. The first group received a control diet, whereas the second group received Cd-exceeding diets for 100 days. After 100 days, the livers were collected. A total of 551 differentially expressed mRNAs, 24 differentially expressed miRNAs, and 169 differentially expressed lncRNAs were identified (p < 0.05, |log2FC| >1). Differentially expressed genes (DEGs) were analyzed by gene ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses. We found that under Cd exposure, DEGs were enriched in the adenosine 5'-monophosphate–activated protein kinase pathway, which is involved in autophagy regulation, and the peroxisome proliferator–activated receptor pathway, which is involved in lipid metabolism. In addition, the apolipoprotein A4 gene, which has anti-inflammatory and antioxidant effects, the anti-apoptotic gene ATPase H+/K+ transporting the nongastric alpha2 subunit, and the cholesterol metabolism–associated gene endothelial lipase gene were significantly downregulated. C–X–C motif chemokine ligand 3, cholesterol 7α-hydroxylase, and stearoyl-CoA desaturase, which are involved in the development of fatty liver, were significantly upregulated. These genes revealed the main effects of Cd on the liver of Xiangxi yellow heifers. The current study provides insightful information regarding the DEGs involved in autophagy regulation, apoptosis, lipid metabolism, anti-inflammation, and antioxidant enzyme activity. These may serve as useful biomarkers for predicting and treating Cd-related diseases in the future.


INTRODUCTION
Heavy metals accumulate significantly in the soils of most parts of China, with cadmium (Cd) being the most prevalent heavy metal contaminant. Contaminated areas are mainly located in the central, southern, and southwestern parts of China, and Hunan is one of the most seriously contaminated provinces (1). There are several large mines and smelting areas in the Hunan Province, which is known for its non-ferrous metal mining activities and is a typical area of dense low-temperature deposits in China (2). Heavy metal contamination in soil is considered a potential cancercausing health risk for Hunan residents. Moreover, heavy metal pollution in rivers used for agricultural irrigation in Hunan Province is considered to indirectly affect the health of humans and animals (3). Garden soil, rice soil, vegetables, and rice in the Hunan mining area are contaminated with heavy metals (4). Cd pollution is considered to be the most serious heavy metal pollution in Hunan. The Cd content in both rice straw and endosperm samples in Hunan is higher than that in the Yangtze River delta (5). The western part of Hunan is known as Xiangxi. The Xiangxi yellow heifer is the main farmed animal in the Xiangxi region, which mainly consumes local forage. This is highly significant for the development of the Xiangxi farming industry. The beneficial characteristics of Xiangxi yellow heifers include mild temperament, resistance to heat, medium size, welldeveloped and strong bones, tender meat, and rich nutrition (6). Beef is one of the main food sources for the local people. Studies have shown that metal accumulation in livestock can be transferred to humans through the food chain (7,8).
Heavy metals greatly endanger the health of animals and humans through environmental exposure and the food chain. Cd contamination has caused several vicious incidents, including the Cd rice incident in Hunan, China, and the Itai-Itai disease public health incident in Japan. In China, Hunan, Henan, Anhui, and Jilin provinces have high levels of Cd in the soil. Cd in Hunan farmland exceeds the soil environmental background value by 30-fold (9). Cd is absorbed mainly by the respiratory and digestive tracts. Moreover, Cd adversely affects the liver, kidneys, lungs, immune system, reproductive system, and cardiovascular system (10). When Cd-containing foods and water are mistakenly consumed, Cd is transported from the intestine to the liver, causing functional disorders (11). Studies have shown that Cd interferes with collagen metabolism and extracellular signal transduction in the liver matrix, as well as cell adhesion, growth, and migration (12)(13)(14). Cd also causes morphological and structural changes in the liver, such as blurred hepatic trabeculae, increased cell size, hepatocyte deformation, nuclear consolidation, lipid droplet aggregation, and mono-nuclear cell infiltration (15). Hepatic mitochondrial function has been reported to be impaired in a rat model of Cd exposure (16)(17)(18). Cd exposure inhibits mitochondrial oxidative phosphorylation, resulting in a decrease in ATP. In addition, high concentrations of Cd inhibit basal respiration (19). Moreover, Cd can directly damage the mitochondrial structures in the liver, causing mitochondrial degeneration, vacuolation, and disruption of the mitochondrial membrane structure (20,21). In the presence of severe mitochondrial damage, the process of mitochondrial autophagy is activated. Metallothionein (MT), which binds to Cd, protects organisms from Cd toxicity. Cd decreases the levels of MT and antioxidant glutathione (GSH) and inhibits the activity of antioxidant enzymes (22). Therefore, it can induce oxidative stress (23,24). Cd can also induce systemic inflammatory responses involving innate, adaptive, and mucosal immunity (25). In addition, the increase in Ca 2+ levels accelerates reactive oxygen species (ROS) formation, also resulting in oxidative stress (26,27). Cd induces endoplasmic reticulum stress and DNA mutations by inhibiting the DNA repair system and inducing protein damage (28)(29)(30). Cd damages the adenosine triphosphate system and the tricarboxylic acid cycle, which can lead to an increase in various metabolites related to lipid metabolism and triglycerides and disrupted energy metabolism (18). Cd exposure is a potential risk factor for the development of diabetes mellitus, non-alcoholic steatohepatitis (NASH), and non-alcoholic fatty liver disease (NAFLD). Cd exacerbates hepatocyte lipid storage through dysregulation of autophagy, which aggravates hepatic steatosis and leads to hyperlipidemia (31). Cd is also associated with an increased risk of prostate, breast, lung, and hepatocellular carcinoma (32)(33)(34).
RNA-seq, characterized by low background noise and a large detection range, can provide a large amount of information on gene expression levels, facilitating the study of genetic mechanisms. Currently, it is the preferred method for studying gene expression (35,36). Some studies have claimed that heifers are potential biological carriers of heavy metal contamination in soil (8). An association exists between metal levels in heifer tissues and the risk of human exposure. In addition, the liver is the most prominent target organ of Cd toxicity (37).
Therefore, in this study, the Hunan landmark product Xiangxi yellow heifers were used as experimental animals. Then, RNA-seq was used to study the effects of feeding Cd overload diets on the expression of mRNAs, miRNAs, and lncRNAs in the liver of Xiangxi yellow heifers to understand their transcriptional biology. At present, there are numerous studies on the hepatotoxicity of Cd, but there is a relative lack of information on the hepatic transcriptomics of Cdexposed Chinese Xiangxi yellow heifers. Hunan Province is rich in mineral resources. However, Cd pollution is present. Cd has great potential to endanger the health of the local dominant breed of Xiangxi yellow heifers. Therefore, this study hypothesized that certain transcription factors and metabolic pathways would be affected by Cd exposure during foraging in Chinese Xiangxi yellow heifers. The current experiment aimed to find the pathological changes that occurred in the liver of Xiangxi yellow cattle after feeding Cd-containing forage. This experiment also aimed to analyze the gene expression at the RNA level and the pathways that are significantly enriched in Xiangxi yellow cattle exposed to cadmium using transcriptomics. Our findings provide new insights into the potential health risks of Cd, improve the current understanding of the toxic effects of Cd on the liver, and help to prevent Cd toxicity.

Materials, Animals, and Treatment
The experiments were conducted at the National Origin Farm of Xiangxi Cattle, which is located 10 km south of Huayuan County, FIGURE 1 | The content of cadmium in liver tissues. Cadmium content in the liver tissues of the control group (Control) and cadmium-treated group (Cadmium), as obtained by graphite furnace atomic absorption spectrophotometry (µg/kg). Data were analyzed by independent samples t-test. The significance level was set at * p < 0.05.
Hunan, China. All experimental protocols and procedures used in this study followed Chinese standards for the use and care of research animals (SOP-017). A total of six healthy female crossbred (Angus bull × Xiangxi local cattle) heifers weaned (BW = 78.50 ± 18.45 kg) at 6 months of age were selected from a paternal halfsib family. They were randomly allocated to the control group (control) or cadmiumtreated group (cadmium), with three heifers in each group. Heifers were housed in an open-sided, naturally well-ventilated barn, and the groups were separated using railings. Thus, three heifers within each group were housed together in one pen. The feeding experiments were conducted for 109 days. The first 9 days were the pre-feeding period, during which the test heifers were dewormed and acclimatized to the diet. The following 100 days were the regular experimental feeding periods. During the experiments, the diets of both groups consisted of concentrate, paddy, and corn straw silage ( Table 1). These diets were mixed daily and provided to the cattle a total mixed ration. In China, the standard for agricultural soil cadmium is 0.3 mg/kg (38). The corn silage used for the cadmium-treated group was made from a whole-plant corn species with strong cadmium enrichment ability (Yuqing 23). This corn was planted in soil with excessive cadmium content (3-5 mg/kg) in Yonghe Town, Liuyang City, Hunan Province. The corn silage used for the control group was made from normal whole-plant corn that was planted in soil with a Cd content of 0.2 mg/kg. The Cd concentrations of the diets in the control and cadmium groups were 0.72 and 6.74 mg/kg, respectively. The orts were removed one time daily before the morning feed, and a new feed was delivered two times daily at 7:00 a.m. and 17:00 a.m. Cattle within each group shared one feed trough. All the cattle had free access to water. A part of the liver was placed in 10% formalin and used to prepare histological sections for morphological analysis. The remaining liver tissue samples were frozen at a −80 • C in a refrigerator for RNA-seq.

Preparation of Liver Tissue Sections
Liver tissues were fixed in 10% formalin for 24-48 h, followed by dehydration in different grades of ethanol. Next, 6-µm-thick sections were cut from the prepared liver blocks, affixed to glass slides, rehydrated, and stained with hematoxylin and eosin. The prepared slides were observed under a microscope equipped with a microphotographic system. First, 0.5 g of the liver was dried in a desiccator at 75 • C for 1 h. Then, 3 mL of the digestion solution (HNO 3 :HCLO = 7:3) was added and heated to near evaporation. The residue was dissolved in 1.0% HNO 3 . The samples to be measured and cadmium standard solutions with different standard concentrations were placed in the Z-5000 graphite furnace atomic absorption spectrophotometer. Finally, the determination was carried out separately according to a pre-set procedure (drying at 80-140 • C for 40 s; ashing at 300 • C for 20 s; and atomization at 1,500 • C for 5 s). Finally, a concentration-absorbance curve was automatically generated by the machine.

Library Preparation and Sequencing
Liver tissue samples frozen at −80 • C were used for RNA-seq. A 50-100 mg sample was collected from each group and added to 0.5 ml TRIzol (Invitrogen, Carlsbad, CA, United States) on ice. Then, the samples were mixed using a homogenizer (Cebo, Shanghai, China), and 0.5 mL TRIzol was added. The tissue and TRIzol were transferred together into a 1.5 mL EP tube and left at 37 • C for 5 min. The sample was centrifuged at 12,000 × g for 5 min and the sediment was discarded. Then, 0.2 mL of chloroform was added. The samples were shaken for 15 s and incubated at room temperature for 3 s. The samples were then centrifuged at 4 • C for 20 min at 12,000 × g. When the samples appeared to be stratified, the RNA was mainly in the upper aqueous phase; therefore, the aqueous phase was transferred to a new tube. Immediately after, 100% isopropanol (0.5 mL) was added. The samples were mixed by inversion and left at 37 • C for 10 min. The supernatant was removed by centrifugation at 12,000 × g for 10 min at 4 • C. Then, 1 mL of 80% ethanol was added and the sample was mixed. After centrifugation at 7,500 × g for 5 min at 4 • C, the supernatant was discarded. The samples were then left at room temperature for approximately 5-10 min. Finally, 50 µL of DEPC water was added to solubilize the RNA, and the RNA solution was obtained by centrifugation under a light flow. Total RNA concentration and quality were determined using a NanoDrop spectrophotometer and an Agilent 2,100 Bioanalyzer (Thermo Fisher Scientific, MA, United States).
Two new libraries, lncRNA/mRNA and miRNA, were constructed. In the lncRNA/mRNA library, the Ribo-Zero Magnetic Kit (Epicenter) was used to remove 5.8S rRNA, 18S rRNA, 28S rRNA, 12S rRNA, and 16S rRNA from total RNA. Single-stranded cDNA was synthesized using random primers and reverse transcriptase. DNA polymerase I and RNase H (Invitrogen) were used to synthesize the second-strand cDNA. At this point, the RNA template was removed and The total mapping ratio is the overall mapping ratio of the sequenced reads to the reference genome and the uniquely mapping ratio is the mapping ratio of the sequenced reads to the reference genome. dTTP was replaced with dUTP. The double-stranded cDNA product was then subjected to 3'-adenylation and splice ligation. The products were subjected to several rounds of PCR and thermal denaturation to a single strand and then cyclized by splint oligo to obtain single-stranded circular DNA. Pairedend sequencing (100 bp) was performed using a BGISEQ-500/MGISEQ-2000 system (BGI-Shenzhen, China). Each sample from the lncRNA/mRNA library yielded an average of 13.65 G of data. For miRNA sequencing, 18-30 nt of small RNA (Ladder Marker, Takara) was isolated and purified from total RNA by polyacrylamide gel electrophoresis (PAGE) gel cutting. The 5adenylated, 3-blocked single-stranded DNA junction was ligated to the 3' end of the purified RNA, and the 5' junction was ligated to the 5' end of the RNA. RNA with 3' and 5' junctions attached was synthesized into cDNA by reverse transcription extension using RT primers bearing UMI tags. The cDNA was subsequently amplified by PCR. The PCR products were separated by PAGE to obtain fragments in the range of 110-130 bp. Fragments were purified using a QIAquick Gel Extraction Kit (QIAGEN, CA, United States) and quality-controlled using an Agilent 2,100 Bioanalyzer (Thermo Fisher Scientific, United States). The purified library products were subsequently sequenced at the 50 bp single-end on the DNBseq platform (BGI-Shenzhen, China). The raw data obtained from sequencing were converted into raw sequence data (raw data or raw reads) by base calling and stored in the fastq file format, which contains the sequence of reads and sequencing quality information. The miRNA library yielded an average of 23.77 M data per sample. miRanda (3.3a) and TargetScan (7.1) were used to predict the target genes of the miRNAs, and the David database (6.8) was used for pathway enrichment analysis of the miRNA target genes.

Transcriptome Data Analysis
The sequencing data in the fastq format were filtered using SOAPnuke (v1.5.2) (39) for reads with unknown base N content >5%, reads with splice contamination, and reads without insert tags. The clean reads obtained were then mapped to the reference genome Bos taurus (vGCF_000003205.7_Btau_5.0.1) using Bowtie2 (v2.2.5) (40). mRNA expression was calculated using RSEM (v1. 2.12), which uses the normalization method fragments per kilobase of exon model per million mapped fragments (FPKM). miRNA expression was calculated using reads of exon model per million mapped reads values. Finally, the DESeq2 (41) method was used to determine the differential expression of two groups of genes (three biological replicates per group) based on a negative binomial distribution model with high sensitivity and accuracy. The obtained P was calculated using the likelihood-ratio test to represent the difference between the two groups of samples at the inclusion level. For DEG analysis, a P ≤ 0.05 and | log2FC |>1 were considered significant.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Enrichment Analyses
To annotate the biological functions of the genes, all differentially expressed mRNAs (DEmRNAs) were aligned against the KEGG and GO databases. The number of genes per GO term and the enriched KEGG pathways were calculated. The candidate genes were then compared to all background genes of the species using the phyper function (https://en.wikipedia.org/wiki/ Hypergeometric_distribution) based on a hypergeometric test to obtain the p. The p-values were corrected using the Benjamini-Hochberg algorithm to obtain FDR values. FDR ≤0.05 was used as the threshold value. Finally, the GO terms and KEGG pathways Frontiers in Veterinary Science | www.frontiersin.org that met this condition were recognized as significantly enriched compared with the background genes of this species.
In the histogram of GO enrichment, the X-axis represents the number of genes annotated to a particular term. The Yaxis represents the number of selected genes annotated to a particular term. In the KEGG pathway diagram, the rich factor on the X-axis represents the ratio of the number of DEGs in the metabolic pathway to the number of genes annotated in the pathway, with larger values indicating greater enrichment. The Y-axis represents the pathway to which the gene was annotated. Box plots and variable shear plots were constructed to determine the quality of the RNA-seq data. The box plot demonstrates the distribution of gene expression levels for each sample and displays the dispersion of the data distribution, which was FIGURE 4 | Heatmap of DEmRNAs, DElncRNAs, and DEmiRNAs clustering analysis. Red represents upregulated genes, and blue represents downregulated genes.

RNA Extraction
Liver tissues of Xiangxi yellow heifers were removed from the −80 • C refrigerator. Total RNA was extracted according to the instructions. A small portion of the liver was homogenized in 1 mL of TRIzol. Then, 200 µL of chloroform was added and the sample was shaken vigorously. The sample was centrifuged at 1,200 × g for 3 min and the supernatant was transferred to an absorption column. Next, 400 µL of 96% ethanol was added and recentrifuged at 12,000 × g for 1 min. The supernatant was removed from the absorption column and 500 µL of 3 M sodium acetate trihydrate solution was added. It was then subjected to centrifugation at 12,000 × g for 1 min. Finally, 70% alcohol (400 µL) was added and the mixture was centrifuged at 12,000 × g in triplicate. The supernatant was removed from the absorption column, and 46-50 µL of DEPC was added to the RNA pellet. Total RNA was collected in Eppendorf tubes and centrifuged for 1 min at 10,000 × g. Total RNA content was measured using a NanoDrop spectrophotometer and an Agilent 2,100 Bioanalyzer (Thermo Fisher Scientific, MA, United States).
Quantitative RT-PCR cDNA was synthesized by reverse transcription using the Superscript III reverse transcription kit (ABI Invitrogen, United States). qRT-PCR was performed for genes associated  Supplementary Table S1. The total volume of the qRT-PCR mixture was 20 µL, which contained 10 µL of SYBR (Takara Bio), 1 µL of the forward primer, 1 µL of the reverse primer, 2 µL of the cDNA, and 6 µL of ddH 2 O. The reaction conditions were as follows: 1 cycle at 95 • C for 5 min, 40 cycles at 95 • C for 10 s, 58 • C for 20 s, 72 • C for 20 s, and 1 cycle at 95 • C for 15 s. The magnitude of the fold differences in gene expression levels was calculated using the 2 − CT method. An independent-samples t-test was used to test the difference between the cadmium-treated and control group. Statistical significance was set at p < 0.05.

Statistical Analysis
t-Tests were performed on all experimental data in triplicate using IBM SPSS Statistics 22 and GraphPad Prism 5. All results are expressed as mean ± SEM, and the significance level was set at * p<0.05, * * p < 0.01, and * * * p< 0.001.

Cadmium Accumulation and Pathological Changes in the Liver
The Cd content in the liver of heifers in the Cd-treated group was significantly higher than that in the control group. This finding indicates that Cd accumulated in the liver of Xiangxi yellow heifers fed a high-Cd diet for 100 days (Figure 1).
Histopathological observations reveal that the liver of heifers in the cadmium-treated group had expanded sinusoids compared to the control group (Figure 2).

Analysis of RNA-seq Libraries
To evaluate the toxic effects of Cd on the liver of Xiangxi yellow heifer, the liver was collected for whole transcriptome analysis. As shown in Table 2, a total of 132-140 million raw reads were obtained. After removing reads with an unknown base N content >5% and low-quality reads, 131-141 million clean reads remained, and 13-14 Gb of clean bases were acquired.
Clean reads were mapped to the reference genome. The localization rate was ∼91.8% ( Table 3). Approximately 88.6% of the reads could be uniquely mapped to the Bos taurus genome. The box plot represents the degree of dispersion and volatility of the data distribution ( Figure 3A). Figure 3A shows that the boxes were symmetrical about the middle line, indicating that the data were normally distributed. The smaller height of the boxes was reflected by the smaller degree of fluctuation in the data. Good standardization of the data was indicated by the consistency of the data distribution of the samples in Figure 3A. As shown in Figure 3, mRNA was mainly generated by exon skipping to produce splice isoforms of different mRNAs and regulate gene expression.

DEGs Analysis
Cluster analysis of all DEGs was performed to observe gene expression patterns. Heat maps were generated, with red and blue representing upregulated and downregulated genes, respectively    ( Figure 4). The overall expression pattern of genes in the Cdtreated group was quite different from that of the control group, and the three samples within the same group showed similar expression patterns. This indicates that the samples from both groups were reproducible, demonstrating that the derived data were credible. Volcano plots of the DEGs showed the differences in the distribution of gene expression levels between control and Cd-treated samples. A total of 551 DEmRNAs were identified in the control and cadmium-treated bovine livers. Of these, 318 genes were upregulated, and 233 genes were downregulated. Compared with the control group, 24 differentially expressed miRNAs (DEmiRNAs) were identified. Of these, 13 genes were upregulated and 11 downregulated in the cadmium-treated group. In total, 169 differentially expressed lncRNAs (DElncRNAs), were identified. Of these, 73 genes were upregulated and 96 genes downregulated (Figure 5).

Functional Analyses of DEGs
A total of 551 significant DEGs in the liver, were subjected to GO and KEGG pathway analyses to determine the enrichment of functions and action pathways by DEGs. Based on the results of the functional classification of DEGs in the categories of biological process (BP), cellular component (CC), and molecular function (MF), 92 GO terms were obtained. The terms of interest were included in the first 30 enriched GO terms (including oxidation-reduction process, cholesterol biosynthetic process, cell cycle, cyclin B1-CDK1 complex, cyclin-dependent protein kinase holoenzyme complex, and oxidoreductase activity) ( Figure 6A). Table 4 lists the DEmRNAs involved in these functionally related events.
The KEGG pathways that were significantly enriched (padj <0.05) were metabolic, peroxisome proliferator-activated receptor (PPAR) signaling, AMPK signaling, and steroid biosynthesis pathways (Figure 6B; Table 5). Table 6 shows the DEmRNAs involved in these signaling pathways in Cd-exposed livers of Xiangxi yellow heifers and the positions of these DEmRNAs in the pathways (Figure 7).

DISCUSSION
Hunan is rich in mineral resources, but the problem of Cdcontaminated soil is severe. Cd is a toxic heavy metal with high mobility in soil, and can greatly endanger the health of humans and animals when it enters their bodies through diet and water. The Xiangxi yellow heifer is a special breed in Hunan that feeds on local forage. Therefore, they are highly susceptible to Cd contamination. Chronic Cd exposure is mainly deposited in the liver and can induce liver inflammation (43), and even acute and chronic liver lesions (44). Moreover, redox imbalance, mitochondrial dysfunction, fatty acid oxidation inhibition, and apoptosis are considered the main adverse effects of Cd exposure (26,45,46). In this study, the whole transcriptome of liver tissues of Xiangxi yellow heifers fed high and low Cd concentrations was determined, and DEGs were analyzed to study the mechanisms of Cd toxicity in the heifer liver. The critical DEmRNAs involved in Cd detoxification between the control and cadmium-treated groups were analyzed using KEGG and GO analyses.
The predominant toxic effect of Cd exposure is ROS production, which promotes the production of mitochondrial phagosomes that engulf damaged mitochondria (47). Cd exposure has been reported to increase the expression of adenosine 5'-monophosphate (AMP)-activated protein kinase (AMPK) (48). AMPK is also involved in autophagy regulation (49). This is consistent with the KEGG enrichment results of the current study. The effect of Cd on mitochondrial autophagy in the liver of Xiangxi yellow heifers has been demonstrated. SCD is enriched in the AMPK signaling pathway and mediates fatty acid metabolism in the liver and adipogenesis (50). SCD-1, a common isoform of this enzyme (51), is a crucial gene in hepatic steatosis and plays a straightforward role in intracellular fat accumulation in NAFLD. In addition, excessive accumulation of liver fat leads to liver inflammation (52).
Research has found that the mRNA levels of SCD-1 are also significantly higher in rats administered with CdCl 2 in their drinking water (53). NAFLD is associated with elevated mRNA levels of SCD-1. High Cd exposure significantly upregulated SCD expression in our study. This demonstrates that the entry of Cd into the bodies of humans and animals can seriously disrupt the metabolic system and even cause metabolic diseases. In the current study, metabolic pathways were enriched with the most DEmRNAs. The endothelial lipase gene (LIPG) can hydrolyze phosphatidylcholine (PC) in HDL and release free fatty acids (FFAs) (54). In the present study, LIPG expression was significantly downregulated, suggesting that HDL hydrolysis was inhibited in the liver. Thus, cholesterol in the surrounding tissues could not be transported to bile acids in time, and fat accumulated in the liver. Cholesterol 7α-hydroxylase (CYP7A1) is in the metabolic pathway and belongs to the cytochrome P450 enzyme (CYP450) family. It is correlated with inflammatory cytokines and lysosomal function (55,56). CYP7A1 overexpression strongly induces autophagy (57), and CYP7A1 is required for autophagy induction. In the present study, the expression levels of CYP7A1 mRNA were significantly altered, as shown by qRT-PCR. This indicates that Cd induces autophagy in the liver. The PPAR signaling pathway was also significantly enriched in the liver. PPARs are a class of nuclear receptors that bind to fatty acids and their derivatives (58). PPAR-α has been reported to have transcriptional regulatory effects on genes involved in gluconeogenesis and lipid transport, and PPAR-γ controls glucose metabolism, adipogenesis, and adipocyte differentiation. Part of PPAR-γ inhibits inflammatory responses by competitively inhibiting inflammatory pathways and inflammatory mediator production (59). PPAR-β promotes hepatic fatty acid oxidation and limits inflammation (60)(61)(62). Evidence suggests that PPARγ is closely related to the toxicity of heavy metals and that Cr inhibits fibroblast differentiation by reducing PPAR-γ (63). SCD and CYP7A1 genes enriched in the PPAR signaling pathway are closely related to cellular inflammation and adipogenesis in humans and animals. In this study, the mRNA expression levels of SCD and CYP7A1 were significantly upregulated in heifers under chronic Cd exposure. These findings show that this pathway is indeed involved in the toxic effects of heavy metals. In addition to the DEmRNAs in these three significantly different pathways, there were several DEGs related to liver physiological functions after Cd exposure (CACNG4, CXCL3, ATP12A, APOA4, BMP8A, BMP8B, and IGF-1). CACNG4 regulates voltage-gated calcium channels, and studies have shown that Cd exposure increases intracellular Ca 2+ concentration (64). Elevated Ca 2+ channel activity leads to increased CACNG4 expression (65,66). The elevated levels of CACNG4 expression in this study suggest an imbalance in cellular Ca 2+ concentrations in the liver owing to Cd exposure. CXCL3 is a pro-inflammatory cytokine that is induced alongside MMP12 in mice undergoing oxidative stress (67). CXCL3 exerts its anti-apoptotic effect mainly by counteracting the reversal of the cationic gradient during apoptosis (68). The significant upregulation of CXCL3 in this study demonstrates that Cd causes liver inflammation. ATP12A is a non-gastric, functionally active H + /K + -ATPase that is expressed mainly in the human and porcine respiratory tract, human prostate tissue, and bone marrow mono-nuclear cells. (69)(70)(71). In the present study, ATP12A expression was measured in the livers of heifers. The marked downregulation  of ATP12A expression levels suggests that Cd reversed the cationic gradient in the liver, leading to apoptosis. APOA4, a lipid-binding protein, is involved in the synthesis of cholesterol, steroids, and other lipids, and improves glucose homeostasis (72). The anti-inflammatory and antioxidant activities of APOA4 are crucial for maintaining hepatic homeostasis (73,74). In our study, APOA4 was downregulated, indicating that Cd injury affected fat metabolism, inflammation, and oxidative stress in the liver. Furthermore, BMP8A regulates ROS homeostasis via Nrf2 (75), and BMP8B is involved in the regulation of growth traits in heifers (76). In addition, IGF-1 modulates cell proliferation and apoptosis through the PI3/AKT and MAPK pathways (77). Overall, these seven genes were not enriched in significantly different pathways in the transcriptome data. Nevertheless, they indicated that Cd exposure is involved in hepatotoxicity regarding several crucial metabolic pathways and life processes, including inflammation, apoptosis, oxidative stress, and lipid metabolism. These findings provide insight into the mechanisms underlying Cd-induced hepatotoxicity in mammals. Studies have demonstrated that miRNAs and lncRNAs can regulate the cell cycle, apoptosis, DNA damage, and inflammation in the mechanisms underlying Cd toxicity (78)(79)(80)(81). Partial DElncRNAs and DEmiRNAs were selected for further analyses. They were verified by RT-PCR and matched with the data obtained by sequencing, which proved the accuracy and reliability of the sequencing data. The target mRNAs of the DElncRNAs and DEmiRNAs were enriched in the AMPK signaling pathway, glycerophospholipid metabolism, metabolic pathways, and peroxisomes. Therefore, the DElncRNAs, DEmiRNAs, and DEmRNAs identified in this study may be used as biomarkers of Cd toxicity in the future, which is of major research significance.

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/, PRJNA801770; https://www.ncbi.nlm.nih. gov/, PRJNA801802.

ETHICS STATEMENT
The animal study was reviewed and approved by the Tab of Animal Experimental Ethical Inspection, JLU. Written informed consent was obtained from the owners for the participation of their animals in this study.