High-Throughput Sequencing Reveals Transcriptome Signature of Early Liver Development in Goat Kids

As a vital metabolic and immune organ in animals, the liver plays an important role in protein synthesis, detoxification, metabolism, and immune defense. The primary research purpose of this study was to reveal the effect of breast-feeding, weaning transition, and weaning on the gene expression profile in the goat kid liver and to elucidate the transcriptome-level signatures associated with liver metabolic adaptation. Therefore, transcriptome sequencing was performed on liver tissues, which was collected at 1 day (D1), 2 weeks (W2), 4 weeks (W4), 8 weeks (W8), and 12 weeks (W12) after birth in Laiwu black goats at five different time-points, with five goats at each time point. From 25 libraries, a total of 37497 mRNAs were found to be expressed in goat kid livers, and 1271 genes were differentially expressed between at least two of the five time points. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses revealed that these genes were annotated as an extracellular region fraction, exhibiting monooxygenase activity, positive regulation of T cell activation, mitotic spindle mid-region assembly, cytokinesis, cytoskeleton-dependent cytokinesis, regulation of cytokinesis, regulation of lymphocyte proliferation, and so on. In addition, it mainly deals with metabolism, endocrine, cell proliferation and apoptosis, and immune processes. Finally, a gene regulatory network was constructed, and a total of 14 key genes were screened, which were mainly enriched for cell growth and development, endocrine, immune, and signal transduction-related pathways. Our results provide new information on the molecular mechanisms and pathways involved in liver development, metabolism, and immunity of goats.


Introduction
The liver is a central organ that plays a key role in nutrient metabolism, including the metabolism of proteins, fats, and carbohydrates [1][2][3]. Furthermore, liver also plays a crucial role in immunity, mainly due to the fact that the blood contains large amounts of gastrointestinal antigens circulating in the hepatic vasculature, helping to form an immune barrier against pathogens [4].
During embryonic and neonatal periods, the liver functions primarily as a key hematopoietic organ [5]. During postnatal development, this organ transitions to mature metabolic function, a process that dramatically changes liver gene expression and promotes a shift from a hepatic microenvironment to one that is extensively involved in metabolism and detoxification [6][7][8][9]. At present, most studies on the molecular mechanisms related to the early development and functional transformation of the liver mainly focus on the fetal The experimental animals used in this study were Laiwu black goat, which is a local goat breed with both cashmere and meat produced in the central mountainous area of Shandong Province. It has germplasm characteristics, such as high reproduction rate, excellent meat quality, and strong disease resistance [27]. All goat kids were obtained from Shandong Fengxiang Livestock Breeding Technology Co. Goat kids lived with their dams and suckled from their mother as the only food source for the first 30 days, which included the colostrum stage (D1), the middle nursing milk stage (W2), and the late nursing t milk stage (W4). Goats were supplemented with solid concentrate starting at 30 days, and fermentation of the solid concentrate was the only food source for the goat kid after 60 days. Based on this feeding strategy, we defined age after solid concentrate supplementation as the mixed-feeding stage (W8) and the starter-feeding stage (W12). All goat kids had ad libitum access to water and concentrate. On day 1, five newborn goat kids were sacrificed 6 h after suckling colostrum, and all other groups of goat kids were sacrificed at their respective cut-off dates to collect liver tissues (five goat kids in each time point). Overall, a total of 25 healthy goats (third parity) were selected as experimental animals. The collected liver samples were snap-frozen in liquid nitrogen.

Library Preparation and Sequencing
RNA-seq transcriptome strand libraries were set up from 5 ug of total RNA using the Illumina TruSeq TM Total Stranded RNA Library Preparation Kit (San Diego, CA, USA). The paired-end transcriptome libraries were sequenced by Shanghai Majorbio Bio-Pharm Biotechnology Co., Ltd. (Shanghai, China) using the Illumina NovaSeq6000 system.

