EEB5300_SOMK

Differential Gene Expression Analysis Using Milk Somatic Cell Data from Dairy Cows with Subclinical Mastitis

Overview

Mastitis, the inflammation of the mammary gland, remains one of the most economically devastating diseases in dairy cattle worldwide [2]. Subclinical mastitis (SCM), unlike its clinical counterpart, occurs without overt symptoms but significantly elevates somatic cell count (SCC) in milk, reflecting an immune response often unnoticed until substantial productivity losses are incurred [5]. It compromises milk quality, reduces yield, and silently impacts animal welfare and herd longevity [4]. Because of its asymptomatic nature, SCM often escapes early detection, making genomic approaches invaluable for uncovering molecular markers associated with its onset and progression [5].

Motivation

This project investigates what genes may be downregulated or upregulated due active subclinical mastitis. Our motivation stems from the urgent need to integrate functional genomics into dairy breeding programs. While traditional phenotypic selection based on SCC and milk yield remains standard, it is limited by low heritability and delayed expression of traits. Molecular markers identified through RNA-Seq can serve as early indicators of susceptibility, improving selection accuracy and accelerating genetic progress [6].

Moreover, understanding host gene expression in naturally occurring infections, rather than induced models, enhances the translational value of our findings [7]. This project not only contributes to mastitis biology but also showcases the power of integrating public RNA-Seq datasets, advanced bioinformatics tools, and biological interpretation to reveal hidden patterns of disease resistance.

Related Work/Citations

Our project is a subset of previous work conducted by a group from the Sherbrooke Research and Development Centre in Quebec, Canada in collaboration with the Department of Animal Science at Laval University in Quebec City, Quebec. In this publication, entitled "Gene co-expression in response to Staphylococcus aureus infection reveals networks of genes with specific functions during bovine subclinical mastitis", differential gene expression analysis between healthy cows and mastitis affected cows revealed 4,077 differentially expressed genes (DEG). 1,616 genes were upregulated and 2,461 genes were downregulated [1]. The data used in this project was collected as part of an effort to improve understanding of the genetic changes in the mammary gland and the molecular mechanisms underlying the staphylococcus aureus mastitis, as well as reveal candidate discriminant genes with potential regulatory roles in response to S. aureus infection.

(1) Wang M et al., "Gene co-expression in response to Staphylococcus aureus infection reveals networks of genes with specific functions during bovine subclinical mastitis", Journal of Dairy Science, 2023;106(8):5517-5536

(2) Bradley, A. J. (2002). Bovine mastitis: An evolving disease. The Veterinary Journal, 164(2), 116–128. https://doi.org/10.1053/tvjl.2002.0724.

(3) The Gene Ontology Consortium. The Gene Ontology resource: enriching a GOld mine. Nucleic Acids Research, 49(D1):D325–D334, 2021. https://doi.org/10.1093/nar/gkaa1113

(4) Seegers, H., Fourichon, C., & Beaudeau, F. (2003). Production effects related to mastitis and mastitis economics in dairy cattle herds. Veterinary Research, 34(5), 475–491. https://doi.org/10.1051/vetres:2003027

(5) Halasa, T., Huijps, K., Østerås, O., & Hogeveen, H. (2007). Economic effects of bovine mastitis and mastitis management: A review. Veterinary Quarterly, 29(1), 18 31. https://doi.org/10.1080/01652176.2007.9695224

(6) Kommadath, A., et al. (2014). Gene co-expression network analysis of the bovine mammary transcriptome reveals novel insights into the transcriptional regulation of host responses to E. coli infection. BMC Genomics. https://doi.org/ 10.1186/1471-2164-15-1005

(7) Loor, J. J., Everts, R. E., Bionaz, M., Dann, H. M., Morin, D. E., Oliveira, R., ... & Lewin, H. A. (2011). Nutrition-induced ketosis alters metabolic and signaling gene networks in liver of periparturient dairy cows. Physiological Genomics, 43(7), 388 398. https://doi.org/10.1152/physiolgenomics.00195.2010

(8) Lippolis, J. D., & Reinhardt, T. A. (2005). Proteomic survey of bovine neutrophils. Veterinary Immunology and Immunopathology, 108(1–2), 53–62. https://doi.org/ 10.1016/j.vetimm.2005.06.011

(9) Andrews, Simon. “FastQC- Per Sequence GC Content.” Babraham Bioinformatics , 3 Jan. 2023, www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/5%20Per%20Sequence%20GC%20Content.html.

Software Packages

Shifu Chen. 2023. Ultrafast one-pass FASTQ data preprocessing, quality control, and deduplication using fastp. iMeta 2: e107. https://doi.org/10.1002/imt2.107

Andrews, S. (2010). FastQC: a quality control tool for high throughput sequence data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc

Ewels, P., Magnusson, M., Lundin, S., & Käller, M. (2016). MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 32(19), 3047–3048. https://doi.org/10.1093/bioinformatics/btw354

Kim, D., Paggi, J. M., Park, C., Bennett, C., & Salzberg, S. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nature Biotechnology, 37(8), 907–915. https://doi.org/10.1038/s41587-019-0201-4

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., ... & Durbin, R. (2009). The Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078–2079. https://doi.org/10.1093/bioinformatics/btp352

Okonechnikov, K., Conesa, A., & García-Alcalde, F. (2016). Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics, 32(2), 292–294. https://doi.org/10.1093/bioinformatics/btv566

Anders, S., Pyl, P. T., & Huber, W. (2015). HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics, 31(2), 166–169. https://doi.org/10.1093/bioinformatics/btu638

Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550. https://doi.org/10.1186/s13059-014-0550-8

Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2. https://doi.org/10.1093/biostatistics/kxw041

Research Goals

We initially aimed to identify mastitis-related genes that were differentially expressed between healthy and affected cattle using historic RNA-seq data from milk somatic cells. Milk, a non-invasive sampling medium, reflects the host immune response directly at the infection site. Somatic cells—primarily neutrophils, macrophages, and epithelial cells—serve as effective biosensors for intramammary infection [8]. By analyzing transcriptomic changes in these cells, we sought to characterize the genetic landscape associated with early mastitis.

This study also served as a subset replication of prior work, with the goal of developing proficiency in differential gene expression (DGE) analysis. As the project progressed, we refined our focus to identifying the most significantly up- and downregulated genes and interpreting their biological roles. To support this, we performed Gene Ontology (GO) enrichment analysis, applying methods from the UCONN CBC RNA-Seq Tutorial to our dataset.

Data

For this project, we used historical data that was prepared and sequenced at the Animal Genomics and Epigenomics Laboratory of the Sherbrooke Research and Development Centre, Agriculture and Agri-Food Canada in 2023. RNA was isolated from the milk somatic cells using the RNeasy Mini Kit according to the manufacturer's protocol. The quantity and integrity of total RNA were checked with Agilent BioAnalyser 2100 and LabChip GXII instruments, respectively. Then, the NEBNext RNA First Strand Synthesis and NEBNext Ultra Directional RNA Second Strand Synthesis Modules from NEBNext Ultra RNA Library Prep Kit were used for cDNA synthesis. The original data contains 17 RNA-seq libraries from Staphylococcus aurus affected cows and 9 RNA-seq libraries from healthy cows. Our subsetted data contains 3 libraries from each disease-status (S.aureus and healthy) and they were chosen because they contained the highest number of paired-end reads (>20M) with at least 100bp in length in each library. This raw data is organized as fasta files and is stored on UCONN's Xanadu cluster. The bos taurus genome is used as the reference genome and was downloaded from NCBI's database, along with the RNA-seq libraries.

Workflow/Analysis

01. RNA Library Selection

The top 3 libraries with the highest coverage for both healthy and mastitis affected cows were selected. It is important to note that for RNA-seq studies, 3 libraries is the minimum amount of coverage necessary to get decent results. The accession list was made manually:

Accessions
SRR24593481
SRR24593474
SRR24593469
SRR24593488
SRR24593486
SRR24593484

The top three accession numbers are samples from healthy cows, and the bottom three accessions are samples from mastitis-affected cows.

02. Downloading the genome, transcripts and annotation files from NCBI

The SRA Toolkit was used to download raw sequencing reads (FASTQ files) from the NCBI Sequence Read Archive (SRA). The parallel software was used to run mutliple jobs at once to ensure efficient computational usage.

module load parallel/20180122
module load sratoolkit/3.0.1

The complete script is called 01_fasterq_dump.sh and can be found here. The script get_genome.sh downloads the Bos taurus reference genome from NCBI using wget and can be found here. It also uses samtools to generate FAI index files for both the genomic and coding sequence (CDS) FASTA files. The indexes enable efficient, random access to specific regions of the genome or transcript sequences during downstream analysis.

03. Metadata Construction

The metadata table was made to include important information, such as the accession number, associated bioproject, layout, biosample number, disease-status, organism, and the RNA-seq library name:

Run	Assay Type	LibraryName	Bioproject	Layout	Biosample	Disease_Status	Organism
SRR24593481	RNA-Seq	EIAMK3196CT_negSA_RNAseq	PRJNA878880	PAIRED	SAMN35038917	Healthy	Bos Taurus
SRR24593474	RNA-Seq	EIAPI920CT_negSA_RNAseq	PRJNA878880	PAIRED	SAMN35038919	Healthy	Bos Taurus
SRR24593469	RNA-Seq	EIASR52CT_negSA_RNAseq	PRJNA878880	PAIRED	SAMN35038924	Healthy	Bos Taurus
SRR24593488	RNA-Seq	EIAME32q3SA_posSA_RNAseq	PRJNA878880	PAIRED	SAMN35038928	Subclinical Mastitis	Bos Taurus
SRR24593486	RNA-Seq	EIAME9q3SA_posSA_RNAseq	PRJNA878880	PAIRED	SAMN35038930	Subclinical Mastitis	Bos Taurus
SRR24593484	RNA-Seq	EIAMK2700q4SA_posSA_RNAseq	PRJNA878880	PAIRED	SAMN35038932	Subclinical Mastitis	Bos Taurus

04. Quality Control Using Fastp, FastQC, and MultiQC

We chose to use Fastp, FastQC, and MultiQC to ensure comprehensive and efficient quality control of our paired-end RNA-seq reads. Fastp was selected as a versatile preprocessing tool capable of performing adapter trimming, quality filtering, and base correction in a single step with minimal configuration. Its speed and multithreading support made it ideal for handling large datasets, and it automatically generated detailed quality reports in both HTML and JSON formats. The complete script is called fastp_trim.sh and can be found here.

module load fastp/0.23.2
module load parallel/20180122

Following preprocessing, we used FastQC to independently assess the quality of each read file, providing insights into per-base quality scores, sequence duplication, and overrepresented sequences. To streamline the interpretation of results across multiple samples, we applied MultiQC, which aggregates the FastQC and Fastp reports into a single, easy-to-navigate summary. This combination of tools allowed us to validate the integrity and cleanliness of our data before proceeding to downstream analyses. The complete scripts are called fastqc.sh and multiqc.sh and can be found here and here, respectively.

module load fastqc/0.11.7
module load parallel/20180122

module load MultiQC/1.9
multiqc -f -o fastp_multiqc fastp_reports
multiqc -f -o fastqc_multiqc fastqc_reports

The multiqc output gave us an html file that nicely demonstrated the quality of each RNA-seq library. Below are some visuals from a compilation of quality status checks.

Results from quality control indicate high overall data quality, characterized by minimal adapter contamination, negligible per-base N content, and consistently high per-base sequence quality with PHRED scores averaging around 36 across all trimmed reads. These metrics confirm the RNA-seq libraries are suitable for reliable alignment and downstream analysis. One metric of concern was the many poor quality per sequence GC content across the majority of the libraries, as seen below:

Typically, the average per-sequence GC content follows a normal distribution because each sequencing read represents a random sample from the transcriptome. In a homogeneous, uncontaminated sample, the GC content across reads tends to cluster around a central mean and form a symmetrical, bell-shaped curve characteristic of a normal distribution. Poor per sequence GC content in RNA-seq studies are generally less concerning than in DNA-seq studies, especially if they are consistant across samples [9]. Some explanations for the poor GC content could be a contaminated sample, as half of the RNA-seq samples came from cows affected with Staphylococcus aureus, which could have caused the difference we are seeing.

Most importantly, key quality metrics — including basic statistics, per-base sequence quality, per-sequence quality scores, per-base N content, and adapter content — all passed quality control.

05. Indexing and Alignment

Once trimming and quality checks were completed, the next step was to index and align the RNA-Seq reads to the appropriate reference genome. Indexing RNA-seq reads using hisat2-build involves creating a genome index from a reference genome sequence, which allows HISAT2 to efficiently align RNA-seq reads to the genome during the mapping process. It is specifically designed for mapping reads from spliced transcripts, which is essential in eukaryotic systems like bovine mammary tissue where exon-exon junctions are common.

The cleaned, paired-end reads were aligned to a pre-indexed reference genome using HiSat2 with strand-specific parameters suitable for the dataset. The output files were initially in SAM format and were converted to sorted BAM files using Samtools, which helped prepare the data for further analysis. An index file was also generated to make the files compatible with other downstream tools.

After this step, alignment success was monitored by checking mapping percentages and ensuring that most reads aligned uniquely to the genome. The alignment procedure was tailored to preserve strand information and optimize accuracy. This step was critical to ensure that gene counts calculated in the next stage would reflect true biological expression levels.

The full script name is align.sh and is found here.

06. Post-Alignment Quality Control

Quality control post-alignment was done using samtools and qualimap.

Mapping statistics were generated with samtools, which provided essential insights into the quality and efficiency of our sequence alignments. These statistics typically include the total number of reads, the percentage of reads mapped to the reference genome, the percentage of proprely paired reads, and the distribution of reads across different genomic regions. Below is a chart demonstrating the percentage of reads mapped to the reference genome:

This ensured us that the majority of the reads in the libraries we chose were able to be mapped appropriately to the reference genome, and it is important because it indicates the samples are of good quality.

