Di Wang,Lei Liu,Serafino M.A.Augustino,Tao Duan,Thomas J.Hall,David E.MacHugh,5,Jinhuan Dou,Yi Zhang,Yachun Wang and Ying Yu*

Abstract Background:Mastitis in dairy cows caused by Staphylococcus aureus is a major problem hindering economic growth in dairy farms worldwide.It is difficult to prevent or eliminate due to its asymptomatic nature and long persistence of infection.Although transcriptomic responses of bovine mammary gland cells to pathogens that cause mastitis have been studied,the common responses of peripheral blood leukocytes to S.aureus infection across two consecutive generations of dairy cattle have not been investigated.

Keywords:Dairy cow,Disease resistance,Mastitis,Peripheral blood leukocyte,Staphylococcus aureus,Transcriptome,Two generations

Background

Bovine mastitis is an inflammation-driven disease of the mammary gland in cows.It normally occurs in response to infection by one of a number of pathogenic microorganisms includingEscherichia coliandStaphylococcus aureus[1,2].It is widely recognized that the high incidence of bovine mastitis,coupled with associated animal welfare problems and the use of antimicrobials,means that this disease is one of the major challenges facing the dairy industry in the twenty-first century[3].S.aureusis an important type of Gram-positive bacteria and is defined as a“contagious pathogen”[4].S.aureuscan adapt to the mammary gland environment and establish subclinical infections[1,5].The symptoms ofS.aureusmastitis are usually less severe than mastitis caused by infection withE.coli,and may even be asymptomatic;however,the infection can persist for long periods and can also exhibit marked resistance to antibiotics[6].S.aureustypically spreads among cattle within dairy herds during the process of milking[7],making prevention,control and elimination particularly difficult.Therefore,bovine mastitis caused by this bacterial pathogen represents a substantial economic burden to dairy production and also poses a significant public health risk through milk consumption.

Global functional genomics technologies,such as RNA sequencing(RNA-Seq),which is enabled by high-throughput sequencing(HTS),are ideally suited for investigating the complex host-pathogen interaction underlying mastitis disease caused byS.aureusinfection and may also provide relevant data for elucidating the molecular mechanisms associated with resilience to mastitis disease.RNA-Seq enables expression quantification analysis of thousands of genes simultaneously[6,8],thereby facilitating identification of individual genes that exhibit the largest expression changes in response to a biological perturbation such as infection by a bacterial pathogen.Several transcriptomics studies have investigated the bovine mammary gland tissue response toS.aureus[9]andE.coli[10],and also the hepatic tissue response toE.coliinfection and lipopolysaccharide(LPS)challenge[11,12].However,to-date,no studies have been carried out on the peripheral blood transcriptome forS.aureusinfection in cows and their offspring simultaneously.

In the current study,using RNA-Seq,we characterized the transcriptomes of peripheral blood leukocyte(PBL)samples collected fromS.aureus-positive motherdaughter pairs andS.aureus-negative control motherdaughter pairs.The main objective of this work was to identify candidate blood-based transcriptional biomarkers for bovine mastitis caused byS.aureusand to provide new insights into host-pathogen interaction and the genetics of disease resistance.

Methods

Animal selection and sampling

All procedures for collection of animal blood and milk samples were approved by the Animal Welfare Committee of China Agricultural University,Beijing,China.All experiments were conducted according to the regulations and guidelines established by this committee(permit number:DK996).

In the present study,a lactating dairy cow and her lactating mother were defined as a mother-daughter pair.A total of 38 mother-daughter pairs were selected based on pedigree information and somatic cell counts(SCC)from a lactating herd of Holstein cows(n≈1,200)at a dairy farm near Beijing,China.The SCC data for three consecutive months were measured and recorded.The Dairy herd improvement(DHI)records were provided by the Dairy Data Centre of China(www.holstein.org.cn).For the present study,milk SCC values of less than 100,000 cells/mL from both mothers and daughters were considered to be healthy[13].Conversely,SCC values larger than 100,000 cells/mL were considered for diagnosis of subclinical mastitis.Detailed information on cow selection is provided in Fig.1a and Table 1.

ForS.aureusidentification,a total of 30 mL of fresh milk was collected and mixed from all of the four lactating quarters of each cow.Subsequently,based on the bacteriological culture and PCR results,fourS.aureuspositive mother-daughter pairs were selected and separated into theS.aureusmother group(SM,n=4)and theS.aureusdaughter group(SMD,n=4).Similarly,fourS.aureus-negative mother-daughter pairs were separated into the healthy mother group(HM,n=4)and the healthy daughter group(HMD,n=4)(Fig.1a).

S.aureus isolation and identification

