Preliminary study on gene regulation and its pathways in Chinese Holstein cows with clinical mastitis caused by Staphylococcus aureus

Bovine mastitis is an inflammatory response of the udder tissue due to physical trauma or microorganism infections, most significant among these being bacterial infections (10). Despite considerable efforts to control bovine mastitis and explain its causes, it remains the most costly and common disease of dairy cattle worldwide, having a negative economic effect on the dairy industry due to reduced milk production, increased treatment costs, reduced fertility, and increased culling of affected animals. Mastitis can be classified as clinical or subclinical in form according to the degree of inflammation (1). The incidence rate of clinical mastitis (CM) ranges from 13%–40% yearly, varying by country and type of cattle housing, and the average cost of CM can amount to $744 for treatment per case (17). The disease can be easily diagnosed based upon visible symptoms, which are udder inflammation showing redness in the affected part or complete udder, warmth, swelling, pain upon touch, milk clots, discoloration, and changes in the consistency of milk (27). Severe cases of CM can be fatal. More importantly, the recurrent nature, a frustrating aspect of CM, brings challenges in the therapeutic and management approaches to bovine CM (17).

Staphylococcus aureus, a common mastitis pathogen widespread in the natural environment of dairy farms, is a frequently isolated bacterium responsible for bovine mastitis (5). This bacterial disease causes a substantial reduction in both milk production and quality, and consequent considerable economic losses in the dairy industry. The bacterium can propagate and persevere in various ways once it adheres to host tissues or prosthetic materials, and expresses many virulence proteins involved in evading the host defences, hence facilitating microbial colonisation of the mammary glands of animals (35). Because of the resistance of S. aureus to antibiotics, infection with this pathogen has remained recalcitrant to eradication. Staphylococcus aureus remains a predominant pathogen causing the most virulent forms of bovine mastitis and poses the greatest challenge to dairy production worldwide.

Sequencing of RNA (RNA-seq) has become a powerful, effective and well-known tool to characterise transcriptome profiles in different tissues or organs. It has also been widely used for exploring differentially expressed genes (DEGs) associated with complex traits and diseases. Published studies on CM have investigated bovine mammary gland tissue (29) or bovine mammary epithelial cells (32) challenged with S. aureus, and have also identified a set of genes in dairy cows with a subclinical mastitis response to S. aureus infection (31). Although the immune response and inflammatory processes are known as key mechanisms for bovine mastitis, limited studies have been carried out on the functional genes and pathways underlying clinical mastitis with S. aureus infection. In this regard, this study aims to gain a deeper knowledge of the host’s transcriptional response to S. aureus infection in developed CM using RNA-seq. The present study will also provide insight into CM pathogenesis with S. aureus infection by identifying some novel functionally relevant genes and biological processes as specific key actors defining the course of the disease.

Staphylococcus aureus isolation and identification. A total of 20 lactating Chinese Holstein dairy cows kept on the private Yajie dairy farm located in Luoyang, China, were recruited to the present study. To detect pathogenic bacteria in cattle with mastitis, milk samples were collected from each cow and examined. Considering Escherichia coli, Streptococcus agalactiae and S. aureus as the most common bacteria associated with mastitis, the three corresponding indexes were used for subsequent detection in this study. In brief, agar with the addition of 10% bovine blood was used for the isolation of pathogens. The samples were inoculated onto plates and incubated at 37°C for 24 h. The grown colonies were identified by morphological (size and colour) and physiological properties (Gram staining and haemolysis). Suspected colonies were then transferred into semisolid slope agar to be identified and subcultured on eosin-methylene blue medium, Group B Streptococcus medium or Baird-Parker agar, as appropriate for the particular colony. Then, triple sugar iron, indole, methyl-red, Voges-Proskauer and citrate utilisation tests were performed to identify E. coli. Catalase, hippurate hydrolysis, Christie-Atkins-Munch-Peterson and sorbitol tests were carried out to identify Streptococcus agalactiae. The identification of S. aureus was carried out based on growth characteristics and biochemical tests according to Hogan et al. (15) and PCR amplification for the S. aureus protein A gene (spa) (3). To investigate whether S. aureus infection affects clinical mastitis, only the milk samples identified as positive for S. aureus exclusively were selected for subsequent research.

Sample collection. Seven dairy cows were selected based on the method described above. Three of them were the mastitic group (M, n = 3), because they had the clinical features of mastitis, such as flakes and clots in milk, udder inflammation (hyperaemia, pain and udder enlargement). Four healthy cows formed the non-mastitic control group (NM, n = 4), because of these cows’ somatic cell counts (SCC) below 200,000 cells/mL at the monthly check-up and freedom from any clinical signs of mastitis.

