Milk production by the sow is a major limiting factor in the growth and survival of her litter. Knowledge of porcine milk composition, as well as understanding genetic factors underlying its variation, is a matter of ongoing interest. In this study, we performed the first exhaustive characterization of the porcine milk transcriptome derived from whole milk samples. The goal was to highlight differences in samples collected during early and mid-lactation and compare transcriptome profiles of dams across multiple parities.
A limited number of milk transcriptome studies in livestock have been conducted [9, 21-25], and to date, RNA-Seq of porcine milk has not been reported. Non-invasive sampling of the transcriptome of milk-producing cells via RNA secreted into porcine milk provides a powerful window into the biology of swine lactation. Lemay et al.  showed that RNA extracted from milk provides a better representation of RNA from mammary epithelial cells than RNA from whole mammary tissue.
Total RNA was isolated from 130 fresh whole milk samples (65 colostrum and 65 mature milk) from dams across four parities. In most milk transcriptome studies, milk is fractionated, and RNA is extracted from somatic cells, milk fat, or whey. Total RNA concentrations tend to be higher in the milk fat and somatic cells than in the whey fraction, while RNA integrity of somatic cells is higher than those of milk fat and whey [27, 28]. Low RIN values in this study (average RIN = 4.0) are likely due to the presence of small amounts of cytoplasmic material in milk fat globules , bacteria and small RNAs (miRNA) in the fat fraction , and degraded and/or free RNAs. Each milk fraction has its own place in research settings. The advantages and disadvantages of each RNA source has previously been summarized . In this study, we chose to utilize whole milk samples in order to capture the broader transcriptomic signatures of porcine colostrum and milk. We were able to process the samples much more quickly than had we fractionated the milk, and our sample represents the entirety of what is being ingested by the growing piglet.
Libraries were sequenced to an average depth of 46 million reads per library. A depth of 40 million reads is considered sufficient for reliable detection of major splice isoforms for abundant and moderately abundant transcripts . When generating our sequence data, we targeted a depth of 50 million reads per library. However, there was considerable variation in sequence depth across libraries. Some of this variation can be attributed to technical aspects of NGS technology, such as the stochasticity of sequencing, RNA quality, and library preparation.
A total of 44,234 transcripts were identified in this study, of which approximately 95% are annotated in the current swine genome build. Transcripts corresponded to 17,740 unique gene loci, including 16,515 known porcine genes. The number of expressed genes is comparable to those reported in similar studies in sheep  and goat . A smaller number of expressed genes (~13,500) was reported in the buffalo milk transcriptome . This discrepancy is likely to due to the swine, sheep, and goat reference genomes being more complete and of higher quality.
In general, gene expression values covered a wide range of magnitudes (Figure 7), and the gene expression profiles of colostrum and mature milk were not highly correlated (Pearson correlation coefficient 0.222). There was a large overlap (12 out of 15) in the top fifteen most abundantly expressed genes in colostrum and mature milk (Table 3). Among the top expressed genes were CSN3, CSN2, CSN1S1, LALBA, FASN, EEF1A1, TPT1, SAA3, and WAP, which have been previously identified among the top expressed genes in milk samples from other species [23-25, 32, 33].
As expected, many of the top expressed genes were related to biosynthesis of milk proteins. Expression levels of CSN2, CSN3, CSN1S1,LALBA, and WAP, which encode for the synthesis of the main milk proteins casein and whey, increased from early to mid-lactation stages. A similar gene expression pattern has been identified in a previous swine study , as well as in goat , cattle , and sheep . High expression of the EEF1A1 gene is also related to high levels of milk protein synthesis, as EEF1A1 is one of the most abundant protein synthesis factors . Additionally, ribosomal proteins RPLP0 and RPS2 were among the top expressed genes in colostrum and exhibited a slight decrease in expression during mid-lactation. These genes were also found to be highly expressed in early lactation in buffalo .
In addition to milk protein synthesis genes, genes associated with milk fat were among the top expressed genes, and their expression was nearly constant across lactation stages. Milk fat composition is known to influence piglet growth and development . The FABP3 gene, which is involved in the uptake and transport of fatty acids, has been linked to milk fat synthesis in cattle . FASN is directly involved in most of the short and medium-chain fatty acids in milk , and XDH is involved in the formation of the lipid droplet in milk .
The 20 most significant DET for milk type and interaction are given in Tables 4 and 5, respectively. Many of the most significant DET were associated with genes involved in milk fat synthesis. Transcript MSTRG.335722.71455, which is associated with insulin induced gene 1 (INSIG1), was found to be up-regulated in mature milk samples. INSIG1 is known to regulate the expression of sterol regulatory element-binding protein 1 (SREBP1; also denoted SREBF1), which is central to milk fat synthesis [20, 39]. Although not one of the most significant DET, transcript MSTRG.239932.53831, associated with SREBP1, was identified as a DET for the milk type x parity interaction. Transcripts MSTRG.189750.42732 (THRSP gene), MSTRG.108970.23990 (SP1 gene), MSTRG.283058.62377 (ANXA7 gene) are other milk fat synthesis genes among the most significant DET. In general, expression of milk fat synthesis genes was up-regulated in mature milk samples compared to colostrum, which agreed with expression patterns observed across bovine lactation stages . Our results highlight that the transition from swine colostrum to mature milk is marked by a shift from high protein contents to high fat and lactose contents .
In this study, DEG were identified by aggregating P-values across transcripts associated with each gene via the Lancaster method, rather than using gene read counts directly. Using this approach not only maintains both transcript and gene-level resolution, but also bypasses issues of different variances and directions of change across constituent transcripts. This method outperforms other gene-level methods and provides a coherent analysis between transcripts and genes .
One of the major aims of this study was to evaluate DEG between lactation stages across multiple parities. Progeny born to multiparous sows generally exhibit superior growth performance compared to those born to primiparous sows. However, colostrum and milk composition profiles, are highly similar across parities .
Glucose transport is a major precursor to lactose synthesis, which is synthesized in the Golgi vesicle of mammary secretory alveolar epithelial cells during lactation . Glucose-6-phosphate transporters SLC37A2 and SLC37A3, glucose cotransporter SLC5A8, and glucose transporter SLC2A5 were identified as DEG for the milk type main effect. Glucose transport across the plasma membrane of mammalian cells is carried out by two distinct processes one of which involves glucose transporters from the GLUT gene family (encoded by SLC2A genes) and the other which involves glucose transporters from the SGLT family (encoded by SLC5A genes). Both the SLC2A5 and SLC5A8 genes were up-regulated in colostrum. Crisá et al.  identified significant up-regulation of members of the SLC2A gene family and polysaccharide and glycosamino-glycan binding molecular function to be enriched in goat colostrum samples compared to mature milk.
Members of the SLC35 gene family encode nucleotide sugar transporters localizing at the Golgi apparatus and/or the endoplasmic reticulum. These transporters transport nucleotide sugars pooled in the cytosol into the lumen of these organelles, where most glycoconjugate synthesis occurs . Currently, the SLC35 gene family is comprised of 31 genes which are divided into 7 subfamilies, SLC35A to SLC35G . GDP-fucose transporters SLC35C1 and SLC35C2 were identified as DEG for the milk type main effect, with SLC35C2 up-regulated in colostrum and SLC35C1 down-regulated. SLC35A5 and SLC35G1 were identified as DEG for the milk type x parity interaction, with both genes being up-regulated in mature milk in the first three parities and up-regulated in colostrum during the fourth parity (Figure 8). Crisà et al.  identified 3 DEG from the SLC35 family that were up-regulated in goat colostrum compared to mature milk, as well as enrichment of glycosaminoglycan binding molecular in colostrum. Consistent with this result, we also identified the enrichment of glycosaminoglycan binding molecular function, with 29 of the 124 porcine genes associated with the GO term being present in our milk type DEG set. These results support the hypothesis that oligosaccharide metabolism decreases over the course of lactation as oligosaccharide concentrations in milk have been shown to peak in colostrum and decrease over the remainder of lactation in human , goat , and pig .
Forty-two genes from the PI3K-Akt pathway were found to be differentially expressed. Although not statistically significant after FDR-correction, this pathway had a nominal P-value of 0.025 in the enrichment analysis. The PI3K-Akt pathway is a key signaling node for lactogenic expansion and differentiation of the luminal mammary epithelium, as numerous signaling pathways that regulate lactogenic development converge on PI3K-Akt, including the insulin-like growth factor 1 receptor (IGF1R), RANKL and RANK, integrins, and PRLR-to-JAK2-to-STAT5A pathways . In general, expression of DEG in the PI3K-Akt pathway was up-regulated in colostrum. The janus kinase 1 (JAK1) gene, a key component of the PI3K-Akt pathway, had an FDR-corrected P-value of 0.02 but was filtered out of the DEG list based on the log2 fold change threshold. JAK1 was more highly expressed in mature milk, and the difference in expression level between colostrum and mature milk was more pronounced in the later parities (Figure 9). Using a novel mammary gland-specific JAK1 knockout mouse model, Sakamoto et al.  demonstrated that JAK1 is essential to involution, the mammary gland remodeling process at weaning where the milk-producing epithelial cells are replaced with adipocytes . The potential effects of involution on mammary development and milk yield during the next lactation are not well understood. Ford et al.  noted that mammary glands in sows that were suckled during lactation were larger than non-suckled glands at the end of the involution process, suggesting a possible beneficial effect on redevelopment during the next gestation.
Several GO terms related to immune response, particularly leukocyte activation, were significantly enriched in the DEG for the milk type main effect (Table S10). The majority of genes in these pathways were up-regulated in colostrum. This finding is consistent with the widely accepted theory that immunoglobulins and immune cells are transferred to the neonate via colostrum. In pigs, the epitheliochorial nature of the placenta prohibits transfer of maternal immune cells and immunoglobulins to the fetus, and thus, the piglet relies on the successful absorption of colostral components to acquire maternal immunity . Proinflammatory cytokines play an important role in the development of the neonatal immune system by mediating the early local and systemic responses to microbial challenges . A total of 37 DEG were associated with cytokine secretion, including interleukin 17 receptor C (IL17RC), interleukin 27 receptor A (IL27RA), tumor necrosis receptor superfamily members 1b and 15 (TNFRSF1B and TNFRSF15), and tumor necrosis factor superfamily member 4 (TNSFF4). Several other genes in these gene families were shown to be up-regulated in early porcine lactation by Palombo et al. .
Antimicrobial proteins naturally present in colostrum and milk can kill and inhibit a broad spectrum of bacteria . Milk is also known to exert chemotactic activity on neutrophils , an important innate host defense against microorganisms. The chemokine superfamily encodes secreted proteins involved in immunoregulatory and inflammatory processes. The CXC chemokine ligand 16 (CXCL16), which encodes a chemokine antimicrobial protein , was up-regulated in colostrum samples (Figure 10). Palombo et al.  identified significant up-regulation of CXCL2 and CXCL10 in day 1 postpartum swine milk samples compared to colostrum. Many of the main chemokines (CXCL2, CXCL8, CXCL9, CXCL10, CXCL11, CXCL13, CXC14, and CXCL16) were expressed in our samples, and expression patterns for CXCL10 and CXCL2 were consistent with the findings of Palombo et al. . Our results differed from previously published results in that CXCL8 was the most abundantly expressed chemokine in both colostrum and mature milk, and CXCL3 was not expressed in our samples. One factor contributing to this discrepancy was the use of the improved reference genome (Sscrofa 11.1), where many of the gaps and misassemblies present in the Sscrofa 10.2 genome build were resolved and the annotation was significantly improved. These results suggest that chemokine ligands may play an important role in the transition from colostrum to mature milk in swine, likely helping prompt recruitment of neutrophils.