Bacteriological culture of milk samples was carried out according to National Mastitis Council standards[14].A volume of 3 mL milk was mixed into trypticase soy broth containing 7.5% NaCl and cultured at 37°C for 18-24 h.After that,a total of 10 mL culture was placed into Baird-Parker agar plates with tellurite and 5% egg yolk and cultured at 37°C for 18-24 h.Two suspected colonies from each sample with surrounding clear zones were transferred to trypticase soy agar plated for DNA collection.Following this step,PCR amplification and sequencing of theS.aureusthermonuclease gene(nuc)[15]was performed(Fig.1b and c).The PCR reaction was performed in 25μL,containing 3μL of genomic DNA(30-50 ng/μL),1μL of each primer(10μmol),12.5μL of Taq?Mix(1.25 units/25μL reaction)and 7.5μL of ddH2O.PCR was performed using the following thermocycler program:94°C for 10 min;35 cycles of 94°C for 30 s,59°C for 30 s and 72°C for 30 s;72°C for 7 min.

Blood collection,RNA extraction and RNA-sequencing

A 20-mL blood sample was obtained from the caudal vein from each animal for buffy coat(leukocytes)collection using 15 min centrifugation at 3,000 r/min.

TRIzol reagent(Invitrogen,Carlsbad,CA,USA)was used to isolate total RNA from leukocytes according to the manufacturer’s protocol.The RNA quality was checked on a 1% agarose gel and quantified using a Qubit RNA Assay Kit and a Qubit 2.0 Fluorometer(ThermoFisher Scientific,Waltham,MA USA).RNA integrity was assessed with the BioAnalyzer 2100 System(Agilent Technologies,Santa Clara,CA,USA).All 16 RNA samples had an RNA integrity number(RIN)larger than 7.0(Supplementary Table S1).The 28S:18S rRNA ratios of all samples were larger than 1.7(Supplementary Table S1).An equivalent amount(4μg)of total PBL RNA purified from each animal was used to construct RNA-Seq libraries with the NEBNext? Ultra?RNA Library Prep Kit for Illumina?(NEB,Ipswich,MA,USA).Finally,the libraries were sequenced using 150 bp paired-end reads with the Illumina HiSeq X Ten System(CapitalBio Technology,Beijing,China).

Quality control for raw data

