Comprehensive Transcriptome Profiling of Dairy Goat Mammary Gland Identifies Genes and Networks Crucial for Lactation and Fatty Acid Metabolism

Milk fatty acids secreted by the mammary gland are one of the most important determinants of the nutritional value of goat milk. Unlike cow milk, limited data are available on the transcriptome-wide changes across stages of lactation in dairy goats. In this study, goat mammary gland tissue collected at peak lactation, cessation of milking, and involution were analyzed with digital gene expression (DGE) sequencing to generate longitudinal transcript profiles. A total of 51,299 unigenes were identified and further annotated to 12,763 genes, of which 9,131 were differentially expressed across various stages of lactation. Most abundant genes and differentially expressed genes (DEGs) were functionally classified through clusters of euKaryotic Orthologous Groups (KOG), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. A total of 16 possible expression patterns were uncovered, and 13 genes were deemed novel candidates for regulation of lactation in the goat: POLG, SPTA1, KLC, GIT2, COPS3, PDP, CD31, USP16/29/37, TLL1, NCAPH, ABI2, DNAJC4, and MAPK8IP3. In addition, PLA2, CPT1, PLD, GGA, SRPRB, and AP4S1 are proposed as novel and promising candidates regulating mammary fatty acid metabolism. “Butirosin and neomycin biosynthesis” and “Glyoxylate and dicarboxylate metabolism” were the most impacted pathways, and revealed novel metabolic alterations in lipid metabolism as lactation progressed. Overall, the present study provides new insights into the synthesis and metabolism of fatty acids and lipid species in the mammary gland along with more detailed information on molecular regulation of lactogenesis. The major findings will benefit efforts to further improve milk quality in dairy goats.


INTRODUCTION
Goat milk contains many macro-and micro-nutrients that are essential for the growth and development of a newborn, and has been recognized as beneficial for human health (German et al., 2002;Casado et al., 2009). The nutritional value of goat milk is mainly attributable to its fat and protein fractions, which are critically important to provide both energy and essential nutrients for the humans (Shi et al., 2015). One of the major goals of goat farming is to improve milk quality through the alteration of milk composition, according to specific needs of target groups such as infants or immune-compromised individuals (Wickramasinghe et al., 2012). To achieve this goal, a thorough understanding of milk components and their regulatory factors is required.
The mammary gland has been widely-used as a research model to investigate the mechanisms of milk formation and secretion (Shennan and Peaker, 2000). However, it was not until the mid-2000s that the first study utilizing largescale transcriptome profiling was published (Rudolph et al., 2003). That original study, although broad in its focus, and a subsequent one underscored that biosynthesis and secretion of lipid in mammary gland tissue is regulated by complex gene networks (Rudolph et al., 2007;Andres and Djonov, 2010). Subsequent work in dairy cows focusing on a curated list of 45 genes associated with lipid synthesis (triacylglycerol and phospholipids) and secretion reported the first longitudinal profiles from late pre-partum to the end of subsequent lactation (Bionaz and Loor, 2008). These observations allowed for the development of a network of genes participating in the coordination of milk fat synthesis and secretion.
Despite the substantial body of work in bovine, research on the mammary gland of small ruminant species, particularly in goats, is limited. For instance, most published studies on goats have focused on functional verification of candidate genes associated with lipid metabolism using data from human, mouse and cattle, including SCD1 (Yao et al., 2017), SREBP1 (Xu et al., 2016), FASN (Zhu et al., 2014), INSIG (Li et al., 2019), and ELOVL7 . Although those efforts have enhanced our knowledgebase, economicallyimportant traits in livestock such as milk lipid synthesis and composition are genetically affected by some causal genes together with many genes with minor effects, which together result in complex regulatory networks (Goddard and Hayes, 2009;Georges et al., 2019). Extensive and detailed molecular research in dairy cows has underscored that the ruminant mammary gland is a very dynamic organ with remarkable plasticity (Akers, 2017), especially during the transition from pregnancy to lactation. However, there is a paucity of data regarding changes occurring in the goat mammary gland during the physiological transition from established lactation to dry-off to involution. Therefore, for a better understanding of regulatory mechanisms underlying milk fatty acid metabolism and lactation in dairy goats, a systematic screen of gene expression profiles is needed.
Development of high-throughput sequencing has made it feasible to generate a whole genome sequence in the domestic goat (Mardis, 2008;Metzker, 2010). In addition, mining the differentially expressed genes (DEGs) in mammary tissue across discrete stages of lactation could vastly expand our knowledge of the molecular regulation of mammary gland function in small ruminants. Therefore, in the present study we used Illumina digital gene expression (DGE) profiling sequencing to characterize transcript profiles in mammary gland tissue across 3 stages of the lactation cycle in Xinong Saanen dairy goats. Beyond expanding our mechanistic knowledge, identifying novel promising genes responsible for the genetic variation could contribute greatly to our understanding of lactation and milk fatty acid metabolism. It could help in the future to develop a molecular breeding program with the aim of improving fatty acid composition in dairy goat milk.

