Study-Aware Meta-Analysis Identifies Recurrent Pathway Programs and Graded Candidate Stability in Bovine Heat-Stress Transcriptomes

1. Introduction

Heat stress is a recurring constraint on cattle production because the physiological work required to dissipate heat competes with lactation, growth, reproduction, and immune function. Across dairy-cow experiments, heat exposure reduces dry-matter intake, milk yield, and feed efficiency, while consequences extend beyond production to inflammatory control, disease susceptibility, and animal welfare [1,2,3]. Genetic selection and molecular phenotyping may complement environmental cooling and nutritional management, but thermotolerance is polygenic and biologically distributed. Molecular candidates intended for later thermotolerance studies therefore need heat-response evidence that extends beyond one tissue, one exposure protocol, or one processing strategy [4,5].

Bulk RNA sequencing has revealed several parts of the bovine heat-stress response. Mammary studies have described heat-shock, involution, metabolic, and milk-synthesis programs [6,7,8]. Liver studies have implicated mitochondrial energy metabolism, protein folding, inflammation, and nutrient partitioning [9,10]. Skin transcriptomics has linked thermal status to local metabolic, antioxidant, and immune responses [11], while recent blood and multi-omics work has extended the response to systemic immunity, metabolism, and intergenerational exposure [12,13,14]. These studies are complementary, but reported differentially expressed gene lists do not preserve a common effect scale or its uncertainty. Differences in reference annotation, preprocessing, filtering, tissue composition, exposure definition, and statistical power can therefore masquerade as biological disagreement. The unresolved question is not whether cattle respond transcriptionally to heat, but which effects recur on a common scale after study design and small-sample uncertainty are preserved.

Cross-study synthesis introduces a second problem: only a few eligible in vivo studies are available. Standard random-effects meta-analysis estimates an average effect and between-study variance, but normal-approximation intervals can be too narrow when the number of studies is small and study precisions differ. Modified Knapp–Hartung inference uses a $t$ distribution and prevents the adjusted standard error from becoming smaller than its conventional counterpart, yielding more reliable error control in small-$k$ settings [15]. Prediction intervals add a different piece of information by estimating the range of true effects expected in a comparable future context rather than only the uncertainty around the mean effect [16]. These distinctions are especially relevant for transcriptome-wide meta-analysis, where thousands of genes are tested and the number of contributing studies varies by gene.

Here, we asked at which level of biological organization the bovine heat response remains reproducible across heterogeneous in vivo studies. Five public datasets were reprocessed against one bovine reference, while heat-versus-control effects were estimated within each original study design. REML random-effects models and modified Knapp–Hartung inference quantified gene-level evidence. HSTSI then combined inferential support, effect magnitude, study coverage, direction consistency, heterogeneity, and gene-level jackknife stability to prioritize candidates. Weight perturbation, alternative score definitions, and five complete study deletions tested that ranking. Directional enrichment and pathway-level study deletion evaluated whether coordinated biological programs were more recurrent than individual gene ranks. This design separates formal gene-level inference, graded candidate stability, and pathway recurrence, thereby defining both the shared response and its context dependence.

2. Materials and Methods

2.1. Dataset Search, Eligibility, and Contrast Definition

NCBI GEO/DataSets, NCBI SRA, and EBI BioStudies were searched through 11 June 2026. Search terms combined heat stress, thermal stress, heat load, high temperature, or hyperthermia with Bos taurus, cattle, bovine, or dairy cow, together with RNA-seq or transcriptome terms. After title and metadata de-duplication, 23 plausible GEO series underwent eligibility assessment; five were included and 18 were excluded. Primary-analysis eligibility required an in vivo Bos taurus bulk RNA-seq experiment, contemporaneous heat-stress and non-heat comparison groups, at least three biological units per group, and publicly accessible raw reads. In vitro, organoid, single-cell or single-nucleus, developmental-only, technical library-preparation, and longitudinal before–after studies without a concurrent non-heat group were excluded. Complete search strings, candidate records, evidence links, and final decisions are provided in Supplementary Table S1.