Trimmomatic software version 0.38[16]was used to filter out the adapter sequence and low-quality bases/reads with the default parameters(http://www.usadellab.org/cms/?page=trimmomatic).Further quality assessment of the sequence reads was then undertaken using FastQC version 0.11.8[17].After these QC procedures were completed,the sequence read data were used for the downstream analyses in the computational workflow.

Reads alignment and abundance estimation

TheBos taurusARS-UCD1.2 reference assembly(FASTA format)and annotated gene model(GTF format)were downloaded from the Ensembl database(ftp.ensemble.org/pub/release-96/gtf/bos_taurus).The QC-assessed sequence read data for each sample were aligned to the reference genome using STAR with the basic options[18].Transcript abundance was quantified using featureCounts[19]in the R subread package[20]under the default setting,and read counts were calculated to estimate the transcript expression levels.Genes with read counts more than ten in at least two samples were defined as expressed genes.The rlog-normalized read count was calculated with DESeq2[21],which was then used to performed differential expression analysis.

The transcriptional responses toS.aureusinfection were investigated by comparing differentially expressed genes between the infected and non-infected healthy control groups.We firstly analyzed gene expression levels in the S group(n=8,4 animals each in the SM and SMD groups)and H group(n=8,4 animals each in the HM and HMD groups),and compared the differentially expressed genes(DEGs)between the two groups.Detection of DEGs was also performed for the comparisons of the mother groups(SM vs.HM)and the daughter groups(SMD vs.HMD).

Functional enrichment and annotation

Gene set enrichment analysis was performed for all the expressed genes detected between the S and H groups,using the GSEA software package(v4.0.3)with the Human NCBI Gene ID Molecular Signature Database(MSigDB)version 7.0 and Hallmark database version 7.0[22-24].Gene sets were considered significant whenP＜0.05 andFDR＜0.25[25,26].

Fig.1 Identification of S.aureus isolated from milk samples from Holstein cows.a Workflow of sample selection.A total of 38 mother-daughter pairs were selected from 1,200 lactating dairy cattle according to pedigree and SCC records.After S.aureus identification from milk samples,four S.aureus-positive and four S.aureus-negative pairs were used for subsequent studies.b Specific PCR and electrophoresis map of the S.aureus thermonuclease(nuc)gene(279 bp).SA:abbreviation of S.aureus.+:S.aureus-positive control;-:S.aureus-negative control.SM:S.aureus mastitis mother;HM:healthy mother;SMD:S.aureus mastitis daughter;HMD:healthy daughter.c Sequence alignment of partial nuc gene.The upper lines were amplified from samples collected from the cows’milk;the lower line is nuc sequence amplified from one of the S.aureus-positive samples.Identity=95.68%

The DEGs(P＜0.05 and|FC|＞2)in the comparisons of S vs.H,SM vs.HM and SMD vs.HMD were used to perform the following functional enrichment analysis.Kyoto Encyclopedia of Genes and Genomes(KEGG)pathway and Gene Ontology(GO)enrichment analyses of gene clusters were performed using clusterProfiler[27],with results exhibitingP＜0.05 considered significant.An interaction network analysis was also conducted using Ingenuity Pathway Analysis(IPA;Qiagen,Redwood City,CA,USA)[28].

Bacterium S.aureus S.aureus S.aureus S.aureus S.aureus---------1month 3.65S.aureus 3.02 2.91 3.30S.aureus 3.22±0.14 3.24 3.13 3.32S.aureus 2.88 3.14±0.08 2.92 3.13 3.23 3.12 3.1±0.06 3.08 3.46 3.18 3.18 3.23±0.07-2months 3.75 2.65 2.98 3.29 3.17±0.20 3.11 3.05 3.06 2.86 3.02±0.05 2.97 3.17 3.22 3.25 3.15±0.05 3.18 3.32 3.14 3.22 3.22±0.03 PP -3months 3.63 3.49 3.38 3.42 3.48±0.05 2.96 3.03 3.31 2.83 3.03±0.088 3.06 3.12 3.34 3.88 3.35±0.16 3.01 3.27 3.17 3.11 3.14±0.047-1month 4.28 3.59 4.10 3.52 3.87±0.16 3.51 4.15 2.85 1.73 3.06±0.45 2.88 3.89 3.46 3.95 3.55±0.21 3.53 5.11 2.60 4.33 3.89±0.47-2months 4.06 2.75 3.97 3.81 3.65±0.26 3.46 3.60 3.26 4.95 3.82±0.33 3.02 4.40 3.81 4.31 3.89±0.27 3.66 2.57 3.30 3.76 3.32±0.23 FP -3months 3.93 3.79 2.14 2.10 2.99±0.44 3.67 3.84 2.76 2.89 3.29±0.24 2.78 3.82 3.78 4.69 3.77±0.34 3.91 3.42 3.31 3.48 3.53±0.11-1month 26.80 34.40 35.80 24.60 2.39 41.60 35.60 24.70 37.70 34.9±3.14 38.00 30.10 56.00 38.20 40.58±4.74 38.80 25.10 35.80 36.60 34.08±2.65-2months 28.60 42.60 37.50 31.30 35±2.72 30.4±50.90 37.70 25.50 38.60 38.18±4.49 47.40 35.30 58.80 40.80 45.58±4.38 40.30 46.30 42.70 42.20 42.88±1.09 Table 1 Basic information and bacterial culture of the 16 Chinese Holstein milk samples MY -3months 28.80 31.00 36.70 31.30 31.95±1.45 58.00 42.70 25.70 31.30 39.43±6.18 51.80 38.10 58.30 42.10 47.58±3.97 35.70 52.00 42.90 37.30 41.98±3.19-1month 6.42 8.42 6.09 4.27 6.30±0.74 7.49 4.01 2.96 3.54 4.50±0.88 2.55 2.19 1.64 1.82 2.05±0.17 1.11 5.36 2.29 1.91 2.67±0.81-2months 6.19 8.54 6.43 5.55 6.67±0.56 5.58 3.28 1.49 3.18 3.38±0.73 3.63 1.57 0.88 1.44 1.88±0.52 0.60 1.36 0.68 2.29 1.23±0.34 SCS-3months 5.49 6.28 4.93 3.63 5.08±0.48 4.05 3.21 2.99 2.99 0.22 3.01 1.16 2.60 3.61 0.45 0.82 0.82 2.26 1.26 1.29±0.30-1month 1071.00 4289.00 852.00 241.00 1613.25±787.25 2249.00 202.00 64.00 145.00 73.00 57.00 39.00 44.00 53.25±6.58 2.60±27.00 514.00 61.00 47.00 162.25±101.72 SCC×1000/mL -2months 911.00 4641.00 1076.00 585.00 1803.25±823.94 596.00 121.00 102.00 113.00 233±104.84 665±457.92 3.31±155.00 37.00 23.00 34.00 62.25±26.90 19.00 32.00 20.00 61.00-3months 561.00 972.00 382.00 155.00 517.5±149.63 207.00 116.00 283.00 99.00 176.25±37.03 101.00 28.00 76.00 153.00 89.5±22.54 22.00 22.00 60.00 30.00 33.5±7.82 33±8.48 Sample Age Days of lactation 482 266 529 291 171 299 125 33 258 222 173 244 145 212 271 250 8 5 8 5 6 3 3 2 7 6 5 5 3 4 3 3 SM1 SM2 SM3 SM4 Mean±SD SMD1 SMD2 SMD3 SMD4 Mean±SD HM1 HM2 HM3 HM4 Mean±SD HMD1 HMD2 HMD3 HMD4 Mean±SD Note:SM-S.aureus infected mastitis mothers;SMD-S.aureus infected daughters;HM-Healthy mothers;HMD-Healthy daughters;the four SMD animals were the descendants of the four SM animals,respectively;the four HMD animals were the descendants of the four HM animals,respectively;except HMD1 and HMD3 are half-siblings of one sire,the other animals are descendants of different sires;SCC:Somatic cell counts per milliliter of milk sample;SCS=log2(SCC/100,000)+3;MY:milk yield,FP:milk fat percentage;PP:milk protein percentage;-indicates samples without S.aureus infection

Validation of RNA-Seq results through reverse transcription quantitative real-time PCR(RT-qPCR)

To quantitatively assess the reliability of our sequencing data,the expression levels of seven(SM vs.HM comparison)and ten(SMD vs.HMD comparison)randomly selected genes were estimated with RTqPCR in the same samples used for RNA-Seq(Supplementary Table S1).The RNA samples were reverse transcribed into cDNA using the PrimeScript?RT reagent kit according to the manufacturer’s instructions(Takara Bio,Dalian,China).Real-time detection of specific PCR products was performed with the Sybr Green I Master Mix Kit(Roche Diagnostics,Mannheim,Germany)on the LightCycler 480 II(Roche Diagnostics Ltd.,Basel,Swetzerland),according to the manufacturer’s protocol.The oligonucleotide primers used for the RT-qPCR analysis are provided in Supplementary Table S2.The thermocycler program used was as follows:one cycle of preincubation at 95°C for 10 min,45 cycles of amplification(95°C for 10 s,60°C for 10 s,and 72°C for 10 s).Duplicate RT-qPCR assays were performed on each cDNA sample and relative expression values were calculated using the 2-△△Ctmethod with the bovine glyceraldehyde-3-phosphate dehydrogenase gene(GAPDH)as the internal reference[29,30].The log2|FC|values were calculated from RT-qPCR data to directly compare to the RNA-Seq results.

Statistical analysis

A linear regression analysis was conducted using the Prism software package(version 8;GraphPad Software,San Diego,CA,USA)to evaluate gene expression levels among experimental groups.For RT-qPCR data,the significant differences between groups were examined with the Student’st-test.

Results

Transcriptomic data and detection of expressed genes in peripheral blood leukocyte samples

Fig.2 Comparison of gene expression profiles between mothers and daughters.a Venn diagram showing genes only expressed in S.aureuspositive mother-daughter pairs(purple circle),genes only expressed in S.aureus-negative mother-daughter pairs(yellow circle),and genes in intersection are common to both groups.b Linear regression analysis of expression levels for 201 genes exclusively expressed in the SM and SMD groups.c Linear regression analysis of expression levels for 162 genes exclusively expressed in the HM and HMD groups.d Significantly enriched(P＜0.05)KEGG pathways of genes exclusively expressed in the S.aureus-positive groups(SM and SMD).e Significantly enriched(P＜0.05)KEGG pathways of genes exclusively expressed in the S.aureus-negative groups(HM and HMD)

In the present study,an average of 22,435,225(ranging from 18,801,274 to 25,620,906)paired-end reads were generated as raw data from the 16 RNA-Seq libraries.After filtering,an average of 20,955,429 reads(ranging from 16,864,277 to 23,829,973)were retained for each library.The average uniquely mapping rate was 94.11%(Supplementary Table S3).There were 10,287 detectable genes expressed in both the SM and SMD groups,and 10,248 gene expressed in the HM and HMD groups(Fig.2a).A total of 10,086 detectable expressed genes were shared across the four groups.A total of 201 genes were identified that were exclusively expressed in the two consecutive generations of animals naturally infected withS.aureus(SM and SMD),which may be associated with susceptibility toS.aureus-induced mastitis.Conversely,162 genes were exclusively expressed in the two consecutive generations of healthy control animals(HM and HMD),which may be associated with resilience toS.aureus-induced mastitis.

Furthermore,the correlation between the gene expression profiles(log2|FC|of 4000 randomly selected genes)in the two generations of SM vs.HM and SMD vs.HMD was significant(R2=0.3410,P＜0.0001,Supplementary Fig.S1).This result indicates that gene expression changes due toS.aureusinfection challenge were moderately conserved between the two generations.The linear regression analyses were also performed using the normalized read counts of the 201 exclusively expressed genes(R2=0.9991,P＜0.0001,Fig.2b)in the two generations of animals naturally infected withS.aureus(SM and SMD)and the 162 exclusively expressed genes(R2=0.9921,P＜0.0001,Fig.2c)in the two generations of healthy controls(HM and HMD).

Transcriptome changes in peripheral blood leukocytes from S.aureus-positive cattle compared to S.aureusnegative control cattle

It is well established thatS.aureus,as an intracellular pathogen of mammals,has evolved a wide range of mechanisms for immunosuppression and immunoreaction[31,32].To investigate host-pathogen interaction in bovine PBL infected withS.aureus,the transcriptomes of infected cattle and non-infected healthy control animals were compared.Although no significant DEGs were detected using an FDRPadj.threshold of 0.05,a total of 17,13 and 10 DEGs were observed for the three comparisons(Table 2,3 and 4),respectively,using the criteria of FDRPadj.＜0.1 and|FC|＞1.2[33-35].Furthermore,a total of 301,283 and 260 DEGs(P＜0.05,log2|FC|＞2)(Supplementary Table S5,6 and 7 and Fig.S2)were used for functional enrichment analysis.The expression heat maps of the 301,283 and 260 genes of the three comparisons revealed different transcriptional profiles betweenS.aureus-positive andS.aureus-negative samples(Fig.3).To validate the RNA-Seq results,seven and ten genes were selected for RT-qPCR validation in the SM vs.HM and SMD vs.HMD comparisons(Fig.4),respectively.Results from this analysis showed that the gene expression patterns obtained using RNASeq were consistent with the results generated from RTqPCR(Fig.4a and c).The correlation coefficients between the RNA-Seq and RT-qPCR results for the genes in the SM vs.HM and SMD vs.HMD comparisons were 0.97 and 0.90(P＜0.0001),respectively(Fig.4b and d).

Table 2 Differentially expressed genes in the S vs.H comparison(FDR Padj.＜0.1 with|FC|＞1.2)

Table 3 Differentially expressed genes in the SM vs.HM comparison(FDR Padj.＜0.1 and|FC|＞1.2)

Functional enrichment and annotation

The GSEA analysis performed between S and H groups revealed that a total of 38 gene sets were upregulated in S group compared to the controls,of which six gene sets were significant atFDR＜0.25 andP＜0.05.Significantly enriched gene sets(FDR＜0.25 andP＜0.05)in the S group are shown in Supplementary Fig.S3 and the enrichment information is summarized in Supplementary Table S4.

KEGG pathway analysis of 201 genes exclusively expressed inS.aureus-positive animals revealed nine significantly enriched KEGG pathways(P＜0.05)(Fig.2d).Whereas genes exclusively expressed inS.aureusnegative animals revealed four significantly enriched KEGG pathways(P＜0.05)(Fig.2e).Five genes(EPOR,IL9,IFNL3,CCL26,andIL26)exclusively expressed inS.aureus-negative animals were enriched in the cytokinecytokine receptor interaction KEGG pathway.

KEGG pathway analysis of the 301 DEGs in the S vs.H comparison revealed 15 significant enriched pathways(P＜0.05)(Fig.5a).Most of these pathways were associated with the immune response or inflammation.For example,the IL-17 signaling pathway was enriched by three DEGs(CCL20,IL13,andMMP3).Furthermore,ten KEGG pathways were significantly enriched(P＜0.05)by DEGs in the SM vs.HM comparison and,among these,S.aureusinfection was the most significantly enriched pathway(Fig.5b).Moreover,DEGs in the comparison of SMD vs.HMD were significantly(P＜0.05)enriched in ten pathways,includingS.aureusinfection,phagosome.Importantly,six DEGs were significantly involved in theS.aureusinfection pathway(Fig.5c).The top ten biological processes identified using GO enrichment analysis are shown in Supplementary Fig.S4.

To construct interaction networks,the DEGs in each comparison were further analyzed using the IPA softwaretool.Interestingly,the PI3K family is shared across all networks generated for the three comparisons and interacted with 13,6 and 11 molecules in comparisons of S vs.H,SM vs.HM and SMD vs.HMD,respectively(Figs.6,7 and 8).Moreover,the IL12 complex and kinase AKT both play key roles in the SM vs.HM and SMD vs.HMD comparison(Figs.7 and 8).Collectively,these functional analysis results suggested that the immune responses were activated in PBL fromS.aureus-infected dairy cattle.

Table 4 Differentially expressed genes in the SMD vs.HMD comparison(FDR Padj.＜0.1 and|FC|＞1.2)

Fig.3 Transcriptomic changes of PBL in the three comparisons.Cluster analyses and heat maps of DEGs in(a)S vs.H,(b)SM vs.HM and(c)SMD vs.HMD.The different columns represent different samples and the different rows denote different DEGs.Red and blue show increased and decreased expression,respectively

Fig.4 RT-qPCR validation of RNA-Seq results.a and c RT-qPCR confirmation results for the randomly selected DEGs from the SM vs.HM comparison and the SMD vs.HMD comparison,respectively.b and d Regression analysis of the log2|FC|values between the RNA-Seq and RT-qPCR validation of the SM vs.HM and the SMD vs.HMD comparisons,respectively.The GAPDH gene was used as an internal reference control gene.*means P＜0.05.**means P＜0.01

Fig.5 KEGG pathway enrichment of DEGs.a KEGG pathways of DEGs in the comparison of S vs.H and expression differences of the DEGs involved in the enrichment of the IL-17 signaling pathway.b KEGG pathways of DEGs in the comparison of SM vs.HM and expression differences of the DEGs involved in the enrichment of the S.aureus infection pathway.c KEGG pathways of DEGs in the comparison of SMD vs.HMD and expression differences of the DEGs involved in the enrichment of the S.aureus infection pathway

Finally,the linear regression results between normalized read counts and SCC showed that significant correlations were observed between the expression levels of DEGs(CXCL9,SOCS1,LOC508858,andCYP2E1)and SCC(Fig.9),suggesting these genes might be served as potential molecular biomarkers of mastitis caused by infection withS.aureus.

Discussion

Fig.6 IPA interaction network generated from the DEG list in the S vs.H comparison.Red denotes gene nodes that were increased in expression,and green denotes nodes that were decreased in expression in S samples compared to H samples.The red circle denotes a commonly shared node in the three comparisons

In recent years,healt h traits,including mastitis resistance,are increasingly being incorporated into the breeding goals of the modern dairy industry[36-38].Mastitis is a low heritability trait[39,40]and it is therefore difficult to achieve significant improvements using conventional breeding strategies.Therefore,it is important to investigate the complex host-pathogen interaction underlying mastitis disease caused byS.aureusinfection from the perspective of the peripheral blood transcriptome.In previous genome-wide association studies(GWAS),a number of candidate genes(TRAPPC9[41],mTORC1[42],JAK2andSTAT5A[43])were observed to be associated with mastitis-related traits such as SCC.However,transcriptome profiles ofS.aureusinfection in both cows and their offspring have not been studied in detail.In the current study,for the first time,transcriptome profiles of PBL in response toS.aureusinfection were investigated in two generations of dairy cattle.

Previous studies have identified many genes as potential expression biomarkers for bovine mastitis caused by infection withS.aureusmastitis[6,44].However,the concordance among these results is not high,possibly as a consequence of environmental factors and varying genetic backgrounds.To address this issue,eight mother-daughter pairs(fourS.aureus-infected pairs and four uninfected mother-daughter pairs)from the same farm were used as experimental subjects for the current study.Three well-defined animal selection criteria were used for comparative analyses of the transcriptomes of PBL sampled from two generations of dairy cattle infected withS.aureus-induced mastitis:1)availability of mother-daughter pairs,2)three consecutive months of SCC score recordings,and 3)unambiguousS.aureusdetection and identification.