Read Mapping and Transcriptome Assembly
SeqPrep (https://github.com/jstjohn/SeqPrep, accessed on 18 January 2021) and Sickle (https://github.com/najoshi/sickle, accessed on 18 January 2021) were used as default parameters for cropping and the quality control of raw-paired-end reads. The clean reads were then aligned with the goat reference genome (Assembly: GCA_001704415.1) in the direction mode using the HISAT2 software [28]. The associated reads for each example were assembled by the StringTie [29] in a reference-based approach (Goat reference genome annotation: Capra hircus ARS1). The expression level of each transcript was calculated with the fragment-per-kilobase-per-million-reads (FPKM) method.

Differential Expression Analysis and Functional Enrichment
Essentially, DESeq2 [30] was performed for the differential expression analysis, and differential expression genes (DEGs) with fold change >2 or <−2 and a Q value ≤ 0.05 were considered to be significantly differentially expressed genes. In addition, the mean FPKM values for at least one group were not <1.
The K-means (k = 4) clustering of the relationship between samples and transcripts was performed based on the expression levels of the transcripts, and the clustering results were presented in heatmaps. The expression scales of the 25 groups of samples were normalized by the Z-Score, and the corresponding heatmaps were constructed using MEV 4.9.0 [31].
In addition, we performed a functional enrichment analysis, including GO and KEGG, to determine which DEGs were significantly enriched in the GO terms and the pathways when compared to the transcriptome-wide background at a Bonferroni-corrected p-value ≤ 0.05. Goatools [32] and KOBAS [33] were used to perform the GO function enrichment and KEGG pathway analysis.

Construction of the Protein-Protein Interaction (PPI) Network
Across the 25 libraries, in accordance with the KEGG results, four clusters of DEGs were significantly enriched in 34 pathways containing a total of 148 genes. All 148 genes were used for the PPI network construction. Their interactions were predicted after in the database of STRING [34]. The protein-coding genes regulatory network was constructed using the Cytoscape [35].

Quantitative Real-Time PCR (qRT-PCR) Validation
A total of 15 DEGs were selected randomly for RT-qPCR analysis. The primers for controls were designed based on the sequences of the selected DEGs. The primers used in this study are shown in Table 1. cDNA was synthesized following the instructions of the Evo M-MLV RT Kit with gDNA Clean for qPCR II (AG11711, ACCURATE BIOTECHNOLOGY (Hunan) Co., Ltd., Changsha, China). RT-qPCR was performed using the SYBR Green Premix Pro Taq HS qPCR Kit (AG11701, ACCURATE BIOTECHNOLOGY (Hunan) Co., Ltd., Changsha, China) and the Roche LightCycler ® 96 System (Roche, Pleasanton, CA, USA). There were reaction volumes of 20 µL, 0.2 µL of forwarding primer, 0.2 µL of reverse primer, 10 µL of 2 × SYBR Green Pro Taq HS Premix, 2 µL of cDNA, and 7.6 µL of the RNase free water. The reaction conditions of qPCR are divided into two steps, the first step is 95 • C for 30 s, and the second step is 95 • C for 5 s, then 60 • C for 30 s, and this cycle was repeated 40 times. All reactions were performed in triplicate. β-actin (ACTB) was used as an internal reference gene. After amplification, the gene expression was calculated using the 2 −∆∆CT ± SEM [36]. Correction for multiple tests were determined with one-way ANOVA, least-significant difference (LSD), and Tukey's multiple range tests. Table 1. Primer sequences for qRT-PCR testing of randomly selected DEGs.

Gene Name
Up Down

Identification and Analysis of DEGs
As listed in Table S4 and Figure 1, 1271 differential genes were used in the analyses. Cluster plots of differential expression patterns between the groups revealed that 1271 DEGs were categorized into four clusters: 424 DEGs in the first cluster, 248 DEGs in the second cluster, 217 DEGs in the third cluster, and 382 DEGs in the fourth cluster. The first cluster showed several DEGs between the D1 and other groups, and their expression levels decreased with age in Laiwu black goat ( Figures 1A and 2A). The second cluster showed a high expression of DEGs in the W2 and W8 groups and a low expression in the other groups ( Figures 1B and 2B). The third cluster showed a high expression of DEGs in the D1 and W4 groups and an upregulation of the low expression in the other groups ( Figures 1C and 2C). The fourth cluster showed several DEGs between the W12 and other groups, and their expression levels increased with age in the Laiwu black goat ( Figures 1D and 2D). The fourth clusters of DEGs were analyzed by GO and KEGG enrichment analysis.

GO and KEGG Enrichment Analysis
GO and KEGG enrichment analyses were performed on DEGs to further confirm their involvement in the regulation of biological processes. The results of the GO analysis are depicted in Table S6. During the enrichment analysis of the first cluster of DEGs, 711 terms were obtained, all of which were not significantly enriched. In the second cluster of DEGs, 804 terms were obtained, all of which were not significantly enriched. In the enrichment analysis of the third cluster of DEGs, 690 terms were obtained, of which 34 terms were significantly enriched, including 19, 0, and 15 terms for the cellular component (CC), molecular function (MF), and biological process (BP) categories, respectively (FDR < 0.05). Microtubule (GO:0005874) and mitotic spindle midzone assembly (GO:0051256) revealed the highest enrichment significance in the CC and BP terms, respectively. In the enrichment analysis of the fourth cluster of DEGs, 1006 GO terms were obtained, of which 25 GO terms were significantly enriched, with 1, 1, and 23 GO terms in the CC, molecular function (MF), and BP terms, respectively (FDR < 0.05). The extracellular region     Notably, in the first cluster of DEGs in the KEGG pathways, endocrine correlation was recorded. In the second cluster of DEGs in the KEGG pathways, three pathways were associated with amino acid metabolism and another one with endocrine relevance. In the third cluster of DEGs in the KEGG pathways, three pathways were associated with cell growth and death, one with the endocrine system, and one with cancer. Among the fourth cluster of DEGs in the KEGG pathways, 14 were associated with immunity and disease, four with amino acid and lipid metabolism, five with signal transduction, and one with the digestive system. These results indicated differences in the transcript levels in the liver of goat kids. DEGs of the first cluster, as well as the second cluster, were enriched to the same endocrine-related pathway. However, the Q values were different, and the upregulated genes were not identical in both comparisons.

Screening of the Key Regulatory Genes
To further identify the key regulatory genes that contributed to differential enrichment in the early developmental, metabolic, and immune-related pathways in the goat kid liver, PPI network analysis was performed (Figure 3). Based on the KEGG pathway enrichment outcomes, these four classes of DEG were significantly enriched in 34 pathways containing 148 genes.
All 148 genes were used for the PPI network construction, and 690 interactions between 148 genes were obtained. We sorted all genes in the network according to the degree of interaction and filtered out the 14 most important genes as potential key regulatory genes: CD4, MYC, CCND1, ESR1, ZAP70, CD5, CDK1, IL7R, PLK1, AURKA, CCNA2, CD274, PPARG, and SELL. According to the KEGG enrichment analysis, 38 pathways were significantly enriched, of which three pathways were associated with cell growth and death, four with the endocrine system, 24 with immunity and cancer, and seven with signal transduction. In addition, 15  third cluster of DEGs in the KEGG pathways, three pathways were associated with growth and death, one with the endocrine system, and one with cancer. Among the fou cluster of DEGs in the KEGG pathways, 14 were associated with immunity and dise four with amino acid and lipid metabolism, five with signal transduction, and one w the digestive system. These results indicated differences in the transcript levels in the li of goat kids. DEGs of the first cluster, as well as the second cluster, were enriched to same endocrine-related pathway. However, the Q values were different, and the upre lated genes were not identical in both comparisons.

Screening of the Key Regulatory Genes
To further identify the key regulatory genes that contributed to differential enr ment in the early developmental, metabolic, and immune-related pathways in the g kid liver, PPI network analysis was performed (Figure 3). Based on the KEGG pathw enrichment outcomes, these four classes of DEG were significantly enriched in 34 pa ways containing 148 genes. All 148 genes were used for the PPI network construction, and 690 interactions tween 148 genes were obtained. We sorted all genes in the network according to the gree of interaction and filtered out the 14 most important genes as potential key regulat

Validation of DEGs by qRT-PCR
To validate the accuracy of the sequencing and analysis outcomes, we selected differently expressed mRNAs to validate the qRT-PCR results (Figure 4). qRT-PCR validation results were found to be roughly consistent with the RNA-seq results.

Validation of DEGs by qRT-PCR
To validate the accuracy of the sequencing and analysis outcomes, we selected d ferently expressed mRNAs to validate the qRT-PCR results (Figure 4). qRT-PCR valid tion results were found to be roughly consistent with the RNA-seq results.

Discussion
Liver development involves significant changes in gene expression, which is mediated by transcriptional and post-transcriptional control [37]. This organ is an important hub for several physiological processes, including macronutrient metabolism, regulation of immune system support, blood volume, endocrine control of growth signaling pathways, lipid and cholesterol homeostasis, and the breakdown of exogenous compounds, such as drugs [38].
During early postnatal development in male rats, decreased β2-adrenergic receptormediated glucose mobilization may lead to lowered transcription of the receptor-encoding genes in the liver [39]. It is reported that liver growth factor promotes neurogenesis and neuronal survival, migration of new neurons, and striatal growth. The transcription factors HNF3, MAZ, and Sp1 are required for elevated levels of S-protein gene expression in hepatocytes [40]. Research shows that liver X receptor α induces monounsaturated fatty acid synthesis in goat epithelial cells via the Srebp-1-dependent control of stearoyl-CoA desaturase 1 [41]. As is reported, antioxidant protection may play a role in chicken embryonic development, and that functional changes in postnatal liver development, correlating with chicken liver maturation [42].
Presently, various studies have focused on the metabolism and cellular differentiation of the liver. Studies on the transcriptome of the animal liver based on different developmental periods, as well as dietary structure, have been reported [43][44][45][46]. To understand the signaling pathways and key regulatory genes of the goat kids' liver and functions related to immune regulation, this study used high-throughput transcriptome sequencing to detect the gene expression levels in the liver bodies of the goat kids at five different times points. More than 93% of Q30 was achieved in 25 samples. Compared with the reference genome, all genes identified in the library could be assigned to the reference genome of the goat, with the lowest assignment rate of 96.94%. In previous mammary gland and goatskin studies, the mapping did not exceed 80%, which indicates the high sequencing quality of our study [47].
To investigate the function of DE mRNAs with different expression profiles, we analyzed the DE mRNA expression patterns and identified four clusters. The DE mRNA expression profiles were related to the developmental stages in goats. We found that four clusters of target genes were enriched in the KEGG signaling pathway. DE mRNA was expressed at low levels at birth and at high levels in the liver of weaned goats (Cluster 1). In addition, they were significantly enriched in the PPAR signaling pathway, which is a widely studied nuclear receptor that regulates and controls metabolic changes in humans and animals. The PPAR is a nuclear receptor that mediates nutrient-dependent transcriptional activation and regulates the metabolic networks through energy homeostasis. Through different mechanisms, nutrients can regulate PPAR, thereby ultimately contributing to the prevention of various metabolic disorders [48]. Peroxisome proliferator-activated receptors include PPAR-α, PPAR-β/δ, and PPAR-γ, which play important roles in glucose and lipid metabolism [49]. Studies have shown that at knockout of PPAR-γ, there is an interaction between glucose and lipid metabolism in adipose tissue, muscle, and liver, which indicates that PPAR-γ plays an important role in maintaining systemic glucose and lipid homeostasis [50]. Mammalian milk is usually rich in lipids, and fatty acid metabolism is an important energy source for neonates in the early postpartum period [13,18]. The liver regulates lipid metabolism by associating with adipose tissue, while related transcription factors may contribute to fatty acid synthesis through the transfer of fatty acid transport protein (FATP) and fatty acid-binding protein (FABP) [38]. The liver is able to use fatty acids as an internal source of energy to other organs via an oxidative pathway, essential for organisms that experience extreme fasting or consume very low levels of dietary carbohydrates [51]. In our study, FATP, FABP, PPARG, and fatty acid transport-related genes upregulated, and PPARG is the key gene. We accordingly hypothesized that FATP and FABP act on fatty acid transport during lipid metabolism by activating PPARG through transport-related transcription factors acting on adipocytes, contribute to the prevention of various metabolic disorders, and provide energy. We recorded high expression of DE mRNA in the W2 and W8 groups for DEGs and low expression in the other groups (Cluster 2). Notably, the second group of DEGs was also significantly enriched in the PPAR signaling pathway. In our study, genes enriched downstream of the PPAR signaling pathway were also involved in cholesterol metabolism. In addition to lipid homeostasis, the liver is critical for the homeostasis of cholesterol in the body. Excess cholesterol from the diet and newly synthesized cholesterol is known to lead to inappropriate cell membrane dynamics [51]. We speculate that the PPAR signaling pathway may be activated by cholesterol in the diet during the mid-lactation and mixed milk feeding phase or by self-synthesized cholesterol. DE mRNA was highly expressed in the D1 and W4 groups and poorly expressed in the other groups (Cluster 3). The third group of DEGs was significantly enriched in five KEGG pathways (FDR < 0.05), namely cell cycle, progesterone-mediated oocyte maturation, cellular senescence, viral carcinogenesis, and the p53 signaling pathway. Most of these pathways were associated with cell growth and development. CCNB1, CDK1, and RRM2 have been reported to be enriched in the p53 signaling pathway and are believed to act as potential biomarkers and therapeutic targets for hepatitis B virus-associated hepatocellular carcinoma [52]. CDK1 was among the significantly enriched genes in the p53 signaling pathway in our study, and CDK1 is also a key regulatory gene. At the heart of the cell cycle, CDKs propel cells through the various phases of the cell cycle. In the G1 phase, cyclin D-CDK4/6 (early) and cyclin E-CDK2 (late) prepare to enter the S-phase, at which point cyclin A-cdk2 takes over and coordinates replication, followed by activation of cyclin A/B-CDK1 to facilitate the transition from the G2 to the M phase [53]. CDK1 can replace other CDKs, and it has been demonstrated to be sufficient to power the mammalian cell cycle [54]. CDK1 functions outside mitosis as a general translational activator. It allows protein synthesis to directly adapt to the rate of cell proliferation [55]. We hypothesized that the goat kids liver cells grow rapidly at 2 and 8 weeks of age. There were several DEGs between the W12 and other groups, and their expression levels increased with age for the Laiwu black goats (Cluster 4). Among all of the significantly enriched pathways of Cluster 4 genes, most of them are related to immunity, and the top five most-enriched pathways were hematopoietic cell lineage, primary immunodeficiency, Th1 and Th2 cell differentiation, toxoplasmosis, and intestinal immune network. Among them, CD4 gene was identified as the top one potential key regulatory gene of PPI network. CD4+ T cells are essential for achieving an effective immune response to pathogens, and they are activated after interacting with antigen MHC complex and differentiate into Th1 and Th2 subtypes, mainly according to the cytokine environment in the microenvironment [52]. The balance between Th1 and Th2 cells is important for normal immune response, and Th2 cell differentiation is usually promoted by Th1 response to maintain immune balance [56]. In our study, CD4 genes were up-regulated with the increase of goat kid ages. We found that CD4 genes were significantly enriched in the Th1 and Th2 cell differentiation signaling pathways and that the PPI regulatory network was a key regulatory gene. As an important immune organ, the liver is involved in regulating T cell differentiation and proliferation, and the proportion of CD4+ T cells in the liver of adult animals is significantly higher than that in the peripheral blood [57,58]. The gradually increased expression of the CD4 gene found in this study indicates that the amount of CD4+ T cells in goat liver and their differentiation are gradually activated and enhanced over time during the neonatal stage of goat kids. We accordingly hypothesized that CD4 genes play an important role in the normal immune responses by acting in the Th1 and Th2 cell differentiation signaling pathways.
A total of 690 interactions among 148 genes were obtained in the PPI regulatory network constructed from the KEGG pathway enrichment results. A total of 14 key regulatory genes were screened in order of the degree of interaction: CD4, MYC, CCND1, ESR1, ZAP70, CD5, CDK1, IL7R, PLK1, AURKA, CCNA2, CD274, PPARG, and SELL. MYC has been found to be unregulated in almost half of the human solid tumors and leukemias, and increased levels of MYC has been attributed to retroviral transfection or transduction of circulating mammalian cells that reduces the need for growth factors, prevents cell cycle secretion, accelerates cell distribution, and increases the number of cells [59]. The oncoprotein MYC is an important transcription factor that regulates the growth, proliferation, and expression of genes involved in different cellular signaling pathways [60]. MYC dysregulation has been causally linked to the development, persistence, and progression of cancer, although possible mechanisms by which oncoproteins play a role in these processes include increased cell proliferation, inhibition of cell death, metabolic regulation, promotion of angiogenesis, and regulation of stem cell formation [61][62][63]. We thus speculated that MYC genes play an important role in early development-related mechanisms in the liver of goat kids. Thus, we believe that the functions of these key genes deserve further investigation.

Conclusions
In this study, the transcriptome of the goat livers at five developmental stages during the juvenile age were sequenced, and DEGs functional enrichment analysis showed a major concentration in the pathways related to early development, substance metabolism, and immunity. This study provides important information for facilitating the study of the role of the liver in the early development, metabolism, and immunity at a young age. Our results extend the current goat mRNA database and contribute to the better understanding of the mechanisms of early development from the perspective of ruminant liver development.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13050833/s1, Table S1: Reads of the 25 Libraries; Table S2: Summary of sequencing read alignment to the goat reference genome; Table S3: Statistics of known and novel transcripts in each group; Table S4: The expression of 1271 DEGs in each group; Table S5: GO annotation for differentially expressed genes between groups; Table S6: KEGG results for differentially expressed genes between groups.

Informed Consent Statement: Not applicable.
Data Availability Statement: The original contributions presented in the study are publicly available. All data can be downloaded from the Gene Expression Omnibus (GEO) database, and the accession number is: GSE194190.