Five studies met these criteria: GSE108840, GSE226351, GSE227323, GSE244620, and GSE289946 (Table 1). Contrasts were defined from GEO SOFT records, SRA run metadata, and the corresponding articles before differential-expression analysis [6,9,10,11,12]. GSE108840 contributed 48 post-dry-off mammary libraries from 12 cows and retained time in the model. GSE226351 contributed the thermoneutral and heat-stressed liver groups; its pair-fed group was excluded from the primary contrast. GSE227323 contributed postnatally cooled and heat-stressed calf liver samples. GSE244620 contributed control and heat-stress case skin samples. GSE289946 contributed dam whole-blood samples; offspring and samples without heat-stress information were excluded. The resulting primary dataset contained 107 libraries (Supplementary Tables S2 and S3).

2.2. Read Processing and Transcript Quantification

FASTQ files were checked against the provider-supplied integrity information and evaluated with FastQC 0.12.1 and MultiQC 1.35. Adapter and quality trimming used Trim Galore 2.2.0 with a Phred threshold of 20 and a minimum retained length of 36 nucleotides. All studies were quantified against the Ensembl release 110 Bos taurus ARS-UCD1.2 cDNA reference. Salmon quantification used automatic library-type detection, selective alignment validation, sequence-bias correction, and GC-bias correction [17]. Both single-end and paired-end libraries were handled according to their run metadata.

2.3. Gene-Level Import and Study-Specific Differential Expression

Transcript estimates were summarized to Ensembl gene identifiers with tximport 1.38.2 using type = salmon, ignoreTxVersion = TRUE, and countsFromAbundance = no [18]. This passes estimated counts and effective lengths to DESeq2 so that gene-level offsets retain transcript-length information. Concordance with a parallel transcript-to-gene aggregation was evaluated using library-level Pearson correlations and median absolute log-scale differences (Supplementary Figure S1; Supplementary Tables S4 and S16).

Differential expression was estimated separately for each study with DESeq2 1.50.2 [19]. Genes were retained when counts were at least 10 in at least 20% of the included libraries for that study. GSE108840 used the design ~time + analysis_group; the other four studies used ~analysis_group. Unshrunken maximum-likelihood log₂ fold changes and their Wald standard errors were used as the meta-analysis inputs. In every model, the reported coefficient was heat relative to control, so a positive log₂ fold change denotes higher expression in the heat-stress group. No pooled all-library differential-expression model was used.

2.4. Random-Effects Meta-Analysis and Small-Sample Inference

For gene $g$ in study $i$, the DESeq2 log₂ fold change $y_{gi}$ and its standard error $s_{gi}$ were entered into a univariate random-effects model,

$$y_{gi}={\mu }_g+u_{gi}+{\epsilon }_{gi}, u_{gi}\sim N(0,{\tau }_g^2), {\epsilon }_{gi}\sim N(0,s_{gi}^2),$$

where ${\mu }_g$ is the mean cross-study effect and ${\tau }_g^2$ is the between-study variance. Models were fitted by REML with metafor 5.0-1 [20]. A gene required finite effects and positive standard errors in at least two studies. Conventional REML standard errors and normal-approximation tests were retained for calibration comparisons.

Primary inference used a modified Knapp–Hartung standard error. For each gene, the standard error was

$$s_{g,\mathrm{m}\mathrm{K}\mathrm{H}}=\mathrm{m}\mathrm{a}\mathrm{x}(s_{g,\mathrm{R}\mathrm{E}\mathrm{M}\mathrm{L}},s_{g,\mathrm{K}\mathrm{H}}),$$

and two-sided tests and 95% confidence intervals used a $t$ distribution with $k_g-1$ degrees of freedom [15]. The primary plug-in 95% prediction interval was

$${\widehat{\mu }}_g\pm t_{0.975,k_g-1}\sqrt{s_{g,\mathrm{m}\mathrm{K}\mathrm{H}}^2+{\widehat{\tau }}_g^2}.$$

Because prediction-interval conventions are unstable when $k$ is small, classification was also repeated for genes with $k_g\ge 3$ using the more conservative $t_{0.975,k_g-2}$ critical value. Benjamini–Hochberg correction was applied across all genes with estimable meta-analytic effects [21]. $I^2$, ${\tau }^2$, and the Cochran $Q$ test were recorded as heterogeneity descriptors. Models that did not converge under default Fisher scoring were refitted with a step adjustment of 0.5, a maximum of 10,000 iterations, and a convergence threshold of ${10}^{-8}$. Persistent non-convergence was prespecified as an exclusion criterion. Fisher and signed Stouffer combinations of study-level $P$ values were secondary diagnostics, while inferential conclusions were based on the effect-size model (Supplementary Table S18).