Ethics Statement
The present experiment was conducted in accordance with approved guidelines by the Animal Care and Use Committee of the Northwest A&F University. The committee approved all procedures and experiments with live animals.

Goat Mammary Gland Collection and RNA Extraction
Mammary gland tissue was collected based on methods described previously (Shi et al., 2015). Briefly, a total of nine healthy Xinong Saanen dairy goats (approximately 3 to 4 years old from 2 to 3 parities) were selected from the experimental farm at Northwest A&F University (Shaanxi, China). All goats were managed similarly and were fed a mixed diet comprised of corn, soybean meal, bran, rapeseed meal and mineral-vitamin premix. Approximately 2 g mammary gland tissue was collected from the mid-region of the right mammary gland of each goat at peak lactation (3 goats, 100 days postpartum; L group), dryoff (cessation of milking, 3 goats, 310 days postpartum; D group) and non-lactating/non-pregnant period (involution, 3 goats; NP group). All tissue samples were obtained under sterile conditions on the same day, harvested within 20 min for each and immediately frozen in liquid nitrogen until RNA extraction.
Total RNA extraction, mRNA purification, and cDNA library construction were conducted by LC Sciences (Houston, TX, United States). In brief, total RNA was obtained from mammary gland tissue using a total RNA purification kit (LC Sciences, Houston, TX, United States), treated with RNAase-free DNAase, and re-purified with the RNA easy kit (Qiagen, Valencia, CA, United States) following the manufacturer's instructions. Total RNA quantity and quality were analyzed with the RNA 6000 Nano LabChip Kit in the Bioanalyzer 2100 system (Agilent Technologies, CA, United States).

cDNA Library Construction and Illumina Sequencing
The mRNA was purified from total RNA using poly-T-oligoattached magnetic beads (Invitrogen, United States). Following purification, the mRNA was fragmented into small pieces using divalent cations at a high temperature. After phosphatase and polynucleotide kinase (PNK) treatment, the cleaved mRNA fragments were reverse-transcribed. The 200-300 bp fragments were purified using 6% polyacrylamide Tris-borate-EDTA to create the final cDNA library. One end sequencing of 1 × 36 bp was then performed on an Illumina HiSeq 2000 platform following the vendor's recommended protocols. The application of DGE sequencing to uncover DEGs has advantages such as low cost, rapid analysis, and high accuracy (t Hoen et al., 2008).

Data Filtering and Expression Level Determination
The raw reads were first determined by fastQC 1 and filtered by removing potential contamination, reads with unknown base more than 1 N (base not available), and low quality sequences (parameter -q 36 -p 90) using the Fastx_toolkit (0.0.13) (Lohse et al., 2012). All short reads were assembled as transcripts using Trinity software (Haas et al., 2013). As a subset of the transcript, unigene was the longest transcript of each cluster that matched the same reference gene. HTSeq v0.6.1 was used to count the reads numbers mapped to each unigene. By aligning sequences with a goat database constructed previously (BioProject ID: PRJNA243005), the gene expression level was recorded according to the sum frequencies of compared reads. Reads per kilo base of exon model per million mapped reads (RPKM) were used for expression normalization, and the gene expression RPKM values were categorized into three groups: high (≥500 RPKM), medium-to-high (10 to 500 RPKM) and low (≤10 RPKM) (Wickramasinghe et al., 2012;Canovas et al., 2014). Assembled unigenes were annotated to the reference dataset using BLAST. Databases including a manually Annotated and Reviewed Protein Sequence Database (SWISS-PROT), NCBI non-redundant protein sequence database (NR), Pfam database, KOG, GO, and KEGG were performed with an E-value cutoff of 1E-10 according to previously published studies (Ogata et al., 1999;Shi et al., 2015).

