Bovine mastitis is a multi-etiological disease in dairy cattle, characterised by inflammation of the mammary gland. It affects animal health, milk production and quality and causes substantial economic losses to dairy farmers (1). Additionally, treating bovine mastitis significantly contributes to the bigger problem of antimicrobial resistance (AMR), particularly through the potential use and misuse of antibiotics in livestock production (2). Various culture-dependent and culture-independent high-throughput research studies have reported the association of hundreds of microbial species, including Streptococcus agalactiae, Streptococcus uberis, Staphylococcus aureus, Streptococcus pyogenes, Streptococcus dysgalactiae, Trueperella pyogenes, Escherichia coli, Klebsiella pneumoniae, Klebsiella oxytoca, Enterobacter aerogenes, and Pasteurella spp. 16S rRNA and whole metagenome shotgun sequencing (WMGS) have significantly contributed to the discovery of microbial species associated with bovine mastitis (3–6). A 16S sequencing-based study reported the top 10 genera of causative pathogens in the milk of mastitic quarters from 65 cows, and also reported less-known bacteria, such as Sneathia sanguinegens and Listeria innocua, which are difficult to identify using culture-based diagnostics (7). Another 16S study highlighted significant variation in taxonomic profiles among different udder quarters (healthy, mastitis-affected, and quarters with undetermined status) at the genus level (8). Similarly, another study investigated the impact of bacteria causing subclinical mastitis on the structure of the cow’s milk microbiome, reporting that Firmicutes and Proteobacteria were predominantly present in subclinical mastitis (9). Burakova et al. (10) investigated the relationship between milk microbiome composition and bovine mastitis before and after antibiotic treatment using 16S rRNA sequencing. This study linked the increased abundance of the genera Hymenobacter and Lachnospiraceae NK4A136 group with the development of subclinical and clinical mastitis. In contrast, a reduced abundance of Ralstonia, Lachnospiraceae NK3A20 group, Acetitomaculum, Massilia, and Atopostipes in mastitic milk was linked to their potential role in maintaining udder health. In the 16S sequencing, specific regions of rDNA are targeted to characterise the community composition. However, it has several major limitations, including low taxonomic resolution (primarily at the genus level), PCR bias, limited functional insight, and a lack of viability assessment of the bacterial community (11–13). To overcome these limitations, the WMGS approach was employed to gain a deeper understanding of the functional aspects of microbial species associated with bovine mastitis. A milk WGS metagenome study on healthy and clinical mastitis (CM) cows reported 363 unique bacterial species and strains in CM samples (4). Another metagenomic analysis was performed on the subclinical mastitis milk samples of Kankrej, Gir (Bos indicus), and crossbred (Bos taurus × B. indicus) animals. A total of 56 different species, with varying abundances, were detected in the subclinical mastitis milk samples (6). Recently, a longitudinal study on the udder microbiome of Norwegian Red dairy cows was conducted using a shotgun metagenomic approach to gain insights into pathogen-driven microbial adaptation and succession. This study revealed that samples with low somatic cell counts were enriched with beneficial genera, such as Corynebacterium, Bradyrhizobium, and Lactococcus, while Staphylococcus predominated in high somatic cell count milk samples (14). In WMGS, the whole stretch of DNA is sequenced to identify the presence of genes and associated bacteria. However, it fails to distinguish between transcriptionally active and dead bacteria, and has low real-time functional resolution at the species and gene levels. These limitations can be addressed through a metatranscriptomics approach.