2.5. HSTSI Construction and Sensitivity Analysis

HSTSI was calculated for every gene with an estimable random-effects model. Six components were scaled to $[0,1]$: modified Knapp–Hartung FDR support, $S_q=\mathrm{m}\mathrm{i}\mathrm{n}[-{\mathrm{l}\mathrm{o}\mathrm{g}}_{10}(q_{\mathrm{m}\mathrm{K}\mathrm{H}})/10,1]$; effect magnitude, $S_{\beta }=\mathrm{m}\mathrm{i}\mathrm{n}\left(\right|\widehat{\mu }|,1)$; study coverage, $S_k=\mathrm{m}\mathrm{i}\mathrm{n}(k/5,1)$; direction consistency, $S_d=\mathrm{m}\mathrm{a}\mathrm{x}(n_{\mathrm{u}\mathrm{p}},n_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}})/(n_{\mathrm{u}\mathrm{p}}+n_{\mathrm{d}\mathrm{o}\mathrm{w}\mathrm{n}})$; low heterogeneity, $S_h=1-I^2/100$ after bounding $I^2$ to $[0,100]$; and jackknife support, $S_j=0.7p_{\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{e}\hspace{0.25em}\mathrm{d}\mathrm{i}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}+0.3p_{P_{\mathrm{m}\mathrm{K}\mathrm{H}}<0.10}$. Missing $I^2$ was assigned 0.5, and unavailable jackknife support was assigned zero. The score was

$$\mathrm{H}\mathrm{S}\mathrm{T}\mathrm{S}\mathrm{I}=100(0.25S_q+0.20S_{\beta }+0.20S_k+0.15S_d+0.10S_h+0.10S_j).$$

HSTSI is a prioritization score and was not used as a statistical test.

For genes represented in at least three studies, the gene-level jackknife removed each available study effect in turn and contributed only to $S_j$. Weight sensitivity used 500 perturbations with random seed 20260710; each base weight was multiplied by an independent value from a uniform 0.75–1.25 distribution and then renormalized. We recorded top-200 frequency and the rank distribution for each gene. Definition sensitivity compared the base score with nine alternatives: equal component weights, six leave-one-component-out scores, a hyperbolic-tangent transformation of effect magnitude, and an effect cap of two rather than one. Rank correlation and overlap across the top 100 to top 500 genes were calculated for each alternative (Supplementary Tables S6–S8 and S19).

Separately, five complete leave-one-study-out analyses removed each study from the entire analysis and repeated random-effects fitting, modified Knapp–Hartung inference, multiple-testing correction, the nested gene-level jackknife component, and HSTSI ranking. Stability was summarized by top-200 overlap and Jaccard index, Spearman rank correlations for all genes and the full-data top 200, direction agreement, and candidate-specific rank distributions. Tier A required effects from all five studies, direction consistency at least 0.80, and complete jackknife direction agreement. It also required weight-perturbation and leave-one-study-out top-200 frequencies of at least 0.80. Tier B required effects from at least four studies, direction consistency at least 0.75, jackknife direction agreement at least 0.80, weight frequency at least 0.50, and leave-one-study-out frequency at least 0.60. Remaining genes in the top 200 were assigned Tier C. The tiers describe ranking stability; they do not define statistical significance or low heterogeneity (Supplementary Figure S3; Supplementary Tables S9 and S10).

2.6. Directional Functional Enrichment

Ensembl identifiers were mapped to bovine Entrez identifiers with org.Bt.eg.db 3.22.0. For genes with a finite meta-analytic estimate, the directional ranking statistic was ${\widehat{\mu }}_g/s_{g,\mathrm{m}\mathrm{K}\mathrm{H}}$. When more than one Ensembl identifier mapped to the same Entrez identifier, the entry with the largest absolute statistic was retained; exact ties were resolved by a deterministic infinitesimal Entrez-ID offset. Gene-set enrichment analysis (GSEA) used the complete ranked list, gene-set sizes from 10 to 500, and the fgsea implementation through clusterProfiler 4.18.4 [22,23]. Gene Ontology biological processes and bovine KEGG pathway–gene mappings were tested separately with seeds 20260710 and 20260720, respectively. Over-representation analysis (ORA) was performed for mapped up- and down-prioritized genes in the top 200, using all tested genes with bovine Entrez mappings as the explicit universe. All enrichment families used Benjamini–Hochberg correction. Finally, each study was deleted in turn, the modified Knapp–Hartung ranking statistic was rebuilt, and GO and KEGG GSEA were repeated with fixed deletion-specific seeds. For every full-data FDR-significant term, we recorded direction agreement and FDR recurrence across the five deletions (Supplementary Table S15).