RNA isolation, cDNA library construction, and RNA sequencing. A 10 mL whole blood sample was collected by venepuncture from each animal. Peripheral blood leukocyte (PBL) enrichment was achieved by depleting the red blood cells and centrifuging for 15 min at 1,000 × g. The leukocytes were mixed with 1 mL of Trizol (TransGen Biotech, Beijing, China), and the total RNA was extracted from the mixture according to the manufacturer’s protocol. RNA quality and quantity were determined by the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) in 1% agarose gel. A total of 4 μg purified total RNA for each sample was used to construct RNA-seq libraries. Overall, the mRNA for each library was enriched using magnetic beads with Oligo (dT) and then fragmented into short segments. Then, these short segments were utilised for synthesising the cDNA with the addition of deoxyribonucleotide triphosphate buffer and DNA polymerase I. The obtained fragments underwent size selection and PCR amplification enrichment, and were subsequently used to construct the library for each sample. Sequencing of each library was performed using the HiSeq 2500 platform (Illumina, Santa Clara, CA, USA).

Sequencing data analysis. For the raw data, Trimmomatic v0.40 software was used to remove the low quality reads, reads containing poly-N sequences, and adaptor sequences (6). After filtering, the clean reads of each library were aligned to the bovine reference genome (ARS-UCD1.2) using HISAT v2.2.1 software (19). The count matrix and transcripts per million values for each gene or transcript were analysed using StringTie v2.1.3 software (20) and the edgeR R package (24), respectively. Differentially expressed gene analysis between the M and NM groups was performed using the DESeq2 software (22). A false discovery rate (FDR) corrected P value ≤ 0.05 and fold change > 2 were defined as the selection threshold for the DEGs. Gene ontology (GO) enrichment analysis (http://www.geneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed for the upregulated and downregulated DEGs using KEGG Orthology-Based Annotation System (KOBAS) software with an adjusted cutoff value of P < 0.01. The STRING database (https://string-db.org) was used to explore the protein– protein interactions (PPI) among the DEGs. In the PPI analysis, the responses of metabolic pathways and platelet activation pathways to CM were considered important and the DEGs located on these pathways were investigated. The PPI network of the DEGs was visualised using Cytoscape v3.8.2 software.

Verification of DEGs by quantitative real-time PCR. Seven randomly selected DEGs (VWF, ITGB3, F2RL3, PRKG1, GDA, HNMT, and GPT2) of biological interest were investigated by quantitative real-time PCR (qRT-PCR) to validate the accuracy of the RNA-seq results. The β-actin gene was set as the reference gene. Primer pairs for VWF (F: 5′-TCGGGATTGGGGACCAGTAT-3′, R: 5′-CAC TGACGAACCCAGAGCAT-3′), ITGB3 (F: 5′-GCT CCCAGATTACCCAAG-3′, R: 5′-TCTTCC ACCTGC CGAACT-3′), F2RL3(F: 5′-CTGCTGCTG TTGCAC TTCTC-3′, R: 5′-CTGAACTCGGCAGACAC-3′), PRKG1 (F: 5′-CTGGAGAAGCGGCTGTCA-3′, R: 5′-TCG TGGAAGGACCTGTAGG-3′), GDA (F: 5′-CCTGGG GCTGGATAGAGAGA-3′, R: 5′-GCCCAC ATACAC CTCCTCAA-3′), HNMT (F: 5′-ACACTG CTGGTC TTATCTCTCC-3′, R: 5′-AGCCCCCTGAAA CCTAGT TG-3′), GPT2 (F: 5′-AGTCCATGAACCCGC AGG-3′, R: 5′-GCACGGATGACCTCAGTGA-3′) and β-actin (F: 5′-ATCACCATCGGCAATGAGCGGTTC-3′, R: 5′-CGG ATGTCGACGTCACACTTCATGA-3′), were designed using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). After amplification, the relative fold change of the differentially expressed genes was calculated using the 2^-ΔΔCT algorithm.

Observations in experimental and control animals and pathogen prevalence in their milk samples. The characteristics of the seven dairy cows are given in Table 1. The results showed that the cows in the M group had higher SCC values (ranging from 456 to 623 × 104/mL) than those in the NM group. Changes in milk appearance (flakes and clots) were observed in the M group. Notably, the cows in the M group had a lower daily milk yield than the cows in the NM group. Samples were positive in 55% (11/20), 20% (4/20) and 30% (6/20) proportions for S. aureus, E. coli, and Streptococcus agalactiae, respectively, and 15% (3/20) were positive for all three pathogens. The three CM cases selected for RNA-seq were positive for S. aureus only (Fig. 1), and negative for E. coli and Streptococcus agalactiae.

Fig. 1