Analysis of DEG Data
The normalized ratio of the gene expression signals was log 2 transformed. The RVM (Random variance model) N 2 -test and Chi-square test were applied to filter DEGs. After analyses for statistical significance and FDR (false discovery rate), DEGs were selected according to P < 0.05 and FDR < 0.05 (Wright and Simon, 2003;Yang et al., 2005;Clarke et al., 2008). Bowtie 0.12.8 was used for abundance analysis of gene expression. Short Time-series Expression Miner (STEM) was used for expression pattern analysis (Ernst and Bar-Joseph, 2006). Using the log normalized data method, the "maximum number of model profiles" was set at 20, and the "maximum unit change in model profiles between time points" was set at 2 (Ernst and Bar-Joseph, 2006). Gene co-expression analysis (Pujana et al., 2007;Xu et al., 2013) was performed to track interactions among DEGs. Pearson correlation analysis was performed for each pair of genes, and significantly correlated pairs were used to construct the network with a threshold of 0.92 using a Perl script (Prieto et al., 2008). The number of correlated genes was used for node gene selection. Based on all DEGs identified when adding up the three pairs of comparisons (L vs. D, L vs. NP, D vs. NP), potential candidate genes regulating lactation and milk fatty acid metabolism were identified.

Dynamic Impact Approach (DIA) Analysis
Data were filtered with the threshold at fold-change >2 and P < 0.001, and were mined through an integrative systems biology approach applying the DIA method . The bovine KEGG pathways and GO biological process category database were used for functional analysis with the DIA. Detailed methodology for data analysis using DIA was described previously (Bionaz et al., 2012b). The estimate of the perturbation in a biological pathway is represented by the "Impact" while the overall direction of the perturbation is represented by the "Flux" (or direction of the impact) . The impact was obtained by combining the proportion of DEGs with the log 2 mean fold-change and meanlog P-value of genes associated with the biological term. The direction of the impact was calculated as the difference of the impact of up-regulated DEGs and downregulated DEGs associated with the biological term (Bionaz et al., 2012b). The analytical flow chart used in the present study is described in Supplementary Figure S1.

Illumina Sequencing and Gene Annotation
A total of 63,256,236 clean reads were obtained with an average of 7,028,470 (range from 3,423,245 to 10,365,063 reads) for each sample. After sequence assembly according to the goat transcriptome sequence published in our previous study (Shi et al., 2015), out of 98,864 transcripts (Supplementary File S2A) the Illumina sequencing uncovered a total of 51,299 unigenes. Among these unigenes, 48,428, 50,442, and 50,684 unigenes were associated with the lactation groups "L, " "D, " and "NP" respectively. The number of unigenes uniquely expressed in the "L, " "D, " and "NP" groups were 94, 225, and 126, respectively ( Figure 1A). Importantly, a total of 12,763 genes were identified through annotation, with 12,239, 12,713, and 12,710 genes in the "L, " "D, " and "NP" group, respectively. Among those identified genes, 3, 24, and 6 genes, respectively, were uniquely expressed in each group ( Figure 1B).

Lipid Biosynthesis and Metabolism-Associated Genes
Gene expression level was evaluated using the RPKM method. Goats in the "L" group expressed a total of 12,239 genes including 65 with the RPKM greater than 500, 1,102 genes with RPKM between 10 and 500, but most genes (11,072) had expression below 10 RPKM. According to the same RPKM classification, out of a total of 12,713 and 10,710 genes in the "D" and "NP" group there were 104 and 77 genes with RPKM greater than 500, 3,470 and 2,262 genes with RPKM between 10 and 500, while 9,139 and 10,371 genes had expression below 10 RPKM.
The most abundant genes were further filtered by focusing on fatty acid biosynthesis and metabolism as well as lipid transport and metabolism pathways or terms using the KOG, GO, and KEGG categories. Of these, the 15 most abundant genes related to lipid transport and metabolism are reported in Supplementary  Table S1. FASN, FABP4, ACBP, EBP, RNPEPL1, DEGS1, and SAPOSIN were detected in all three stages of lactation. The most abundant genes involved in intracellular trafficking, secretion, and vesicular transport were SSR4, lepB, SEC61G, AP2M1, ARF1, CHMP2A, CLTA, RAB11B, RAB1A, SDHD, and SEC61 (Supplementary Table S2).
Within the GO fatty acid metabolism category, EHHADH was the most abundant gene in lactation stages "L, " "D, " and "NP, " with RPKM of 24.134, 34.567, and 34.016, respectively. Wellknown genes involved in fatty acid metabolism such as ACAA2, ACSL and ACAA had a medium-to-high expression level across all lactation stages (Supplementary Table S3). Within fatty acid biosynthetic process in the GO category, FASN was the most abundant gene across lactation stages with 376.222, 120.146, and 537.391 RPKM in "L, " "D, " and "NP, " respectively. In addition, the most abundant genes observed in all three lactation stages were SCD, ELOVL1, DEGS1, NADH, NDUFAB1, ABDH2, MSMO1, PRKAG1, SC5D, STK11, and TECR (Supplementary Table S4).