2.7. Independent Numerical Validation

An independent implementation refitted all 16,756 gene-level random-effects models directly from the study-effect table. Estimates, heterogeneity statistics, confidence and prediction intervals, and adjusted probabilities were compared with the primary analysis. HSTSI values, prediction-interval classifications, study-deletion summaries, enrichment directions, and pathway recurrence were also reconstructed from their component tables. Software session information and the validation criteria accompany the derived results (Supplementary Tables S17 and S22).

3. Results

3.1. Five Studies Provided Complementary Contrasts Across Tissues and Life Stages

The primary analysis combined 107 libraries from mammary gland, liver, skin, and whole blood, spanning adult cows, pre-weaning calves, Holstein populations, and Brangus cattle (Figure 1). Standard tximport produced 22,372 gene rows for each study before study-specific expression filtering. Concordance with the parallel aggregation was nearly complete: the minimum library-level Pearson correlation for log_e$(1+\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t})$ was 0.999973, and the median absolute log-scale difference was zero. The numbers of tested genes were 16,401 for GSE108840, 14,421 for GSE226351, 15,525 for GSE227323, 16,760 for GSE244620, and 13,853 for GSE289946. At within-study FDR $<0.05$, the corresponding differential-expression counts were 58, 161, 4, 12,116, and 284. Because these counts also reflect tissue, contrast definition, biological response, and precision, integration used the effect and standard error from each study rather than the count of differentially expressed genes.

The merged effect table contained 18,255 gene identifiers, of which 16,756 had finite effects from at least two studies and were eligible for random-effects synthesis. Of these, 15,396 were represented in at least three studies and could be evaluated by the gene-level jackknife used in HSTSI. Thus, most meta-analyzed genes supported both an average-effect estimate and repeated deletion of their contributing study effects (Supplementary Table S5).

3.2. Small-Sample Calibration Separated Average-Effect Signals from Transcriptome-Wide Evidence

Conventional REML analysis with normal-approximation tests yielded 19 genes at FDR $<0.05$ (Figure 2A). The same mean effects yielded no gene at FDR $<0.05$ or FDR $<0.10$ when evaluated with modified Knapp–Hartung standard errors and $t_{k-1}$ critical values. The inferential contrast was largest for genes supported by few studies or unequal precision (Figure 2B). Thirty-four genes had primary 95% plug-in prediction intervals that excluded zero. With the alternative $t_{k-2}$ critical value, 15 retained an interval excluding zero, and 23 of the 34 primary exclusions had an estimated ${\tau }^2$ of zero. Point estimates were unchanged, but the strength assigned to them depended on small-$k$ calibration and the prediction-interval convention (Supplementary Figure S2).

An independent refit of all 16,756 gene models reproduced every primary estimate, heterogeneity statistic, interval, and adjusted probability within $1.85\times {10}^{-13}$. The inferential contrast therefore reflected the specified small-sample method rather than numerical implementation.

IL1R2 illustrates a stable positive effect. All five study estimates were positive, and the pooled log₂ fold change was 0.820 (modified Knapp–Hartung 95% confidence interval, 0.379 to 1.261; prediction interval, 0.315 to 1.326; $I^2=6.1\mathrm{\%}$; Figure 2C). Its unadjusted modified Knapp–Hartung $P$ value was 0.00668, but its transcriptome-wide FDR was 0.802. GZMK showed a pooled negative effect (log₂ fold change, $-0.937$; 95% confidence and prediction interval, $-1.552$ to $-0.323$; $I^2=0\mathrm{\%}$), with one imprecise positive study estimate and four non-positive estimates (Figure 2D).

values for the same genes; highlighted points are the 19 genes with conventional REML-z FDR $<0.05$

3.3. HSTSI Prioritized Candidates with Complementary Stability Profiles

HSTSI ranked all 16,756 estimable genes by combining six predefined evidence components (Figure 3A). The top 200 contained 67 genes with positive meta-analytic effects and 133 with negative effects. Based on study coverage, direction, jackknife behavior, weight perturbation, and complete leave-one-study-out performance, 43 candidates met Tier A criteria, 31 met Tier B criteria, and 126 met Tier C criteria. Tier assignment did not require low $I^2$, so heterogeneity remained a separate attribute of each candidate.