Characterisation of transcriptome sequencing data. Seven RNA sequences were generated from the four healthy and three mastitic dairy cows. The results are listed in Table 2. After the quality control of the raw data, we obtained an average of approximately 47.60 Mbp clean reads for each sample. The Q20 and Q30 of each sample ranged from 97.78% to 97.96% and 91.07% to 91.82%, respectively. By aligning clean reads to the bovine reference genome, we obtained an average mapping ratio of 81.10% (from values ranging from 76.37% to 84.97%). These results suggested that the data quality was suitable for subsequent analysis.

Identification of differentially expressed genes. A total of 16,544 genes were detected. Their expression profiles are shown in Fig. 2A. Of these genes, 4,286 were indicated to be DEGs based on the threshold level(log2 (fold change) ≥ 1 and FDR < 0.05) (Fig. 2B). Among them, 3,085 DEGs were classified as upregulated, and the remaining 1,201 as downregulated. Figure 2C displays the hierarchical cluster expression patterns of the DEGs among all the studied groups.

Fig. 2

Gene ontology/pathway enrichment analysis of DEGs. Investigating the relationship between DEGs and the inflammatory response by GO analysis, we found that most upregulated DEGs were annotated to 549 GO terms including 402 biological processes (BP), 70 cellular components (CC), and 77 molecular functions (MF), while most downregulated DEGs were enriched into 297 GO terms consisting of 225 BP, 39 CC, and 33 MF (Fig. 3).

Fig. 3

For the upregulated DEGs, the top three BP were the cellular process, response to stimulus, and biological regulation. The cellular anatomical entity, intracellular component, and organelle were the three most predominant CC, while binding, catalytic activity, and protein binding were the major three MF. These results suggested that the invasion of the pathogen might induce changes in various membrane structures and cellular components.

For the downregulated DEGs, the same three foremost BP as for the upregulated DEGs were observed. The cellular anatomical entity, membrane, and intrinsic component of the membrane were the three most important CC, while binding, ion binding, and catalytic activity were the leading three MF. These results suggested that this pathogenic infection might be related to dynamic changes in gene expression in specific cellular biological processes.

To further explore the DEGs involved in the pathway of the inflammatory response, the KEGG database was employed to identify key candidate genes. We observed that 37 and 22 pathways were significantly enriched with upregulated and downregulated DEGs, respectively (Fig. 4). For the upregulated DEGs, the top three pathways were platelet activation, metabolic pathways, and pathways in cancer. As for the downregulated DEGs, the three principal pathways were metabolic pathways, pathways in cancer, and Cushing syndrome. These results suggest that both upregulated DEGs in platelet activation and downregulated DEGs in the metabolic pathway might play a vital role in the inflammatory response.

Fig. 4

Differentially expressed genes associated with metabolic pathways and platelet activation. The results of the PPI analysis showed that the downregulated DEGs in the metabolic pathway clustered into four PPIs (Fig. 5A), while the upregulated DEGs in platelet activation clustered into one PPI (Fig. 5B). For the downregulated DEGs, the largest PPI network contained 95 genes, 6 of which had high degrees of connectivity and a threshold for log₂ fold change < −3.0. These genes were AKR1C4, PTGS2, HNMT, EPHX2, CMBL, and IDH1. For the upregulated DEGs, the PPI network consisted of 22 genes, 8 of which had high degrees of connectivity and a threshold for log₂ fold change > 3.0. These genes were VWF, GP9, MYLK, GP6, F2RL3, ITGB3, GP5, and PRKG1.

Fig. 5

Confirmation of gene expression with qPCR. The results of RNA-seq validation in NM to M comparisons showed that the gene expression patterns obtained using sequencing were consistent with the results generated from qRT-PCR (Fig. 6).

Fig. 6

Clinical mastitis is one of the most frequent and costly diseases in dairy cows, causing serious inflammatory diseases and economic losses. Currently, there is no optimal treatment or effective drug for CM. Our understanding of the host immune response to S. aureus infection in dairy cows with CM is limited, making it difficult to design novel therapeutics. A bacterial infection usually causes massive alterations in the host transcriptome, leading to an aberrant host metabolism and a modulated immune response which is ideal for bacterial invasion. In this study, we describe a transcription dataset of PBL samples in dairy cows with CM and in age-matched control samples. A total of 4,286 DEGs were observed in the CM group relative to the control group. Of these, 3,085 genes were expressed relatively more intensely, and the remaining 1,201 were expressed relatively less intensely. The number of DEGs identified here through RNA-seq analysis exceeds the number of DEGs previously reported by Wang et al. (31) for the same RNA sample types (17 DEGs consisting of 7 upregulated and 10 downregulated genes). The reason for this difference is that the mastitic individuals in our study had a more serious disease status. It means that a more comprehensive clinical study and a larger scale of transcriptome profiling with more cases are needed to evaluate the pathogenesis of CM with S. aureus infection.