Fig.7 IPA interaction network generated from the DEG list in the SM vs.HM comparison.Red denotes gene nodes that were increased in expression,and green denotes nodes that were decreased in expression in SM samples compared to HM samples.The red circle denotes a commonly shared node in the three comparisons

For dairy cattle,to produce 1 kg of milk,400 to 500 kg of blood must pass through the udder[45].Consequently,PBL represents a valuable sample resource to evaluate the condition of udder health.In our previous study[44],we infected udder quarters withS.aureus(1×106CFU/mL),then conducted mammary biopsies,and collected udder tissues for RNA extraction and transcriptome analysis using RNA-Seq.Mammary biopsies cause trauma and stress;therefore,the use of an accessible and informative tissue such as peripheral blood from naturally infected animals in a production population substantially mitigates animal welfare issues.Thus,PBL samples from two generations ofS.aureus-infected cows were used to perform transcriptome sequencing.

Fig.8 IPA interaction network generated from the DEG list in the SMD vs.HMD comparison.Red denotes gene nodes that were increased in expression,and green denotes nodes that were decreased in expression in SMD samples than HMD samples.The red circle denotes a commonly shared node in the three comparisons

Our study is the first to describe transcriptomic profiles of PBL inS.aureus-infected and uninfected dairy cattle across two generations.More than 160 genes were observed solely in the healthy mother-daughter pairs compared to theS.aureus-infected mother-daughter pairs,an observation that may be relevant to identification of genes underpinning resilience toS.aureusmastitis(Fig.2a).These genes were observed to be significantly enriched in four biological pathways with the cytokine-cytokine receptor interaction pathway directly related to the immune response.Five genes(EPOR,IL9,IFNL3,CCL26,andIL26)that were involved in this pathway might serve as potential molecular markers for breeding programs that enhance resistance toS.aureusinfection and mastitis(Fig.2e).Thus,these immunerelated genes warrant further study in a larger group of animals across two generations.