Post-alignment analysis via Qualimap showed a comparable proportion of exonic and intronic reads. This pattern may reflect partial RNA degradation or residual genomic DNA contamination during sample preparation. Below is the chart demonstrating this: While the majority of transcripts of interest are still detectable, this distribution suggests caution in interpreting expression levels, particularly for low-abundance genes.

Because milk somatic cells can be heterogeneous and partially degraded in the udder, some variability in RNA quality was expected. While not all libraries reached optimal exonic percentages, the data remained usable, and key immune-related expression patterns were robust across samples.

07. Generating Counts of Fragments Mapping to Genes

We used the program htseq to count how many RNA fragments map to each annotated gene in the genome.

module load htseq/0.13.5
module load parallel/20180122

The full name of this script is "htseq_count.sh" and can be found here

08. Pairwise Differential Expression with Counts in R using DESeq2

To identify the differentially expressed (DE) genes, we used the R package DESeq2, a part of the Bioconductor project. For differential expression and visualisation, the following R packages must be installed:

DESeq2
pheatmap
tidyverse
ggrepel
ashr
ggplot2

The first thing we did for this step was make sure we had our counts files in a local folder. From there, we created a DeSeq2 dataset using the HTSeq count files we generated in the previous step. The DeSeq2 performed the normalisation and the model fitting. We then extracted raw differential expression results and shrunk the log2 fold changes in order to stabilise the effect size estimates.

Next, we generated a scatter plot to compare the raw vs shrunken log2 fold changes. Each point is colored by adjusted p-value significance. An identity line was added (grey diagonal) for reference.

This plot reveals that a substantial number of genes were differentially expressed, with a clear distribution of upregulated and downregulated genes relative to mean expression, emphasizing the biological impact of mastitis.

After this, a PCA plot was generated to explore any patterns present across the RNA-seq libraries across conditions.

This PCA plot shows a distinct separation between control and mastitis-affected samples, indicating strong transcriptomic differences and consistent sample grouping.

Finally, we generated a heatmap demonstrating the top 20 most differentially expressed genes.

This heatmap further validates the clustering of gene expression patterns by health status.

The full R code for this part of the analysis is located here.

Limitations

RNA-seq analyses, like most high-throughput technologies, present several challenges that must be addressed to obtain reliable and high-quality results. One key factor in RNA-seq studies is the timing of sample collection; gene expression patterns are only reflective of the moment the sample was taken, making the timing of extraction critical for accurate interpretation. RNA is also highly susceptible to degradation, and tools like Qualimap can provide valuable insights into RNA integrity. Researchers must be cautious about RNA quality when interpreting results. Although we used the minimum of three biological replicates per condition, increasing the number of libraries would have likely enhanced statistical power and improved the robustness of our findings. Additionally, it is essential to recognize that differential expression alone does not necessarily imply functional relevance. Therefore, it is crucial to further investigate the biological functions of the differentially expressed genes, ensuring that the results are interpreted in the context of their potential biological significance.

Discussion

Among the upregulated genes, CXCR2, OSM, and THBD stood out for their roles in immune defense. CXCR2 is involved in neutrophil recruitment, OSM is a cytokine that amplifies inflammatory signaling, and THBD contributes to vascular repair and modulation of inflammation. These genes reflect an activated immune state in response to intramammary infection. Conversely, genes such as GPAM, GSTM3, and HIC1 were significantly downregulated. These genes are typically involved in lipid metabolism, oxidative stress protection, and growth regulation, suggesting that during infection, the mammary gland redirects energy from normal physiological functions toward managing immune challenges. The table below is a good reference for more information regarding gene functions:

When compared to a recent transcriptomic study by Wang et al. (2023), which focused specifically on natural Staphylococcus aureus subclinical mastitis, there was strong agreement in the findings. Both studies identified significant upregulation of immune-related genes such as CXCR2, IL1B, NFKB1, and THBD, and downregulation of genes involved in tissue development and metabolic pathways. Additionally, both found functional enrichment in immune signaling pathways and suppression of processes like cell adhesion and lipid synthesis. The overlap in gene behavior and pathway enrichment between the two studies underscores the reliability of the results and highlights conserved molecular responses to mastitis, regardless of specific pathogen strain.

Subclinical mastitis, although asymptomatic, triggers a robust immune response at the molecular level. This study confirms that milk somatic cells serve as a valuable, non-invasive source for detecting these early changes. Transcriptomic profiling of such samples not only improves our understanding of host-pathogen interactions but also opens new avenues for applied research. These findings can inform the development of early diagnostic tools and support genomic selection strategies aimed at enhancing mastitis resistance. By identifying consistent expression patterns and immune markers, this approach contributes to more sustainable dairy farming practices and aligns with broader goals in animal health and production efficiency.

sdo13001/EEB5300_SOMK

About

Resources

Stars

Watchers

Forks

Releases

Languages