The functional annotation analysis of DEGs indicated that both upregulated and downregulated genes were assorted into categories including cellular anatomical entity, cellular process, binding and membrane-bounded organelle, which all play important roles in bacterial invasion of the host. The upregulated DEGs were found to mainly be involved in platelet activation and the downregulated DEGs to mainly be involved in the metabolic pathways. Numerous studies have shown that regardless of an inflammatory challenge being transient acute, sustained subacute, or repeated transient, such a challenge can directly affect the metabolic function in lactating dairy cattle (7).

The upregulated DEGs in the platelet activation pathway clustered into one PPI network. We observed several genes displaying high degrees of connectivity. Of these, three upregulated DEGs (VWF, ITGB3, and MYLK) were reported by other authors to be associated with inflammation (16, 18, 33). In particular, evidence revealed that the VWF gene is a plasma glycoprotein with a crucial haemostatic role in serving as a molecular bridge linking platelets to subendothelial components following vascular injury (18). Interestingly, several glycoprotein Ib-IX-V complex genes, including GP9, GP6, GP5, and GP1BA, were observed to be receptors for VWF and ascribed high importance in mediating VWF-dependent platelet adhesion to blood vessels (14). In addition, S. aureus was found to bind directly to VWF using its von Willebrand binding protein (11) and form small bacteria-platelet-fibrin microaggregates, which strongly adhere to VWF exposed on the vessel wall (30). The levels and multimer size of VWFs tend to increase during inflammatory and septic states caused by S. aureus (24). An additional gene with inflammatory activity is F2RL3, because it encodes coagulation factor II (thrombin) receptor-like 3, and its methylation contributes to inflammation through platelet activation (9). Jaeger et al. (16) found that ITGB3 gene expression was increased in porcine mammary epithelial cells after a challenge with S. aureus. In addition, PRKG1, a key regulator of adipokine secretion, has an important regulatory function for milk fatty acid metabolism and may be related to feeding intake and feed conversion efficiency in dairy cattle (28). Milk secreted from CM cows was very thin compared with milk from the control group in our study, suggesting that increased expression of PRKG1 in CM cows might be the cause of milk changes. These results suggested that upregulated DEGs in the platelet activation pathway may play a vital role in the inflammatory response of dairy cows infected with S. aureus.

The downregulated DEGs in the metabolic pathway clustered into four PPI networks. The biggest contained 95 genes, 6 of which had high degrees of connectivity and a threshold for log₂ fold change < −3.0. Two downregulated DEGs (PTGS2 and IDH1) were reported to be associated with inflammation (4, 23). It is well known that AKR1C4 belongs to a family of enzymes that regulate the metabolism of progesterone. Evidence showed that the transcript abundance of AKR1C4 was greater in the endometrium of cows with higher postovulatory progesterone, suggesting that increased expression of AKR1C4 may stimulate vascularisation of the endometrium to support blastocyst growth and embryo implantation (2). Some studies indicated that CM cows were at a higher risk of abortion than cows without mastitis during early lactation (12). The expression of AKR1C4 was significantly decreased in CM cows in this study. Future studies are needed to analyse the molecular mechanisms underlying the regulation of AKR1C4 and their effect upon reproductive performance in CM cows. A new study found that HNMT disruption resulted in a robust increase in brain histamine concentration (34). Histamine is a potent inflammatory mediator commonly associated with allergic reactions, promoting inflammatory and regulatory responses that contribute to pathological processes (8). We can speculate that the decreased expression of HNMT in CM cows may be related to the inflammatory response. The PTGS2 gene, also known as cyclooxygenase-2 (COX-2), induced by inflammatory stimuli, hormones, and growth factors, is a more important source of prostanoid formation in inflammation and proliferative diseases (13). The activity of EPHX2 is considered a major determinant of the bioavailability of anti-inflammatory epoxyeicosatrienoic acid metabolites in the body. Research indicated that the inflammation response was attenuated in EPHX2 knockout mice (21). We can speculate that the decreased expression of EPHX2 in the CM cows may be related to the degree of inflammation in the disease course. These results suggest that the downregulated DEGs in the metabolic pathway might play an important role in the inflammatory response in dairy cows infected with S. aureus.

In conclusion, a total of 3,085 upregulated DEGs and 1,201 downregulated DEGs were found in the CM groups compared to the control. Functional analysis of GO showed that both the upregulated and downregulated DEGs were mainly enriched in the cellular process, biological regulation, and regulation of the biological process. Analysis of the KEGG pathway showed that the upregulated and downregulated DEGs mainly participated in the platelet activation pathway and metabolic pathway, respectively. Several key candidate genes are known to be related to the inflammatory response. These findings improve our understanding of the biological processes of the CM response to S. aureus infection, enabling us to identify specific genes for diagnostic markers and possibly for targeted therapy.