EPORencodes the erythropoietin receptor,which is a member of the cytokine receptor family.It can activate JAK2 tyrosine kinase that,in turn,activates a range of intracellular pathways.A previous study observed that theIL9gene emerged as a central node in the network associated withStreptococcus agalactiae-induced subclinical mastitis[46].In this regard,we have also identifiedIL9as a key gene associated withS.aureus-induced subclinical mastitis.Interestingly,sequence variation atIFNL3has been shown to be associated with hepatitis B virus(HBV)infection in a Chinese human population[47].CCL26encodes a secreted protein involved in immunoregulatory and inflammatory processes,which Anna et al.observed to be associated with IL4-mediated signaling pathways in bovine macrophages challengedin vitrowithS.agalactiae[48].The protein encoded byIL26is a member of the IL10 family of cytokines[49],which has been shown to exhibit anti-bacterial activity against a wide range of bacteria includingS.aureus[50],suggesting a broader role in host defenses against bacteria[51].

Fig.9 Linear regression analysis between normalized read counts and SCC.The X-axis indicates normalized read counts of DEGs,and the Y-axis indicates SCC.a-d)represent CXCL9,SOCS1,LOC508858,and CYP2E1,respectively

As expected,a small number of DEGs(Fig.3)were detected using FDRPadj.＜0.1,a less stringent statistical threshold.Chinese dairy cattle represent a relatively outbred animal population with significant variation in genetic background[52].S.aureuscan readily adapt to its host through evasion of almost every facet of the immune system[32];therefore,S.aureus-infected cows may exhibit less severe and milder immune responses,which may account for the relatively small number of DEGs.In this regard,Fang et al.also reported modest differential gene expression in bovine mammary gland induced by a low dose ofS.aureusto simulate naturally infected mastitis[44].It would therefore be instructive to perform a substantially largerin vitrochallenge experiment using varyingS.aureuschallenge doses.