IL1R2 ranked first (HSTSI 71.04, Tier A). It remained in the top 200 in all 500 weight perturbations and all five leave-one-study-out analyses; its worst deletion rank was 4. SDCBP2 ranked second (HSTSI 70.37, Tier B) and had a positive effect in each of four contributing studies. It remained among the top 200 in every perturbation and deletion analysis (worst rank, 12). GZMK ranked third (HSTSI 69.79, Tier B) and remained directionally negative in every estimable deletion, but its top-200 frequency across study deletions was 0.60 and its worst rank was 628. CD8A ranked sixth (HSTSI 67.23, Tier A), with a pooled log₂ fold change of $-0.599$, $I^2=0\mathrm{\%}$, weight-perturbation top-200 frequency of 1.00, and deletion top-200 frequency of 0.80. The per-study heatmap and sensitivity panel distinguish candidates with similar full-data ranks but different dependence on study composition (Figure 3B,C) (Supplementary Table S20).

The broader top 200 was more sensitive to the HSTSI definition than these leading candidates. Across nine alternative definitions, all-gene Spearman correlation with the base ranking ranged from 0.505 to 1.000, and top-200 overlap ranged from 24 to 200 (Supplementary Figure S4). Removing the FDR-support component left the top 200 unchanged because modified Knapp–Hartung FDR values provided little discrimination. Removing effect magnitude produced the largest change (top-200 overlap, 24), and equal component weights retained 51 genes. IL1R2, SDCBP2, and GZMK nevertheless remained in the top 200 under all ten definitions, while CD8A remained in nine.

Canonical heat-shock genes provided a biological benchmark rather than dominating the stability ranking. HSPA1A had a positive mean effect of 0.581 but ranked 935th with $I^2=90.5\mathrm{\%}$; HSPB1 ranked 1,883rd with $I^2=98.5\mathrm{\%}$, and HSP90AA1 ranked 4,442nd with $I^2=90.2\mathrm{\%}$. HSTSI therefore favored effects that recurred with stable magnitude and direction across the included contexts, rather than genes selected solely for established heat-shock biology (Supplementary Table S21).

3.4. Complete Study Deletion Preserved Direction More Strongly than Rank

Removing one study and recomputing every meta-analytic and HSTSI component changed exact top-200 membership but not the direction of the estimable genes in the full-data top 200 (Figure 4). Top-200 overlap was 117 after removing GSE108840, 106 after removing GSE226351, 83 after removing GSE227323, 68 after removing GSE244620, and 101 after removing GSE289946. The corresponding Jaccard indices were 0.413, 0.361, 0.262, 0.205, and 0.338. Across all ranked genes, Spearman correlations with the full-data ranking ranged from 0.449 to 0.858. Direction agreement for the full-data top 200 was 1.00 in every deletion, whereas median deletion ranks ranged from 169.5 to 461.5.

The difference between directional and rank stability was also visible at candidate level. IL1R2 and SDCBP2 remained near the top under all deletions. GZMK and CD8A preserved their negative direction but moved outside the top 200 in some deletions. Among stress-response candidates, MT1E and GPX3 had high weight and deletion frequencies despite high $I^2$, showing that a stable composite rank can coexist with variable effect magnitude. The deletion analysis therefore supports a graded candidate resource rather than an invariant list.

3.5. Pathway Direction Recurred Across All Five Study Deletions

The full meta-analytic ranking mapped to 14,408 unique bovine Entrez identifiers. GO biological-process GSEA returned 215 terms at FDR $<0.05$, and KEGG GSEA returned 51 (Figure 5). Positive GO shifts included electron transport chain (normalized enrichment score [NES] 3.28, FDR $6.09\times {10}^{-10}$), oxidative phosphorylation (NES 3.13, FDR $7.02\times {10}^{-9}$), and translation (NES 2.96, FDR $3.06\times {10}^{-19}$). Protein folding (NES 2.56, FDR $1.07\times {10}^{-6}$) and response to heat (NES 2.53, FDR $2.49\times {10}^{-4}$) also shifted positively. KEGG confirmed the same axis through ribosome, oxidative phosphorylation, proteasome, endoplasmic-reticulum protein processing, and thermogenesis (all FDR $<{10}^{-5}$).

