Integrative microbiome and transcriptome analyses identify parity-associated alterations in rectal fecal microbiota and peripheral immune transcriptional signatures in dairy cows

Increasing parity in dairy cows is accompanied by cumulative physiological, reproductive, lactational, and management exposures that may coincide with changes in metabolic and immune status. However, the relationships among parity group, distal hindgut microbial communities, peripheral immune-related transcription, and their coordinated variation remain incompletely understood. In this cross-sectional observational study, 40 nonpregnant Holstein cows in mid-lactation without recorded disease were randomly selected from an eligible population within a commercial dairy herd and classified into naturally occurring parity groups: primiparous cows (PrC; parity = 1; n = 20) and multiparous cows (MuC; parities = 2–5; n = 20). Rectal fecal 16S ribosomal RNA (16S rRNA) gene sequencing and peripheral whole-blood RNA sequencing were integrated to characterize microbial and immune-related transcriptional differences between the two groups. Multiparous cows had lower Shannon diversity, Chao1 richness, and observed amplicon sequence variant (ASV) richness than primiparous cows, but a higher Simpson dominance index (D), consistently indicating lower α-diversity and greater community dominance in multiparous cows. A one-factor permutational multivariate analysis of variance (PERMANOVA) based on Bray–Curtis distances indicated that parity-group membership was associated with a modest proportion of the variation in rectal fecal microbial composition (R² = 0.0581, permutation p = 0.004). Using linear discriminant analysis effect size (LEfSe), random forest, and analysis of compositions of microbiomes with bias correction 2 (ANCOM-BC2), 13 candidate genera were identified by at least two analytical approaches. Prevotellaceae_UCG-003 was more abundant in PrC, whereas Turicibacter , Clostridium_sensu_stricto_1 , Acetitomaculum , and Lachnospiraceae_NK3A20_group were more abundant in MuC. Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2 (PICRUSt2)-based inference indicated group differences in predicted functional potential related to protein synthesis, amino acid and carbohydrate metabolism, bacterial motility, chemotaxis, and signal transduction; these results represent predicted rather than directly measured microbial functions. Whole-blood RNA sequencing identified 189 differentially expressed genes in the comparison between MuC and PrC, including 116 genes with higher transcript abundance in MuC and 73 with higher transcript abundance in PrC. Overall KEGG analysis of the complete DEG set identified 38 significantly enriched pathways, many of which were related to immune and inflammatory processes. Direction-stratified KEGG analysis further indicated that broad significant immune-related enrichment was mainly associated with genes showing higher transcript abundance in PrC cows, whereas only one KEGG term remained significantly enriched among the MuC-higher genes after Benjamini–Hochberg correction. Integration of differential expression, weighted gene co-expression network analysis, and protein–protein interaction analysis prioritized matrix metallopeptidase 9 ( MMP9 ), mitogen-activated protein kinase 13 ( MAPK13 ), and WD repeat domain, phosphoinositide interacting 1 ( WIPI1 ) as candidate hub genes, all of which showed higher transcript abundance in MuC. Sparse partial least squares integration further identified a multivariate association structure involving selected candidate fecal genera and immune-related genes. This study provides exploratory multi-omics evidence of differences in rectal fecal microbial composition, predicted microbial functional potential, and peripheral blood immune-related transcription between primiparous and multiparous cows. The coordinated microbiome–transcriptome patterns identify candidate microbial and transcriptional features for future investigation. Because this was a cross-sectional observational study, the findings should be interpreted as associations rather than causal effects of parity. Longitudinal studies and validation in independent herds are required.