Indoleamine 2,3-dioxygenase 2 encoded byIDO2is an immunomodulatory molecule with potential effects on various diseases including cancer and autoimmune conditions.A previous study demonstrated thatIDO2functions as a modifier in B cells to control pathogenic inflammation and autoimmunity[53].NR4A1,which encodes the nuclear receptor subfamily 4 group A member 1 protein,was differentially expressed in the comparison of SM vs.HM for the current study.NR4A1has been shown to play a key role as a regulator of the immune response to apoptotic cells[54].

In the present study,immune-related pathways were enriched by different but complementary functional enrichment and annotation methods(KEGG,IPA and GSEA).Three DEGs(CCL20,IL13,andMMP3)were found to be involved in the IL-17 signaling pathway in the S vs.H comparison(Fig.5a).The interleukin 17(IL-17)family are proinflammatory cytokines that play crucial roles in both acute and chronic inflammatory responses[55].Moreover,there are five(BLA-DQB,C1R,C2,FCGR1A,andKRT10)(Fig.5b)and six genes(BLA-DQB,C3AR1,CFI,FCAR,FCGR3A,andLOC10498484)(Fig.5c)enriched in theS.aureusinfection pathway in the SM and SMD groups,respectively.TheBLA-DQBgene encodes the bovine major histocompatibility complex,class I,DQ beta protein,which has been hypothesized to play a key role in the synchronization of the immune response to mastitis-causing bacteria in bovine udder secretory tissues[56].Interestingly,BLA-DQBexhibited contrary directions of expression in SM and SMD compared with the healthy controls,which might be due to the differences in the duration ofS.aureusinfection[57],differences in age[58],parity[59]and lactation period[35]of the cattle.