Identification of Differential Gene Expression Patterns
A total of 29,880 differentially expressed unigenes were detected among different groups, annotated as 9,131 DEGs (Supplementary File S2B). In the "D" and "NP" groups, a total of 46 and 18 unique unigenes were annotated to 6 and 2 DEGs, respectively. Although five uniquely expressed unigenes were detected in the "L" group, they were not annotated to any functional gene (Figure 2). A total of eight unique DEGs across lactation stages are reported in Table 1.
In total, 50 GO terms including biological process, cellular component and molecular function were significantly enriched with 8,311 DEGs (Figure 3). Of these, transcription, DNAdependent and regulation of transcription had the mostabundant DEGs across 25 enriched biological processes. With 15 terms in the cellular component category, most DEGs were related to the nucleus, integral to membrane, and cytoplasm. A total of 10 terms involved in molecular function including zinc ion binding, ATP binding, and protein binding had the most DEGs (Figure 3). In addition, a total of 246 KEGG pathways were enriched 5,048 DEGs (Supplementary File S2C) with Insulin signaling, Jak-STAT signaling, Glycerophospholipid metabolism, mTOR signaling and PPAR signaling among the most-relevant to lactation and milk component synthesis regulation.
The most abundant DEGs involved in fatty acid biosynthesis and metabolism processes, lipid metabolism, protein metabolism, and lactose metabolism were used to mine within GO and KEGG categories for pathways or terms across lactation stages (Supplementary Tables S7-S10). A total of 48 DEGs were identified in terms of fatty acid metabolic processes/pathways, whereas 18 DEGs including ARL5A, MARK, STAM, and EPHA1 were simultaneously revealed by GO and KEGG (Supplementary Table S7). Among a total of 36 identified DEGs associated with fatty acid biosynthetic processes/pathways, SEMA6, RUNX1T1, and LEF1 were detected by both GO and KEGG (Supplementary Table S7). Within GO and KEGG, the well-known genes GPAT3, ADH and ALDH involved in lipid metabolic processes/pathways were among a total of 29 most abundant DEGs (Supplementary Table S8). Highly expressed DEGs such as NR1F1, EIF3S5, MYH, and ABCF2 were in the lipid biosynthetic processes/pathways (Supplementary Table S8).
FIGURE 3 | GO classifications of DEGs in goat mammary gland. Distribution of the GO categories assigned to the goat mammary gland transcriptome. DEGs were classified into three categories: biological processes, cellular components, and molecular functions. L, peak lactation period; D, dry-off period; NP, non-lactating/non-pregnant period.
trafficking, secretion, and vesicular transport with GGA, SRPRB, AP4S1 having the greatest number of correlated genes ( Table 4).
The correlated genes for node gene selection are reported in Supplementary File S2E.