The accurate identification of molecular weapons used by disease-associated pathogens against their hosts is critical for gaining insight into disease initiation and progression (15). RNA-Seq datasets of diseases contain extensive functional information on disease-related pathogens that are expressed simultaneously (16). However, host-centric studies often overlook the functional importance of microbial communities, such as in bovine mastitis. The proportion of opportunistic and commensal microbial populations in the udder tissue constantly shifts as mastitis progresses, which can be utilised to understand in vivo host-pathogen interactions. Therefore, publicly available bovine mastitis dual RNA-Seq studies have been analysed to identify transcriptionally active bacteria associated with bovine mastitis, elucidating their functional role in pathogenesis. We extracted the non-bovine, non-ribosomal sequenced reads and constructed de novo metatranscriptome assemblies (17, 18). The assembled transcripts were further annotated to reveal the taxonomy and functional profile of transcriptionally active pathogens associated with bovine mastitis.

Host-centric RNA-Seq datasets are a valuable resource for exploring in vivo gene expression profiles of pathogens to develop new diagnostic and therapeutic solutions. In this study, non-ribosomal and unmapped reads from the host genome were extracted from all the selected RNA-Seq datasets and used for metatranscriptome assembly and analysis of pathogens. This study aims to identify transcriptionally active pathogens and their associated molecular signatures for bovine mastitis disease. The schematic workflow of the performed study is detailed in Figure 1.

This study primarily focused on identifying transcriptionally active bacteria from publicly available bovine mastitis dual RNA-Seq studies, as these datasets contain extensive functional information on disease-related microbial species within the same biological samples. The microbial content expressed during the host-pathogen interaction, provides in vivo gene expression profile of pathogens, which facilitates the discovery of bacterial virulence factors to develop new diagnostics and therapeutics. Similar attempts have been made to understand the role of microbial communities in human cancer diseases using RNA-Seq data (34, 35).

In this study, rigorous quality control criteria were used to exclude chimeric, low-quality, partial and redundant transcript sequences. A minimum protein length of ≥60 was used to include small proteins from microbial communities that are involved in regulating larger proteins, cell signaling, antibiotic production and toxin regulation, membrane-associated functions, metabolism, and stress response, as well as metal ion homeostasis (36). Moreover, the inclusion of 60 amino acids provided higher confidence in functional gene-to-species assignment compared to shorter lengths and relatively less partial sequences. The average read counts of greater than or equal to 3 for genes across the samples were considered to remove chimeric and rare protein sequences, due to the hypersensitivity of high-throughput sequencing technology. The average read count is a better metric for reducing disparity between small and large numbers of studies. Various isoforms of microbial gene sequences were generated due to the stitching of short k-mers during transcriptome assembly. Therefore, we clustered protein isoforms at 90% sequence identity to reduce functional redundancy and false specificity for cross-species orthology inference. Stringent BLAST similarity search criteria (sequence identity ≥75% and query coverage ≥60%) were used to ensure the high-quality taxonomical and functional assignment (Table 2). This approach is more reliable and accurate than the seed sequence mapping approach, which lacks robustness for assigning gene species links in microbial community analysis.

Even after excluding non-host reads before assembly, a large number of the proteins belonged to the bovine species, which may be due to the presence of conserved domains and motif sequences across the kingdom. In filtered protein sequences, the proportion of bovine sequences varied from 3% to 71%, compared to the low bacterial protein proportion (0.5%–10%) across all projects, except for study PRJNA551141. It has 75% bacterial protein and 5% bovine protein. The enhanced representation of bacterial proteins may be due to the fact that rRNA was removed during sample preparation, as rRNA constitutes 70%–90% of the total RNA in these samples. None of the other projects had used ribo-depletion during their sample preparation.

The microbial diversity, species richness and evenness metric support the functional diversity among the selected RNA-Seq studies. The scatter plot of species richness and Shannon diversity index shows a positive correlation between species richness and Shannon diversity, as indicated by the upward-sloping linear regression line (Figure 3D). This suggests that sample PRJEB43443 demonstrated a relatively high Shannon index with moderate species richness, potentially indicating a more evenly distributed community structure, as supported by its larger point size (Simpson index). In contrast, sample PRJNA551141, while showing the highest richness, exhibited only moderate Shannon diversity, suggesting dominance by a few taxa (Stenotrophomonas and Pseudomonas species) and reduced evenness. Higher or moderate evenness often suggests stable communities, while low evenness may indicate vulnerability to disturbance (PRJNA591729 and PRJNA627642). Low Shannon diversity and high species richness indicate dysbiosis or infection dominance (PRJNA551141). Moderate Shannon diversity and moderate species richness may reflect baseline microbial communities (PRJNA544129), whereas moderate diversity with fewer taxa suggests a balanced microbial environment with minimal colonisation (PRJNA778892), possibly post-antibiotic or in the early stage of infection.