The IPA analysis revealed an interaction network involving the phosphoinositide 3-kinase(PI3K)family in all three comparisons(Fig.6,7 and 8).The PI3K family of enzymes are involved in myriad cellular functions,including cell growth,proliferation,differentiation,motility,survival,and intracellular trafficking.In particular,PI3K proteins have important functions in immunobiology[60].IL12(Fig.7 and 8)encodes an important cytokine that can regulate both innate and adaptive immune responses during infection[61].Kinase AKT(Fig.7 and 8)is involved in the regulation of the development and function of innate immune cells,including neutrophils,macrophages,and dendritic cells[62].

The cytokine-cytokine receptor signaling pathway genes exclusively expressed in uninfected mother-daughter pairs(EPOR,IL9,IFNL3,CCL26,andIL26)and the DEGs involved in immune-related pathways identified in the current study(includingCCL20,IL13,MMP3,andBLADQB),may contain sequence polymorphisms that could serve as genetic markers for susceptibility to mastitis caused byS.aureusinfection.Consequently,identification of SNPs within these genes and systematic evaluation through large-scale association studies using suitable animal cohorts can provide valuable information for future genome-enabled breeding programs to enhance mastitis resilience in dairy cattle.

Conclusions

In summary,the transcriptome profiles of PBL sampled from two consecutive generations of cows with naturally infectedS.aureusmastitis and non-infected healthy control animals were investigated.Many of the genes highlighted in this study as being important in the bovine host response to infection withS.aureusmay represent candidate PBL expression biomarkers for mastitis and may also contain sequence variation that can be leveraged for genomic selection of cattle less susceptible to mastitis disease.

Supplementary information

Supplementary informationaccompanies this paper at https://doi.org/10.1186/s40104-020-00494-7.

Additional file 1:Supplementary Figure S1.Linear regression analysis of the log2|FC|of 4000 randomly selected gene expression values for the SM vs.HM and SMD vs.HMD comparisons.

Additional file 2:Supplementary Figur S2.Volcano plots of DEGs for(A)the S vs.H comparison,(B)the SM vs.HM comparison,and(C)the SMD vs.HMD comparison.

Additional file 3:Supplementary Figure S3.Gene set enrichment analysis(GSEA)plots depicting the enrichment of functional gene sets up-regulated in the S group compared to the control group(FDR＜0.25 andP＜0.05).(A-F):GSEA plot depicting the enrichment of functional gene sets up-regulated in oxidative phosphorylation,heme metabolism,fatty acid metabolism,Kras signaling up-regulation,inflammatory response and interferon-gamma response in the S vs.H comparison(FDR＜0.25 andP＜0.05).

Additional file 4:Supplementary Figure S4.GO Biological Process enrichment of DEGs.(A)Top ten Biological Process GO categories of DEGs in the S vs.H comparison.(B)Top ten Biological Process GO categories of DEGs in the SM vs.HM comparison.(C)Top ten Biological Process GO categories of DEGs in the SMD vs.HMD comparison.

Additional file 5:Table S1.Quality report of 16 RNA samples for RNA sequencing.

Additional file 6:Table S2.Primer pairs of DEGs used for qRT-PCR validation.

Additional file 7:Table S3.Summary of sequence reads aligenment

Additional file 8:Table S4.Six significantly enriched upregulated gene sets.

Additional file 9:Table S5.Differentially expressed genes detected in S vs.H

Additional file 10:Table S6.Differentially expressed genes detected in SM vs.HM.

Additional file 11:Table S7.Differentially expressed genes detected in SMD vs.HMD.

Abbreviations