DIA Analysis of DEGs
The DIA provides a summary of the KEGG pathways in the form of categories and sub-categories ( Figure 5A) that are altered across a given comparison. Details of each pathway are presented in Supplementary File S3. As shown in Figure 5A, KEGG pathway categories were more impacted in stage of lactation "L" versus "NP." Among these pathways, the category "Metabolism" was the most impacted one, followed by the category "Genetic Information Processing." With the exception of the subcategory of pathways within "Metabolism of Terpenoids and Polyketides, " all other subcategories within Metabolism had an impact value >50 in the comparison of "L" versus "NP, " with highest impact value (880) in "Biosynthesis of Other Secondary Metabolites."  Except for inhibition of "Glycan Biosynthesis and Metabolism" and "Nucleotide Metabolism, " compared with "NP, " most of the metabolic pathways were markedly activated in the "L" group including "Carbohydrate Metabolism, " "Energy Metabolism, " "Lipid Metabolism, " "Amino Acid Metabolism, " "Metabolism of other Amino Acid, " "Metabolism of Cofactors and Vitamins, " and "Xenobiotics Biodegradation and Metabolism." According to the impact value, the categories "Environment Information Processing, " "Cellular Process, " "Organismal System" were also altered in the comparison of "L" versus "NP." However, most of their flux values were slightly activated or did not change in the comparisons of "D" versus "NP" or "L" versus "NP." The top 10 overall most-impacted terms are shown in Figure 5B. The most impacted pathway was "Butirosin and neomycin biosynthesis" followed by "Glyoxylate and dicarboxylate metabolism" with flux >200 (Figure 5B). The categories "Galactose metabolism" were highly activated. In contrast, the pathways "Glycosaminoglycan biosynthesis-keratan sulfate, " "Glycosphingolipid biosynthesis -globo series, " and "Glycosphingolipid biosynthesis-ganglio series" were inhibited.