In KO analysis, enriched pathways include the two-component system, amino acid/cofactor biosynthesis, flagellar assembly and secretion systems, and biofilm formation, reflecting bacterial adaptation to host defences, nutrient limitation, and persistence mechanisms. Enrichment of amino acid metabolism further highlights adaptive responses of the pathogenic community to oxidative stress and immune modulation (37–41). The Project-wise KO analysis is provided in Supplementary Figures S2A–E. Various highly abundant sequences were either uncharacterised or hypothetical despite their known genomic sequences (Table 4). Therefore, in silico sequence analysis of hypothetical and uncharacterised proteins was performed to explore their functional information. On average, 62% of the total putative proteins were annotated, which contain domains of catalytic enzymes, transposable elements, retrotransposable elements, reverse transcriptases, prokaryotic membrane lipoproteins, endonucleases, and various other molecular signatures. A summary of hypothetical protein annotation and detailed descriptions of the InterProScan annotation of all hypothetical proteins are provided in Supplementary Table S4 and Supplementary File 4 (ST1).

In our analysis, various abundant hypothetical and uncharacterised proteins contain intrinsically disordered regions, suggesting that these proteins may have roles in host–microbe interactions in vivo conditions, such as bacterial effector protein translocation, evasion of the host immune system, and mimicing host protein functions (42). The presence of intrinsically disordered regions in secreted bacterial effectors is well-reported (43). The biofilm-forming curli protein of E. coli has disordered regions that facilitate its aggregation and contribute to the formation of the extracellular matrix of the biofilm (33). Various intrinsically disordered regions containing hypothetical or uncharacterised proteins exhibit a virulent nature in our analysis. The in silico annotated sequence contains reverse transcriptases, transposable elements and transposon-related domains. Reverse transcriptase converts RNA into double-stranded cDNA, and it is a major component of diversity-generating retroelements (DGRs) in bacteria. Reverse transcriptases in DGRs can create hypervariable regions in target genes, fuelling the faster evolution of bacterial genomes to ensure bacterial survival in harsh environments (44). Transposable elements can also alter the genome by insertion and excision, promoting bacterial evolution. These jumping genes are known to be part of bacterial virulent elements (45). These proteins may have roles in antibiotic resistance and pathogenicity, which need to be explored further.

Based on the number of virulent proteins (≥10), peptidases, secretory proteins, and antibiotic resistance (AMR) genes, a list of 25 transcriptionally active bacterial genera was compiled to highlight pathogenic load associated with bovine mastitis. Pseudomonas species dominated in the constructed metatranscriptome with the highest number of species (205) and expressed genes (3898), including the highest number of virulence (1952), secretory (198), peptidases (231) and antimicrobial resistance genes (44), highlighting a significant role of Pseudomonas species in bovine mastitis pathogenesis. Moreover, P. putida, and P. aeruginosa are emerged as multidrug-resistant pathogens linked with clinical and subclinical bovine mastitis worldwide (46–49). Stenotrophomonas species, a lesser-known pathogen of bovine mastitis, have been shown to express substantial virulence genes (1226) and secretory proteins (170), indicating their active involvement in pathogenesis and potential biofilm formation or host interaction, making them an emerging global opportunistic pathogen (50). Recent research also demonstrated that enterogenic S. maltophilia can migrate from the gut to the mammary gland via the gut-mammary axis to induce mastitis by activating the calcium-ROS-AMPK-mTOR-autophagy pathway (51). Mycoplasma species, despite having a lower virulence gene count (175) out of a total of 243 identified genes, indicate a specialisation in pathogenesis rather than metabolic versatility. Staphylococcus, a known pathogen of mastitis, displayed a lower number of species but a notable number of virulence factors (132). The expression of hypothetical or uncharacterised sequences (112) suggests the presence of novel molecular factors from these species in bovine mastitis pathogenic consortia.