Negative GO shifts included regulation of small-GTPase signaling (NES $-2.70$, FDR $5.64\times {10}^{-8}$), T-cell receptor signaling (NES $-2.11$, FDR 0.0190), and antigen processing and presentation (NES $-2.12$, FDR 0.00471). Regulation of leukocyte-mediated cytotoxicity (NES $-2.09$, FDR 0.0148) and cell-cycle process (NES $-1.80$, FDR $2.06\times {10}^{-4}$) also shifted negatively. KEGG similarly identified Th1/Th2 differentiation (NES $-2.38$, FDR $5.75\times {10}^{-6}$) and cell-adhesion molecule interactions (NES $-1.90$, FDR 0.00121). Disease-named KEGG labels reflected shared stress modules and were not interpreted as disease evidence.

With 128 genes from the top 200 mapped to Entrez identifiers, ORA provided a focused complement to GSEA. Up-prioritized genes were enriched for zinc-ion homeostasis, detoxification, and antimicrobial or humoral responses; gluconeogenesis and hormone biosynthesis were also represented. The down-prioritized set was enriched for lymphocyte and natural-killer-cell immunity, leukocyte cytotoxicity, antigen presentation, adhesion, and cell-cycle terms. Positive humoral stress-defense terms therefore contrasted with negative lymphocyte, cytotoxic, and antigen-presentation programs.

Pathway-level recurrence was high under study deletion. Of the 266 full-data FDR-significant pathways, 259 retained the same direction in all five deletions and the remaining seven did so in four. Full-data NES and median deletion-specific NES had a Spearman correlation of 0.972 (Supplementary Figure S5; Supplementary Tables S23 and S24). Twenty-one pathways remained FDR $<0.05$ in every deletion. Endoplasmic-reticulum protein processing retained positive significance in all five, while small-GTPase regulation, T-cell receptor signaling, cell-cycle process, and Th1/Th2 differentiation retained negative significance in all five. Translation, oxidative phosphorylation, protein folding, proteasome, response to heat, antigen presentation, and leukocyte cytotoxicity preserved direction in every deletion but crossed the FDR threshold in one or two deletions.

4. Discussion

Cross-study reproducibility in bovine heat stress was strongly scale dependent. Modified Knapp–Hartung inference yielded no transcriptome-wide FDR-significant genes. At the same time, a small group of candidates retained direction and showed recurrent prioritization across sensitivity analyses, and 259 of 266 FDR-significant pathways preserved direction in every study deletion. The evidence therefore forms a hierarchy: coordinated programs were the most portable signal, a smaller candidate set supported targeted follow-up, and the wider gene ranking remained context dependent. This scale-dependent recurrence is the principal finding.

The difference between 19 normal-approximation findings and no modified Knapp–Hartung findings reflects uncertainty calibration rather than a reversal of the estimated effects. Both analyses used the same study effects and REML heterogeneity estimates, but only two to five studies informed each average. Modified Knapp–Hartung inference is particularly useful when few studies have unequal precision [15]. Prediction intervals address whether a comparable future true effect could change direction [16]. Thirty-four primary plug-in intervals excluded zero, 15 did so with the alternative $t_{k-2}$ critical value, and 23 of the primary exclusions had ${\widehat{\tau }}^2=0$. Together, these analyses define the gene-level resolution supported by the current evidence.

The study design extends conventional livestock meta-analysis, which typically combines phenotypic outcomes such as milk yield, feed intake, fertility, or physiological indicators across experiments [1]. Here, thousands of genes were treated as outcomes, their effects and precision were re-estimated from uniformly processed raw reads, and each original study design was preserved before integration. Heat effects were defined within study, and cross-study synthesis used their effect sizes and precision. The analysis also differs from candidate-gene and GWAS reviews, which catalogue reported loci without recomputing transcriptomic effects [3,5]. Repeating enrichment after each study deletion further tested whether the biological interpretation itself depended on a single dataset. The resulting effect-size evidence layer is reusable for targeted expression assays, tissue-matched challenge studies, and future integration with physiological phenotypes or thermotolerance GWAS.