DISCUSSION
The development of next-generation sequencing has made it feasible to explore patterns of DEGs in goat mammary glands across stages of lactation. Similar to bovine, this information will help us better understand fatty acid metabolisms and tissue function at the transcriptional level (Bionaz et al., 2012a). One limitation of the current experiment was that mammary gland tissue harvested at the three stages of lactation was not from the same goats, which might overshadow the specific effects associated with lactation. A similar experimental approach for tissue collection and analyses has been utilized in previous work with dairy cows (Li et al., 2016;Cai et al., 2018). However, unlike dairy cows with a large udder mass, repeated mammary biopsies in small ruminants is riskier and prone to inducing cessation of lactogenesis. The fact that peak lactation in the goat occurs at ∼100 days postpartum meant that biopsies collected at dry-off would have taken place roughly 200 days afterward, which was deemed problematic due to the uncertainty about recovery from a biopsy. From a technical standpoint, sampling different animals at the target stages of lactation within the same day increased accuracy of sequencing partly because of the characteristics of temporality and spatiality of mRNA expression.
In the present study, we identified DEGs across three important stages of the lactation period, i.e., established lactation, dry-off, and the non-lactating/non-pregnant period. Data mining revealed key factors involved in the control of mammary gland function, as well as improved our understanding about the role of transcriptional regulation in coordinating lactation. The assembly of 98,864 transcripts was greater than that in macaque (42,702 transcripts) (Gao et al., 2013), bovine (35,945 transcripts) (Shen et al., 2015), and pig (21,331 transcripts) (Zhu et al., 2014). The fact that the latest reference genome was not used in the present study may be a potential reason for the greater number of transcripts detected. However, an advantage of the approach we used is that we could perform a better mapping of the sequences with a more robust unigene annotation. The sequences in this study were assembled according to the goat transcriptome sequences published in our previous study (Shi et al., 2015), i.e., the first study to characterize the complete transcriptome of goat mammary gland, constituting a comprehensive genomic resource available for studies of ruminant lactation. An N50 length is commonly used for assembly evaluation, and a higher number suggests higher-quality assembly (Lander et al., 2001). The N50 length of assembly in our previously published transcriptome sequences (Shi et al., 2015) was higher than other ruminant transcriptome data with a higher overall number of unique transcripts (Wickramasinghe et al., 2011(Wickramasinghe et al., , 2012, suggesting it could provide a comprehensive reference dataset of gene expression profiling for future goat mammary gland research. It is well known that FASN is a crucial gene for de novo fatty acid synthesis in mammary cells (Zhu et al., 2014). Consistent with our previous qRT-PCR analysis, the expression of FASN in mammary tissue of goats was higher during lactation than dryoff (Zhu et al., 2014). This may partly explain the enrichment of short-and medium-chain fatty acids in goat milk (Morris et al., 2007;Silanikove et al., 2010). However, the up-regulation of FASN in the non-lactating/non-pregnant period, despite the shutting down of lactogenesis within mammary cells, was somewhat unexpected. It is possible that this response was related to an additional need of fatty acids for cell proliferation (Mellenberger et al., 1973) associated with mammary gland development during pregnancy (Faucon et al., 2009), which would be important for the turnover of mammary gland (Norgaard et al., 2008).
Although determining differences between lactating and nonlactating periods has helped elucidate regulatory mechanisms that are controlling lactation (Bionaz and Loor, 2008;Li et al., 2012;Wickramasinghe et al., 2012), it is also important to uncover the DEGs that are specifically induced during pregnancy (Faucon et al., 2009). In the present study, in order to avoid the effects of pregnancy (e.g., endocrine status) on gene expression, we examined expression trends of DEGs among the three periods and detected 20 DEGs with lower and stable expression during the non-lactating and non-pregnant period compared with lactation (Figure 4, pattern 4). Although some milk fat-related genes might be important during pregnancy (Linzell, 1959), most of the genes related to milk fat formation had expression levels that were between those observed in the dry-off (pregnant) and the non-lactating (non-pregnant) periods. More functional research is needed for identifying the role of these genes.
A total of two and six DEGs were determined as uniquely expressed in the dry-off and non-lactating periods. In the nonlactating group, bone morphogenetic protein 3/3B (BMP-3B) is a cytokine that belongs to the transforming growth factor β (TGFβ) superfamily, and plays an important role in the antral nervous system as well as in bone formation and remodeling (Takao et al., 1996). The G protein-coupled receptor 64 (GPR64) is well-known as a member of G protein-coupled receptor superfamily, and is crucial for neuronal development regulation (Richter et al., 2013). Despite putative functions in development, their expression in the non-lactating period was relatively low (RPKM were 0.415 and 0.386, respectively). These results are in close agreement with the extensive remodeling of the mammary gland during the non-lactating period (McDaniel et al., 2006;Safayi et al., 2010).
Among genes with higher expression in the dry-off group, protein Myb-Like, SWIRM and MPN domains 1 (MYSM1) are important for the control of B cell development (Jiang et al., 2011). Myb proto-oncogene protein (MYB) promotes proliferation and inhibits differentiation (Hugo et al., 2013) and also has a key role in regulating stem and progenitor cells in the bone marrow, colonic crypts and a neurogenic region of the adult brain (Ramsay and Gonda, 2008). The gap junction protein beta 2 (GJB2) is important for the development of the auditory system (Bozdogan et al., 2015). All these genes had a relatively higher expression level (the RPKM was 0.958, 1.720, and 2.511, respectively) and most of them are beneficial for embryo development, probably playing a role in mammary cells during dry-off to prepare well in preparation for the impending parturition.
Co-expression analysis is a suitable method for regulatory gene selection (Xu et al., 2013). A node gene correlated with the greatest numbers of genes was considered crucial for the regulation of lactation. Our analysis revealed that POLG, SPTA1, KLC, GIT2, COPS3, PDP, CD31, USP 16/29/37, TLL1, NCAPH, ABI2, DNAJC4, and MAPK8IP3 correlated with most genes. In turn, these genes may play roles in regulating mammary cell function in the lactation stages that were examined in the current study.
Replication of DNA and repair in mitochondria is a crucial function of POLG. SPTA1 is the major protein component of the erythrocyte membrane skeleton (Huebner et al., 1985). KLC is a component of Kinesin heterotetramer and acts in binding motor protein to cargo (Hanks et al., 1988). GIT2 plays an important role in cell attachment, spreading and motility (Pahlich et al., 2006). COPS3 is essential for maintenance of cell proliferation in the mouse embryonic epiblast (Buhr et al., 2008). PDP is important for dephosphorylating pyruvate dehydrogenase in the mammalian pyruvate dehydrogenase complex which catalyzes the conversion of pyruvate to acetyl-CoA, a substrate for de novo fatty acid synthesis (Gonsalvez et al., 2007). CD31 is associated with the vascular compartment and is uniquely-positioned to mediate multiple and important cell-cell interactions involving platelets, leukocytes and endothelial cells (Zheng et al., 2005). USP is essential for maintaining the structure and function of neuromuscular junction (Verbiest et al., 2008). TLL1 is necessary for normal septation and positioning of the heart (Thomassen et al., 2009). NCAPH plays an important in regulating cell cycle activity and proliferation (Hardwick, 2008). ABI2 is involved in abscisic acid signal transduction in Arabidopsis (Dean et al., 2001). DNAJC4 is a sterol-regulatory element binding protein (SREBP)-regulated chaperone involved in the cholesterol biosynthesis pathway. MAPK8IP3 is involved in c-Jun N-terminal kinase (JNK) signaling pathway.
All the above genes have a wide number of functions associated with cellular regulation, consistent with the idea that the process of lactation is regulated and balanced by a complex network of genes (Bionaz and Loor, 2008). Consistent with previous studies, the node genes are central for lipid transport and metabolism along with intracellular trafficking, secretion, and vesicular transport. Specifically, PLA2 (Zhu et al., 2011;Shen et al., 2012;Zhang et al., 2013), CPT1 (Bovia et al., 1995), GGA (Gustavsson et al., 2014), and SRPRB (Hill et al., 1988) had the most number of correlated genes. Although the exact role of these genes in fatty acid metabolism during the various stages of lactation remains to be determined, they may represent candidate genes predicted to be beneficial for our greater understanding of milk fat synthesis.
By combining the proportion of DEGs with the log 2 mean fold change and mean -log P-value of genes associated with the biological terms, the DIA is a useful approach to estimate the biological impacts of experimental conditions and the direction of impact . It has been successfully used for studying goat mammary lipid metabolism in isolated epithelial cells . In the present study, we determined "Glyoxylate and dicarboxylate metabolism" as one of the mostimpacted metabolic pathways during lactation, and was closely associated with the biosynthesis of acetyl-CoA, the substrate for de novo fatty acid synthesis. Thus, the present results underscored alterations of milk lipid synthesis during lactation in goat mammary gland as one of the most-important functional aspects occurring in this organ.

CONCLUSION
The comprehensive transcriptome profiling of mammary gland across different stages of lactation in Xinong Saanen dairy goat revealed a total of 51,299 unigenes that were annotated to 12,763 genes. The most highly expressed genes were involved in milk fatty acid synthesis and metabolism. Eight out of 9,131 identified DEGs were uniquely expressed in the non-lactating or dry-off periods. We obtained 16 possible expression patterns among DEGs and detected 13 node genes potentially regulating functional adaptations across stages of lactation: POLG, SPTA1, KLC, GIT2, COPS3, PDP, CD31, USP16/29/37, TLL1, NCAPH, ABI2, DNAJC4, and MAPK8IP3. "Butirosin and neomycin biosynthesis" and "Glyoxylate and dicarboxylate metabolism" were the most impacted pathways, underscoring their importance for functional adaptations of the mammary gland. In addition, PLA2, CPT1, PLD, GGA, SRPRB, and AP4S1 were uncovered as node genes and participate in aspects of fatty acid metabolism. Our findings provide promising candidate genes for elucidating molecular mechanisms controlling mammary fatty acid metabolism during lactation in the dairy goat and should be explored further.

DATA AVAILABILITY STATEMENT
The sequencing data have been submitted to the NCBI SRA, and are accessible through the accession number PRJNA637690.

ETHICS STATEMENT
The animal study was reviewed and approved by Animal Care and Use Committee of Northwest A&F University.

AUTHOR CONTRIBUTIONS
JL was responsible for the experimental design. JZ, WZ, HX, HW, and HS collected the tissue samples and isolated RNA for sequencing. JZ and HS performed the analysis of sequencing data and bioinformatics analysis. CL, JZ, HS, JL, and JJL wrote the manuscript. CL, JL, and JJL supervised the entire experiment, participated in result interpretation, and manuscript revision. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We thank Prof. Cal (Northwest A&F University) for surgical sampling of mammary gland tissue, the goat farm crew for taking care of the animals, and researchers of the Dairy Goat Lab of Northwest A&F University for their useful discussion.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020. 00878/full#supplementary-material FIGURE S1 | Analytical flowchart for digital gene expression sequencing in mammary gland of dairy goat. Nine mammary gland samples were collected in the stage of lactation, dry-off and non-lactation with three biological replicates per stage.
TABLE S1 | The 15 most abundant genes involved in lipid transport and metabolism category (KOG database) in mammary gland.