Among all identified virulent genes, 633 were uncharacterized, and 379 were hypothetical by Uniprot definition. The protein sequence (Sequence id: TRINITY_DN43764_c0_g1_i1.p4; Uniprot id of hit: UniRef100_UPI000DEBC0B7) from Pseudomonas aeruginosa is one of the identified virulent hypothetical proteins that also possesses an AMR gene feature (Drug class: disinfecting agents and antiseptics; mode of mechanism: Antibiotic efflux; ARO ID:3003680; Gene name: TriB). In our study, we identified 84 antimicrobial resistance (AMR) genes, some are reported for bovine mastitis, including tetG, tetC, blaZ, sul1, adeB, acrA, APH(3′)-IIb, mexA, mexB, and mexF. These genes are involved in well-known resistance mechanisms, including efflux pumps, antibiotic-modifying enzymes, and target protection proteins. Other genes like adeB, adeF, APH(3′)-IIb, eptA, kpnH, ompA, oqxA, and oqxB were also reported for bovine mastitis. The presence of these AMR genes in this dataset suggests a broader distribution of resistance determinants among mastitis-associated microbiota (52–54). The association of other identified AMR genes with bovine mastitis has not been found in the literature. In the absence of literature support, it is challenging to speculate on the role of these bacteria in pathogenic consortia mastitis, which are predominantly composed of Pseudomonas species. However, the expression of such a large number of AMR genes from Pseudomonas species suggests a potential role in providing protection against antibiotics, which warrants further in-depth investigation. Most Pseudomonas aeruginosa strains possess a Type III Secretion System (T3SS), which plays a crucial role in pathogenesis and has been linked to elevated somatic cell counts in the milk of cows with mastitis (55). Key T3SS-related genes identified in bovine mastitis isolates include exlA, pcrD, and hopAJ2. Additionally, a range of other virulence and secretion system genes—such as clpV1, flgC, flgG, flgK, flgL, fliC, mucA, mucB, algD, algF, algA, pilV, pilY2, pvdF, pvdS, rhlR, rhlI, rhlA, rhlB, rhlC, tssB1, tssF1, tssG1, vgrG1, as well as components of the xcp (PA3095–PA3105) and hxc (PA0677–PA0687) secretion system are regulated by phosphate availability (56).

The negative correlation of Klebsiella with other genera of pathogenic consortia (Bacillus/ Balutia/ Salmonella/ Leptospira/ Porphyromonas/ Ruminococcaceae/ Lactobacillus) in the non-directional microbial co-occurrence network indicates that when Klebsiella abundance increases, Bacillus/ Balutia/ Salmonella/ Leptospira/ Porphyromonas/ Ruminococcaceae/ Lactobacillus tends to decrease, and vice versa. This could be due to competition, niche exclusion, or antagonistic interactions. These species may not co-exist well in the same environment. Some species, such as Bacillus/Lactobacillus, produce bacteriocins that inhibit the growth of Klebsiella (57). Host factors (such as temperature, oxygen level, or pH) may favour the growth of Klebsiella species after environmental E. coli infection compared to other pathogens, due to similar growth conditions to those of other Enterobacteriaceae. The chance of Klebsiella infection after exposure to environmental mastitis conditions is moderate to high in poor bedding conditions (58). Blautia and Lactobacillus, a genus with probiotic characteristics, may decrease during disease, whereas the genus Bacillus contains both pathogenic and probiotic species; perhaps the abundance of probiotic strains is more prevalent in healthy conditions. It is not easy to speculate and interpret until further multi-bacterial culture-based research investigation. Different Bacillus strains have been reported to exhibit antimicrobial, antioxidant, and immune-modulatory activities in the host. A study by Pinchuk et al. (32) has reported the B. subtilis strains’ anti H. pylori activity, which was attributed to the secretion of aninocoumacin A antibiotic. The antagonistic activity of aninocoumacin A was also documented against enteric E. faecium and Shigella flexneri. An interesting communication by Ripert et al. (59) revealed that the probiotic B. clausii O/C strain protected Vero and Caco-2 cells from the cytotoxic effects of Clostridium difficili and B. cereus toxins.