The positive pathway axis connected protein quality control with cellular energy demand. Heat exposure increases protein-folding and repair requirements, which is consistent with positive shifts in chaperone activity, endoplasmic-reticulum processing, proteasome function, translation, and heat-response genes. Mammary and cell-culture studies have independently reported heat-shock and proteostasis programs [8,24], while bovine myocyte transcriptomic and proteomic data show temperature-dependent remodeling of heat-shock, muscle, and metabolic pathways [25]. In the present analysis, endoplasmic-reticulum protein processing remained positively enriched at FDR $<0.05$ in all five study deletions. Translation, oxidative phosphorylation, protein folding, proteasome, and response-to-heat terms preserved their positive direction throughout. Tissue-specific liver work has also reported suppressed oxidative phosphorylation under particular heat and intake contexts [9]. The combined evidence is compatible with context-dependent mitochondrial remodeling in which individual metabolic genes vary, while energy metabolism remains a recurrent component of the heat response.

The negative pathway axis comprised T-cell signaling, antigen presentation, leukocyte-mediated cytotoxicity, adhesion, and selected cell-cycle and small-GTPase programs. These directions were preserved in every study deletion, and several remained FDR-significant throughout. Conversely, up-prioritized candidates recovered antimicrobial and humoral responses. Heat stress can combine inflammatory signaling with impaired cellular immune competence, depending on tissue, duration, endocrine state, and the compartment measured [2,4,26]. Longitudinal work in dairy-cow blood and peripheral blood mononuclear cells has reported transient inflammatory activation alongside reduced T-cell signaling during chronic heat exposure [13]. Multi-omics analysis has also linked immune regulation to metabolic heat adaptation through CIITA-dependent programs [14]. The present bulk, multi-tissue result is consistent with a selective negative shift in lymphocyte and cytotoxic signaling within broader immune reprogramming.

The pathway results provide context for the leading genes. IL1R2 was positive in all five studies, had low heterogeneity, retained a prediction interval above zero, and ranked within the top four in every study deletion. IL1R2 is a non-signaling decoy receptor that restrains IL-1 activity, so its induction is consistent with feedback regulation of IL-1 signaling [27]. In contrast, CD8A and GZMK were negatively prioritized alongside reductions in T-cell and cytotoxic programs. GZMK marks cytotoxic and inflammatory lymphocyte states in several systems [28]; lower GZMK and CD8A in bulk bovine data may reflect altered immune-cell activity, abundance, trafficking, or a mixture of these processes. SDCBP2 was positive and stable across four contributing studies but remains less developed in bovine heat-stress literature. Its recurrent ranking supports mechanistic follow-up.

Study deletion and HSTSI definition sensitivity determine how the candidate resource should be used. All five deletions preserved the direction of estimable genes in the full-data top 200, but top-200 overlap fell as low as 68 and the median Jaccard index was 0.338. Alternative HSTSI definitions produced top-200 overlaps of 24–200, with effect magnitude exerting the greatest influence. Exact membership is therefore less portable than direction. IL1R2, SDCBP2, and GZMK nevertheless remained in the top 200 under all ten score definitions, while canonical heat-shock genes ranked lower because their magnitudes varied sharply across studies. These profiles support validation panels that combine broadly recurrent candidates such as IL1R2 with context-sensitive candidates such as GZMK, MT1E, or GPX3.

The breadth of the analysis defines its scope. Average effects describe recurrence across the included tissues, breeds, life stages, and exposure designs, while tissue-specific mechanisms and cell-composition effects remain questions for focused studies. Functional annotation covered 128 genes from the top 200, and the five-study evidence base provided greater resolution for pathway direction than for gene-specific heterogeneity or prediction intervals. These properties support the intended use of the resource: pathway maps for shared biology and graded candidates for follow-up. Prospective cohorts, thermotolerance phenotypes, genetic evidence, and functional experiments can now test which candidates generalize within specific biological contexts.

5. Conclusions

Study-aware integration of 107 bovine heat-stress RNA-seq libraries showed that cross-study recurrence was strongest at the level of coordinated biological programs. Translation, protein quality control, mitochondrial energy, and heat-response programs shifted positively, whereas lymphocyte signaling, cytotoxicity, antigen presentation, adhesion, and selected cell-cycle and small-GTPase programs shifted negatively. A graded candidate set led by IL1R2, SDCBP2, GZMK, and CD8A showed greater directional than rank stability. By separating formal significance, candidate stability, and pathway recurrence, the study provides a calibrated route from heterogeneous public transcriptomes to tissue-matched validation of the bovine heat response.