The experiments were conducted at the National Origin Farm of Xiangxi Cattle, which is located 10 km south of Huayuan County, 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 cadmium-treated 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.

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.

Different concentrations (0, 0.1, 0.2, 0.4, 0.6, 1.0, 2.0, 3.0, and 4.0 μg/mL) of HNO₃ solution were prepared by the National Center for Analysis and Testing of Non-ferrous Metals and Electronic Materials.

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₃:HCLO = 7:3) was added and heated to near evaporation. The residue was dissolved in 1.0% HNO₃. 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.

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 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. Paired-end 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 5-adenylated, 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.

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.

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 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 Y-axis 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 generated using the OmicStudio tools (https://www.omicstudio.cn/tool). rMATS (v3.2.5) (42) was used to detect splicing events in the samples.

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).

cDNA was synthesized by reverse transcription using the Superscript III reverse transcription kit (ABI Invitrogen, United States). qRT-PCR was performed for genes associated with Cd exposure using a fluorescent quantitative PCR instrument (Applied Biosystems, United States). Specific quantitative primers (KuMei, Changchun, China) for the 18 transcripts are shown in 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₂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.

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.

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).

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.

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 Cd-treated 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).

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).

Real-time fluorescence quantitative PCR was used to validate the 11 DEmRNAs, four DEmiRNAs, and four DElncRNAs selected for this study. Among the DEmRNAs, stearoyl-CoA desaturase (SCD), voltage-gated channel auxiliary subunit gamma 4 (CACNG4), CYP7A1, insulin-like growth factor 1 (IGF1), and C–X–C motif chemokine ligand 3 (CXCL3) were significantly upregulated. The difference in CXCL3 expression was extremely significant. In addition, LIPG, ATPase H⁺/K⁺ transporting the non-gastric alpha2 subunit (ATP12A), apolipoprotein A (ApoA4), bone morphogenetic protein 8A (BMP8A), and bone morphogenetic protein 8B (BMP8B) were significantly downregulated DEmRNAs. Among the DEmiRNAs, bta-miR-12051, bta-miR-211, bta-miR-222, and bta-miR-11986c were significantly downregulated. Among the DElncRNAs, 107132706, 104972497, 790183, and 104973640 were significantly upregulated. It is worth mentioning that the difference in the expression levels of 104972497 and 790183 was extremely significant (Figure 8). The qRT-PCR results for these DEGs were consistent with the sequencing data.