Pseudomonas, Stenotrophomonas, Comamonas, and Sphingomonas are among the most transcriptionally active pathogens. Pseudomonas and Sphingomonas genera are often reported in association with bovine mastitis in various studies (60, 61). The Stenotrophomonas genus is a lesser-known opportunistic pathogen associated with bovine mastitis. However, various Stenotrophomonas isolates were extracted and reported from bovine milk samples for mastitis cases (62). The Prevalence and survival of Stenotrophomonas species in milk and dairy products have also been reported (63, 64). A study has explored the prevalence of genetic relatedness, antimicrobial resistance, biofilm formation, biofilm genes associated with virulence and integron genes among isolates of S. maltophilia recovered from bovine milk with subclinical mastitis (65). Comamonas is another less-known genus associated with bovine mastitis. However, it was isolated from bulk tank milk (66). Therefore, our study emphasises the inclusion of these prominent pathogens in the list of primary bovine mastitis-causing pathogens.

This is the first in-depth metatranscriptomics study on bovine mastitis, delineating transcriptionally active pathogenic bacterial taxa, their associated molecular signatures, and uncovering a complex pathogenic consortium that goes beyond the well-known causative pathogens. Pseudomonas (P. putida, P. aeruginosa, P. fluorescens, P. monteilii, P. mosselii, and P. plecoglossicida), Stenotrophomonas (S. maltophilia, and S. sp. RIT309, SKA14), Mycoplasma (M. bovis, M. hyorhinis, M. agalactiae, and M. ovipneumoniae), Staphylococcus (S. aureus, S. succinus, and S. sciuri), Escherichia (E. coli CFT073, and NA114), Comamonas (C. aquatica, C. testosteroni) and Sphingomonas (S. sp. IBVSS2, LK11, FARSPH, S. sanguinis, and S. taxi) species were predominant among transcriptionally active pathogens, and should be considered as significant contributors to bovine mastitis pathogenesis alongside established pathogens.

Identified 213 hypothetical and uncharacterized proteins from Staphylococcus, Mycoplasma, and Escherichia, providing a list of unexplored virulence factors that could serve as targets for multi-epitope vaccine development, while the expression of 74 AMR genes, particularly from Pseudomonas and Escherichia species, suggests these pathogens may confer antibiotic protection to the entire pathogenic consortium, necessitating the need to update current treatment strategies. The abundance of uncharacterized proteins presents opportunities for developing novel diagnostics and therapeutics. The negative correlation of beneficial bacteria (Blautia, Bacillus, and Lactobacillus) with pathogenic species indicates that microbiome modulation could be a viable prevention strategy for subclinical mastitis. A balanced microbiome with these antagonistic bacteria could help to reduce pathogen loads and prevent subclinical infections from becoming clinical. However, these findings need to be validated thoughtfully to develop microbiome-assisted prevention strategies for subclinical bovine mastitis. Our study highlights that transcriptionally active bacteria and their expressed genes (AMR, secretory, peptidase and virulence genes) can aid in identifying novel therapeutic targets, early diagnostic biomarkers and expressed candidate genes for developing a multi-epitope vaccine that targets multiple pathogens. Our findings also suggest the potential of microbiome modulation strategies to develop microbiome-assisted prevention approaches for better management of subclinical bovine mastitis.