Milk yield and milk composition of lactating cows are the most important economic traits in dairy cattle. In the past century, genetic selection on milk yield has improved milk production in cattle. With the development of quantitative trait loci (QTLs), genome-wide association studies (GWAS) and RNA sequencing (RNA-seq) technologies, a large number of candidate genes and SNPs associated with milk performance traits have been identified, such as DGAT1 and GHR genes [1–6]. Previous studies described the associations between DGAT1 K232A polymorphism and milk production traits [7, 8]. Blott et al. identified a significant SNP (BFGL-NGS-118998) inside the GHR gene that has an important role related to milk traits [9].

Mammary gland is an important organ to synthesize and secrete milk. The mammary epithelial cell has a remarkable ability to convert blood circulating nutrients into milk components [10]. Thus, almost all the studies related to milk performance traits are based on mammary gland. For example, Cui et al. collected mammary gland samples from four lactating cows after slaughter. They used RNA-seq to generate the bovine mammary transcriptome with high and low milk fat and protein percentage [11]. Finucane et al. compared bovine mammary expression profiles before and after parturition using microarray [12]. However, there are still limitations in analysis of milk performance traits using bovine mammary tissue, such as sampling difficultly, high cost of sampling and tremendous damage to the lactating cows, which resulted in small sample size.

Milk is composed of fat, protein, lactose, minerals, vitamins and water, and all these nutrients derive from blood [10, 13–15]. Some data shows that synthesizing 1L milk requires 400-500L blood flow though the mammary gland. So plenty of blood is essential for synthesizing milk. Previous studies associated with milk traits were only focused on blood physiological and biochemical indexes, and there are very few studies associating gene expression with milk traits in cattle blood. A gene expression profile from blood provides new opportunities to clarify the basic mechanisms of complex traits in cattle milk. Besides, whole blood is a complex mixture of cells and can accurately reflect the physiological condition and health level of cows. Many studies have used blood to diagnose disease in dairy cattle, such as mastitis [16, 17]. Most importantly, blood is easier to sample in comparison with other tissues and involves limited handling of animals. Furthermore, the lactation process requires multiple tissues and organs to complete, such as mammary gland, liver and adipose tissue. Hence, blood has the potential ability to represent milk performance traits more directly and comprehensively. Sandri et al. (2015) analyzed the gene expression profile in the blood related to the gene merit for milk productive traits using microarray [18]. Marcel et al. analyzed porcine peripheral blood mononuclear cells transcription profile of pigs with divergent humoral immune response and lean growth performance [19]. To complete understand the impact of blood transcriptome profiles on milk yield, comprehensive cataloguing of gene expression change within high yielders and low yielders is required. The aim of this study is to compare gene expression profiles in bovine whole blood of high and low milk yield cows, and to investigate potential molecular biomarkers in the blood transcriptome that relate to the productive potential of lactating cows using RNA-seq techniques.

The object of this study is to investigate potential molecular biomarkers in whole blood related to milk production traits in lactating dairy cows, with the aim of putting forward a non-invasive method that identifies the DEGs for milk performance traits. In this study, whole blood genome transcriptome profiles of high and low yield milk cows were investigated using RNA-seq technique. RNA-seq has many advantages over traditional cDNA microarray technologies and it can easily detect low-abundance genes [22, 23]. Marioni et al. demonstrated that RNA-seq and qRT-PCR have a high correlation, and that the Pearson correlation could reach 0.929 [23], which means RNA-Seq is accurate and reproducible.

Among these 100 DEGs, many genes were also detected by other studies associated with milk yield, such as BOLA-DQA1. BOLA-DQA1 is one of the primary histocompatibility complex (MHC) class II molecules that plays an important role in the immune system. It is thought that MHC genes indirectly affect milk production traits by increasing the disease resistance of the cows. In this present study, the expression level of BOLA-DQA1 was higher in high milk yielders than in low milk yielders. Other research has reported that the BOLA-DQA1 gene was associated with resistance to mastitis progression [24].

Splicing event analysis showed that DEGs in high and low milk yield cows were different. In low yielders, cows are more inclined to take the ES pattern but in high yielders A3’SS was more likely. Oxidized low-density lipoprotein receptor (OLR1), which can degrade the oxidized form of low-density lipoprotein and plays important roles in fatty acid transport, was reported as a potential gene for milk-fat percentage and milk fat yield [25]. The splice site analysis in our study also found that AS in OLR1 was significant. There were 5 splice sites of OLR1 in high milk yield including one ES, one A5’SS and three other events. But in low yielders, there was only one splice site on OLR1, which was A5’SS. Also, the expression level of OLR1 was higher in high yielders than in low yielders. Splice events analysis of DEGs can reveal differences between high and low milk yield cows. Among those DEGs containing splice sites, gene ontology enrichment indicated that the PPAR signaling pathway was significantly different between high and low groups. Study showed that PPARA is expressed in heart, liver, adipose tissue, and muscle tissue and is involved in fatty acid catabolism [26]; PPARs play important roles in the regulation of metabolic and inflammatory signaling pathways [26, 27]. It is also reported that PPAR-γ (PPARG) is over-expressed in adipose tissue and macrophages and primarily regulates adipogenesis [27, 28]. PPAR-γ has been reported to significantly increase its expression during lactation in bovine mammary gland [28]. Also, anti-inflammatory properties were observed between PPARA and PPARG [29]. In our study, we found that this PPAR pathway is also in the blood transcriptome, and it was only detected in high yielders. Moreover, these three PPAR pathway DEGs: LPL, OLR1, and CPT1A, were highly expressed in high yielders. Also, the Toll-like receptor signaling pathway existed only in high yielders. The Toll-like receptor pathway is also involved in innate immunity [29]. These results indicate that the regulation of metabolic and immune function is more active in high milk yield cows. The immune system plays a key role in health maintenance, pathogenesis, diseases resistance, and production performance. Further research is required to explore the relation between immune response and milk performance traits in cattle blood transcriptomes.

The present study provides a non-invasive method to identify the DEGs in cattle blood using RNA-seq for milk yield traits. The study revealed 100 DEGs between high yielders and low yielders, and discovered different alternative splicing patterns between the two groups. The enrichment analysis also revealed that specific metabolic and immunological pathways are related to cattle milk yield traits, and could be considered a signature of blood biomarkers selection in dairy cattle. These results provide the valuable resources of biological research in Chinese Holstein cows milk production, but also offer some potential guidelines to understand the relationship between milk production and the immune function.

To further investigate the function of 100 DEGs, we performed the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using DAVID bioinformatics resource tool [21]. GO and KEGG terms with p<0.05 were considered significantly enriched.

A3SS, alternative 3’ splicing site; A5SS, alternative 5’ splicing site; AS, Alternative splicing; BP, biological process; CC, cellular component; DEG, differentially expressed gene; ES, exon skipping; EST, expressed sequence tag; FDR, false discovery rate corrected p values; FPKM, Fragments per kilobase of transcript per million mapped fragments; GO, Gene ontology; IR, intron retention; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; MXE, mutually exclusive exon; RNA-Seq, RNA-Sequencing