S.aureus:Staphylococcus aureus;S:S.aureus-infected cows;H:Healthy noninfected cows;SM:S.aureus-infected mastitis mothers;HM:Healthy noninfected mothers;SMD:S.aureus-infected daughters;HMD:Healthy noninfected daughters;EPOR:Erythropoietin receptor;IL9:Interleukin 9;IFNL3:Interferon lambda 3;CCL26:C-C motif chemokine ligand 26;IL26:Interleukin 26;DEGs:Differentially expressed genes;FDR:False discovery rate;CCL20:Chemokine(C-C motif)ligand 20;IL13:Interleukin 13;MMP3:Matrix metalloproteinase-3;IL-17:Interleukin 17;BLA-DQB:Bos taurusmajor histocompatibility complex,class II,QB beta;C1R:Complement C1r subcomponent;C2:Complement C2;FCGR1A:Fc fragment of IgG receptor Ia;KRT10:Keratin 10;C3AR1:Complement component 3a receptor 1;CFI:Complement factor I;FCAR:Fc fragment of IgA receptor;FCGR3A:Fc fragment of IgG,low affinity IIIa,receptor;E.coli:Escherichia coli;RNASeq:RNA sequencing;HTS:High-throughput sequencing;LPS:Lipopolysaccharide;PBL:Peripheral blood leukocyte;SCC:Somatic cell counts;DHI:Dairy herd improvement;Nuc:Thermonuclease;PCR:Polymerase chain reaction;RIN:RNA integrity number;GSEA:Gene set enrichment analysis;MSigDB:Molecular signature database;KEGG:Kyoto Encyclopedia of Genes and Genomes;GO:Gene ontology;IPA:Ingenuity Pathway Analysis;RT-qPCR:Reverse transcription quantitative real-time PCR;GAPDH:Glyceraldehyde-3-phosphate dehydrogenase;PI3K:Phosphoinositide 3-kinases;IL12:Interleukin 12;AKT:Protein kinase B,also known as PKB;CXCL9:Chemokine(C-X-C motif)ligand 9;SOCS1:Suppressor of cytokine signaling 1;CYP2E1:Cytochrome P450 2E1;GWAS:Genome-wide association study;TRAPPC9:Trafficking protein particle complex 9;mTORC1:Cytosolic arginine sensor for mTORC1 subunit 1;JAK2:Janus kinase 2;STAT5A:Signal transducer and activator of transcription 5A;IDO2:Indoleamine 2,3-dioxygenase 2;NR4A1:Nuclear receptor subfamily 4 group A member 1

Acknowledgments

The authors thank three reviewers for constructive comments and suggestions on the manuscript,Siyuan Mi,Liangyu Shi,Jiangang Qi and all the members of Molecular and Quantitative Genetics Laboratory(China Agricultural University)for sample collection;CapitalBio Technology(Beijing,China)for sequencing services;Gillian P.McHugo for helpful discussions.The support provided by the China Scholarship Council(201906350210)during a visit by Di Wang to University College Dublin is acknowledged.

Authors’contributions

Ying Yu and Di Wang conceived the study.Di Wang,Lei Liu,and Tao Duan collected the samples.Di Wang,Lei Liu,and Thomas J.Hall analyzed the data.Di Wang,Lei Liu,Serafino M.A.Augustino,Tao Duan,David E.MacHugh,Jinhuan Dou,Yi Zhang,Yachun Wang,and Ying Yu wrote and prepared the manuscript.All authors reviewed and contributed to the manuscript.

Funding

This research was financially supported by the NSFC-PSF Joint Project(31961143009),Beijing Natural Science Foundation(6182021),Beijing Dairy Industry Innovation Team(BAIC06),Modern Agro-industry Technology Research System(CARS-36)and the Program for Changjiang Scholar and Innovation Research Team in University(IRT-15R62).

Availability of data and materials

All genomic annotation data defining gene regions are available for download(ftp://ftp.ensemble.org/pub/release-96/gtf/bos_taurus).RNA-Seq data from China Agricultural University is available upon the agreement of China Agricultural University and should be requested directly from the authors.

Ethics approval and consent to participate

All animal management and experimental procedures followed the animal care protocols approved by the China Agricultural University Animal Care and Use Ethics Committee.

Consent for publication

Not applicable.

Competing interests

The authors have declared that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Author details

1Key Laboratory of Animal Genetics,Breeding and Reproduction,Ministry of Agriculture & National Engineering Laboratory for Animal Breeding,College of Animal Science and Technology,China Agricultural University,Beijing 100193,China.2UCD School of Agriculture and Food Science,University College Dublin,Dublin D04 V1W8,Ireland.3Shenzhen Branch,Guangdong Laboratory for Lingnan Modern Agriculture,Genome Analysis Laboratory of the Ministry of Agriculture,Agricultural Genomics Institute at Shenzhen,Chinese Academy of Agricultural Sciences,Shenzhen 518120,China.4Department of Animal Nutrition and Feed Science,College of Animal Science and Technology,China Agricultural University,Beijing 100193,China.5UCD Conway Institute of Biomolecular and Biomedical Research,University College Dublin,Dublin D04 V1W8,Ireland.

Received:4 March 2020 Accepted:7 July 2020