Macroevolutionary changes in gene expression response to an immune stimulus across the diversity of fishes

TL;DR

This study shows that while the immune response to a stimulus is similar in many fish species, the genes involved vary widely, suggesting the need for diverse models to understand immunity.

Contribution

The study reveals macroevolutionary flexibility in the gene regulatory architecture of immune responses across ray-finned fishes.

Findings

The fibrosis response to an immune challenge is conserved across 14 fish species.

Transcriptomic responses to the immune stimulus are inconsistent across species.

Few gene orthogroups show consistent differential expression across all species.

Abstract

Our understanding of the vertebrate immune system is dominated by a few model organisms such as mice. This use of a few model systems is reasonable if major features of the immune systems evolve slowly and are conserved across most vertebrates, but may be problematic if there is substantial macroevolutionary change in immune responses. Here, we present a test of the macroevolutionary stability, across 14 species of ray-finned fishes, of the transcriptomic response to a standardized immune challenge. Intraperitoneal injection of an immune adjuvant (alum) induces a fibrosis response in nearly all jawed fishes, which in some species contributes to anti-helminth protection. Despite this conserved phenotypic response, the underlying transcriptomic response is highly inconsistent across species. Although many gene orthogroups exhibit differential expression between saline versus alum-injected fish in at least one species, few orthogroups exhibit consistent differential expression across species. This result suggests that although the phenotypic response to alum (fibrosis) is highly conserved, the underlying gene regulatory architecture is very flexible and cannot readily be extrapolated from any one species to fishes (or vertebrates) more broadly. The vertebrate immune response is remarkably changeable over macroevolutionary time, requiring a diversity of model organisms to describe effectively.

Figures

a) Species tree inferred from all genes (STAG) with root inferred using duplication events (STRIDE) with branch lengths equivalent to the average number of substitutions per site across gene families. b) Percentage of transcripts from each species assigned to orthogroups (OG). c) Experimental diagram illustrating intraperitoneal injections for the PBS (control) and alum stimulus (immune adjuvant) for four representative species (top to bottom: O. mykiss, G. aculeatus, D. rerio, and A. mexicanus). d) Summed orthogroup (OG) normalized expression for an orthogroup showing a main effect of species (OG0006036) using DESeq2's median of ratios method. Each point represents treatment-specific OG summed expression for each sample within each species, where solid lines represent a significant treatment effect (P < 0.1) and dotted lines represent non-significant treatment effect (P > 0.1). Created with BioRender.com.

a) For each orthogroup summed normalized gene expression, we used a linear model to estimate the effect size (η2) for the effects of injection treatment (PBS versus alum) and the treatment × species interaction on variation in the focal trait. Treatment η2 indicates the extent to which any orthogroup summed normalized expression diverges predictably from PBS to alum. The η2 for the treatment × species interaction measures the extent to which treatment effect depends on the species for every orthogroup. The dashed line is a 1:1 line, for ease of visualization; points falling above this line have a larger treatment effect than interaction effect. Orthogroup expression for PBS and alum injections across all species where solid lines represent a significant treatment effect (P < 0.1) and dotted line represents non-significant treatment effect (P > 0.1). b to d) Standard normal (z-score) expression for b) single-copy ortholog (OG0015527), c) orthogroup with highest effect size of treatment (OG0000016), and d) orthogroup with highest effect size of the interaction of treatment and species (OG0003530).

a) Differential orthogroup expression Wald test statistic correlation between T. ocellicauda and O. mykiss with a linear fit model and c) GO enrichment odds ratio correlation between A. mexicanus and C. viridis with the Pearson's correlation coefficient (PCC) and linear model inferred intercept. Pairwise species correlation distributions for b) differential orthogroup expression and d) GO enrichment.

Heatmap representing the top 30 differentially expressed orthogroups with the dendrogram visualizing the hierarchy of expression data (upper), and the species tree inferred from all orthogroups (lower) with lines illustrating the co-phylogeny between dendrogram and species tree.

a) The distribution of orthogroup differential expression phylogenetic signal categorically grouped by negligible (λ < 0.05), weak (0.05 < λ < 0.2), moderate (0.2 < λ < 0.4), and strong (λ > 0.4) signal. b) Model comparison (ΔDIC) between multivariate generalized linear mixed models with and without a random effect of phylogenetic relatedness. c) Phylogenetic signal for 20 OGs with the highest λ estimate.

Full text

Introduction

Advancements in immunology have largely been driven by research on a few model organisms including mice (Gros and Casanova 2023), nonhuman primates (Bjornson-Hooper et al. 2022), and zebrafish (Liu et al. 2024). While these models have been transformative in our understanding of the immune system, their application to humans relies on the evolutionary stability of the immune process across tens to hundreds of million years of evolutionary divergence. Emerging evidence suggests that the structure and function of the vertebrate immune system vary greatly, at several scales of biological organization. Immune variation can be observed among individuals within a given population. For example, severity of SARS-CoV2 is genotype dependent; individuals homozygous for the e4 allele at ApoE gene display increased mortality (Kuo et al. 2020; Ostendorf et al. 2022). This genetic variation within populations can then be a target for natural selection, enabling rapid evolutionary change within populations through time, driving genetically based immune differentiation between populations and thus species. For instance, the frequency of alleles at major histocompatibility complex (MHC) loci diverges between animal populations subjected to different parasite communities (Wegner et al. 2003; O’Connor et al. 2018). These differences have been thought to contribute to the evolution of reproductive isolation and speciation. As a result, immune genotypes and functions can differ appreciably between even closely related species. Given this documented rapid evolution of immunity within and between closely related species, it seems improbable that immune function should be highly conserved across deeper evolutionary time. At a minimum, we should seek to understand what features of immunity are conserved and what features are divergent across animals (e.g. Gilbertson and Weinmann (2021)).

While vertebrate immune systems mostly overlap across diversity, within fishes, some once considered essential features of adaptive immunity have turned out to be surprisingly dispensable. MHC emerged roughly 500 million years ago at the origin of chordates and is one of the most studied immune genes. MHC is present in nearly all vertebrates, appearing as a fundamental and canonical aspect of adaptive immune function. Yet, there is substantial turnover in the number and identity of MHC paralogs between closely related primates (Lyn Fortier and Pritchard 2025). Sometimes, whole categories of MHC are lost, such as the non-functionalization of MHC Class II in Atlantic cod (Gadus morhua) (Star et al. 2011) and its close relatives (Malmstrøm et al. 2016). Even more dramatically, angler fish have lost many components of their adaptive immune system, enabling an evolutionary transition to extreme size sexual dimorphism (Brownstein et al. 2024) and male reproductive parasitism (Dubin et al. 2019; Swann et al. 2020). It is essential to understand whether these examples of immune function loss (or change) are exceptions, restricted to a few lineages, or are typical of other immune functions and gene families.

Such changes in immune function may prove to be common across vertebrates, but there have been few systematic evaluations of immune function across a diverse swath of vertebrate species. Comparative immunology has been limited in part by the difficulty of studying immune responses in many species simultaneously, especially considering the need for species-specific reagents such as monoclonal antibodies to measure many immune traits. And yet, the advent of genome and transcriptome sequencing has facilitated high throughput comparisons of gain or loss of genes across many species' genomes, or changes in transcriptional responses (e.g. (Brownstein et al. 2024)). Comparative transcriptomics is emerging as an especially promising tool for comparative immunology (Chen et al. 2023), without the need to develop species-specific molecular probes. One can measure differential expression of genes within each of many species, subjected to a standardized immune challenge. We can then compare differential expression among many diverse species, to ask: is the immune response highly conserved across most vertebrates? Or, is the immune response evolving so quickly that each species responds to a given immune challenge with a unique set of transcriptional changes? Almost certainly, the answer lies somewhere between these extremes, and comparative transcriptomics can let us quantify the extent to which immune responses have a shared, or unique, transcriptional basis across a wide range of animals.

Here, we use comparative transcriptomics to test whether a sample of fish species from across the full span of the teleost tree of life exhibit a similar response to a standardized immune challenge. Many species of fish (and other vertebrates) develop extensive fibrosis during an inflammatory immune response (Vrtílek and Bolnick 2021). Fibrosis involves the formation of an extensive extracellular matrix of protein fibers (collagen, fibrinogen), produced by an evolutionarily ancient cell type (fibroblasts). Fibrosis is essential in wound healing, but is also a common feature of inflammation (Wynn and Ramalingam 2012). This fibrosis response can be beneficial: in some fish species (e.g. three-spined stickleback), fibrosis has been found to confer protection against intraperitoneal tapeworm infection (Lohman et al. 2017; Weber et al. 2022). Fibrosis is also widely considered to be a pathological side effect of inflammation and is involved in a variety of human diseases (Henderson et al. 2020). In stickleback, the immune benefits of fibrosis are countered by reduced fecundity (De Lisle and Bolnick 2021) and reduced mobility (Matthews et al. 2023). Consequently, fibrosis responses evolve rapidly within, and diverge between, stickleback populations (Hund et al. 2022), which exhibit microevolutionary changes in fibroblast activation genes (Fuess et al. 2021).

In stark contrast to this rapid evolution of fibrosis in stickleback, another study found that fibrosis responses are highly conserved across teleosts (Vrtílek and Bolnick 2021). Intraperitoneal injection of an immune adjuvant (alum, widely used in vaccines to stimulate inflammation) causes a strong fibrosis response within the body cavity of stickleback (Hund et al. 2022). Extending this assay to 17 species of fish sampled from across the teleost phylogeny, Vrtílek and Bolnick (2021) found that nearly all species initiate fibrosis after alum exposure. This result reveals a puzzling contrast: fibrosis is highly conserved among teleost species exposed to alum spanning over a hundred million years of diversity (Vrtílek and Bolnick 2021), yet fibrosis, as employed in parasite defense, functionally diverges between stickleback populations separated for a mere 10,000 years (Weber et al. 2022). This contrast exemplifies the broader contradiction of rapid immune evolution, and conservation of immune features, in vertebrates. Usefully, we can apply a standardized immune challenge (alum), which elicits a standard phenotypic response (fibrosis) in most fish species. Here, we pose the question: is this similar phenotypic response to alum generated by a standard transcriptomic response, shared across diverse teleosts? Or, do different fish species respond to alum in unique ways, indicating significant changes in the transcriptomic pathways responding to an immune challenge and driving fibrosis?

To answer these questions, we generated reference transcriptomes for 14 species from the Vrtílek and Bolnick (2021) study. We sequenced the head kidney (pronephros) transcriptomes of the experimental fish injected with either saline (controls) or alum solution, using head kidneys archived from the Vrtílek and Bolnick (2021) experiment. Within each species, we test whether each transcript family (orthogroup) exhibits differential expression in response to alum injection (relative to saline controls). Then comparing across species, we test whether this differential response is shared among species, or varies among species. That is, do we observe a treatment by species interaction effect in which an orthogroup is up- or down-regulated in response to alum in some species, but not others (or responds in opposite directions depending on the species). We then examine the most differentially expressed transcripts across any teleosts, to identify what biological functions are responding to alum. For instance, do we see consistent overrepresentation of fibrosis-related or inflammation-related genes, even if the specific orthogroups are inconsistent from species to species? Are the changes in expression, among species, associated with the extent to which species develop fibrosis after alum injection? Considering interpretation limitations, immune responses are temporally dynamic within species (e.g. Soto et al. (2022)) and vary greatly across tissue type (Domínguez Conde et al. 2022). Here we avoid cross-tissue differences by sampling a single immune organ (head kidney) at a single time point for each species to get a transcriptomic snapshot into the transcriptional profile of a conserved inducible immune phenotype. We show that each species exhibits a transcriptional response to alum (even the few species that do not respond with fibrosis). But, this transcriptional response is largely unique to each species, at the level of transcript orthogroups. However, at the level of gene function (e.g. Gene Ontology [GO] category), we see some weakly similar responses across species separated by hundreds of millions of years of evolutionary divergence.

Results

The de novo head kidney-specific transcriptomes assemblies from 14 diverse ray-finned fish were used to infer the phylogenetic relationship among the experimental species and the resultant tree is congruent with the current ray-finned fishes phylogeny (Hughes et al. 2018) (Fig. 1a). Further, 65% to 92% of the assembled transcripts for each species were assigned to an orthogroup (Fig. 1b). The comparative approach identified only 20 single-copy orthologs across all 14 species. In contrast, 7,528 orthogroups had at least one transcript from all species, out of 73,286 total orthogroups identified, with a mean orthogroup size of 10.8 distinct transcripts. This paucity of single-copy orthologs is unsurprising given that our phylogenetic sampling spans at least two major whole genome duplication events.

a) Species tree inferred from all genes (STAG) with root inferred using duplication events (STRIDE) with branch lengths equivalent to the average number of substitutions per site across gene families. b) Percentage of transcripts from each species assigned to orthogroups (OG). c) Experimental diagram illustrating intraperitoneal injections for the PBS (control) and alum stimulus (immune adjuvant) for four representative species (top to bottom: O. mykiss, G. aculeatus, D. rerio, and A. mexicanus). d) Summed orthogroup (OG) normalized expression for an orthogroup showing a main effect of species (OG0006036) using DESeq2's median of ratios method. Each point represents treatment-specific OG summed expression for each sample within each species, where solid lines represent a significant treatment effect (P < 0.1) and dotted lines represent non-significant treatment effect (P > 0.1). Created with BioRender.com.

For many orthogroups, expression varied among the species in our samples (96.76% of orthogroups show a main effect of species, P < 0.1). For OG0006036, the main effect of species showed the highest effect size, while the main effect of treatment and the interaction between species and treatment was non-significant (Fig. 1d; Type II ANOVA, treatment F1,130 = 0.222, P > 0.05; species F13,130 = 297.254, P < 0.001; treatment × species F13,130 = 0.609, P > 0.05). To illustrate, under control saline injections, the mean summed expression for OG0006036 was 23.9 (SE = 2.65) for A. mexicanus while in C. viridis, mean summed expression was 128 (SE = 16.3) (Fig. 1d). Few orthogroups consistently respond to the alum injection treatment (20.34% of orthogroups had a significant main effect of alum treatment, P < 0.1). As an example, OG0011385 log2 fold change (lfc) (comparing alum injected vs control fish in each species) ranged from −0.98 (in S. fasciatus) to −0.05 (in A. mexicanus) indicating a moderately consistent down-regulation of this gene in alum-injected fish, though the magnitude of this change varied among species (Figure S1a; Type II ANOVA, treatment F1,130 = 12.849, P = P < 0.001; species F13,130 = 25.409, P < 0.001; treatment × species F13,130 = 1.571, P > 0.05). The top blastp result for OG0011385 is sap18, a component of the SIN3-repressing complex which can direct the formation of a repressive complex to core histone proteins (Singh et al. 2010). Lastly, when considering the interactive effect of treatment and species, the treatment effect was largely species dependent.

Instead of a uniform transcriptomic response, we observe species-specific transcriptomic responses to alum where the treatment × species interaction effect sizes are consistently larger than the treatment effect size (Fig. 2a). This species-specific transcriptomic response was observed even for a single-copy ortholog (OG0015527) where the interactive effect of treatment and species explained 2.28% of the gene expression variation while treatment only explained 0.28% of the variation in gene expression (Fig. 2a and b; Type II ANOVA, treatment F1,130 = 1.55, P = 0.216; species F13,130 = 32.386, P < 0.001; treatment × species F13,130 = 0.959, P = 0.495). Even for the orthogroup with the highest amount of variation explained by treatment (OG0000016, hemoglobin subunit alpha), treatment only explained 9.4% of the variation while the interactive effect explained a larger proportion of variation (19.3%) illustrated by the species-dependent slopes in response to alum injection (Fig. 2c; Type II ANOVA, treatment F1,130 = 24.98, P < 0.001; species F13,130 = 6.096, P < 0.001; treatment × species F13,130 = 3.943, P < 0.001). OG0003530 had the greatest amount of expression variation explained by the interactive effect of alum injection and species where some species decreased expression in response to alum injection while others increased expression (Fig. 2d; Type II ANOVA, treatment F1,130 = 43.628, P < 0.001; species F13,130 = 98.867, P < 0.001; treatment × species F13,130 = 93.57, P < 0.001) indicating a large species-specific expression response to alum injection. Therefore, the response to alum is highly heterogeneous among fish species indicating there is no single gene regulatory program responding to alum consistently across all species. Similarly, there is no single gene regulatory program associated with high versus low fibrosis within species (Figure S2).

a) For each orthogroup summed normalized gene expression, we used a linear model to estimate the effect size (η2) for the effects of injection treatment (PBS versus alum) and the treatment × species interaction on variation in the focal trait. Treatment η2 indicates the extent to which any orthogroup summed normalized expression diverges predictably from PBS to alum. The η2 for the treatment × species interaction measures the extent to which treatment effect depends on the species for every orthogroup. The dashed line is a 1:1 line, for ease of visualization; points falling above this line have a larger treatment effect than interaction effect. Orthogroup expression for PBS and alum injections across all species where solid lines represent a significant treatment effect (P < 0.1) and dotted line represents non-significant treatment effect (P > 0.1). b to d) Standard normal (z-score) expression for b) single-copy ortholog (OG0015527), c) orthogroup with highest effect size of treatment (OG0000016), and d) orthogroup with highest effect size of the interaction of treatment and species (OG0003530).

The number of distinct transcripts assigned to orthogroups (e.g. gene family size) varied among species, and there was a weak negative correlation between the absolute value of lfc in gene expression (alum vs saline) and the number of transcripts assigned to an orthogroup (Pearson's product-moment correlation = −0.04, t = −11.51, P < 2.2e-16). This indicates copy number minimally influences the estimated differential expression response to injection based on summing the expression for all transcripts present in a species for a given orthogroup.

Despite the lack of conservation at single orthogroups, there are some transcriptome-wide similarities between species' responses to alum. On average across the transcriptome, orthgroups, which are up-(or down-)regulated in one species, tend to be up- or down-regulated in other species. Specifically, there is a correlation between lfc associated with alum, when comparing pairs of species (Fig. 3a and b). For example, the complete vector of log2 fold changes in T. ocellicauda and the corresponding lfc vector for O. mykiss are positively correlated (Pearson's correlation coefficient = 0.12, t = 5.16, P = 2.71e-07) (Fig. 3a). When considering the distribution of species-pairwise correlations for the Wald test statistic used to estimate the treatment effect, we find the correlation estimates are largely positive and the distribution mean is greater than zero (t-test, t = 9.44, P = 4.23e-15) (Fig. 3b). This correlation between orthogroup differential expression varies in strength depending on the pair of species being examined (Fig. 3b). Some species pairs are more strongly correlated, implying more similar transcriptome-wide response to alum, than other pairs of species whose response is uncorrelated (O. niloticus and T. ocellicauda, Pearson's correlation coefficient = 6.94e-05, t = 0.004, P = 0.99), or negatively correlated (A. mexicanus and S. fasciatus Pearson's correlation coefficient = −0.1, t = −4.49, P = 7.45e-06). Yet, overall there is a tendency for a similar general response to alum among species, though this is not observed for individual orthogroups. Further, we might expect species which are more closely related to show similar patterns of orthogroup differential expression, yet there is discordance between the hierarchical grouping of differential orthogroup expression and the inferred species tree (Fig. 4; Icong = 0.9212, P = 4.4676) (de Vienne et al. 2007). The evolution of differential orthogroup expression does not follow species diversification wherein species changes occur independently and may be contained to specific orthogroups with accelerated expression changes (Figs 3 and 4).

a) Differential orthogroup expression Wald test statistic correlation between T. ocellicauda and O. mykiss with a linear fit model and c) GO enrichment odds ratio correlation between A. mexicanus and C. viridis with the Pearson's correlation coefficient (PCC) and linear model inferred intercept. Pairwise species correlation distributions for b) differential orthogroup expression and d) GO enrichment.

Heatmap representing the top 30 differentially expressed orthogroups with the dendrogram visualizing the hierarchy of expression data (upper), and the species tree inferred from all orthogroups (lower) with lines illustrating the co-phylogeny between dendrogram and species tree.

A similar result is observed for the enrichment of GO terms for differentially expressed transcripts (Wald test, P > 0.1) where there are some GO terms which are similarly enriched between species. As an example, A. mexicanus and C. viridis GO enrichment show a positive correlation (Pearson's correlation coefficient = 0.33, t = 29.86, P < 2.2e-16) (Fig. 3c). Overall, species-pairwise correlation estimates for GO enrichment are positive (t-test, t = 26.1, P < 2.2e-16) (Fig. 3d). Therefore, the response to the immune stimulus evokes specific ontology terms which are utilized in mounting the immune response to the alum treatment.

The phylogenetic signal (Pagel's λ) for OG differential expression ranged from 0.002 to 1 with a median value of 0.998, but the 95% confidence intervals (CIs) associated with the λ estimate on average had a range of 0.879; therefore, we present the lower CI (Fig. 5a). For 65% of the OGs, the mixed effects generalized linear model without accounting phylogenetic relatedness had a lower DIC than the model including the random effect of phylogenetic relation, and those OGs with a higher λ had a lower DIC for the phylogenetic models (Fig. 5b). Only 1.6% OGs had a positive posterior variance, lower DIC for the phylogenetic model, and showed a relatively strong phylogenetic signal (λ > 0.4). For the OGs with the highest λ estimate (Fig. 5c), blastp results indicated the OGs with the most phylogenetically conserved differential expression are involved in DNA replication (OG0003653 lysine methyltransferase 5A; OG0011136, MCMBP; OG0011136, MCMBP), mitochondrial function (OG0015672, mitochondrial matrix import factor 23; OG0009765, ATP synthase subunit d; OG0003938, PPARGC1A-related coactivator 1), and cell division (OG0001348, tropomyosin 4a; OG0014197, septin 8a; OG0007722, SKA3). OG0009255, an immune related OG, encodes for macrophage migration inhibitory factor which promotes the movement and recruitment of leukocytes to sites of infection and inflammation (Bernhagen et al. 2007).

a) The distribution of orthogroup differential expression phylogenetic signal categorically grouped by negligible (λ < 0.05), weak (0.05 < λ < 0.2), moderate (0.2 < λ < 0.4), and strong (λ > 0.4) signal. b) Model comparison (ΔDIC) between multivariate generalized linear mixed models with and without a random effect of phylogenetic relatedness. c) Phylogenetic signal for 20 OGs with the highest λ estimate.

Discussion

When comparative biologists observe a phenotype present across most or all members of a clade, they are inclined to conclude that the trait is homologous and conserved (Wiley and Lieberman 2011). Principles of parsimony (likelihood models of trait evolution) would suggest that the trait evolved prior to the most recent common ancestor of the clade and has been maintained as the clade diversified. But, does this phenotypic homology imply that the proximate genetic, regulatory, and developmental processes are also homologous? Or, might homology be superficial? Complex biological systems are often characterized by many-to-one mapping rules. That is, there may be many configurations of parts (different bone lengths, cytokine levels, or gene expression levels) that synergistically yield a single outcome. For instance, the nonlinear properties of four-bar levers imply that many possible fish skull and jaw shapes can yield identical kinematic force transmission for biting or suction feeding (Wainwright et al. 2005). A function (e.g. bite force) could thus be kept constant through evolution (appearing to be homologous), even as its proximal mechanism (skull shape) diverges between species (Alfaro et al. 2005).

The fibrosis immune response is largely conserved across ray-finned fish. Intraperitoneal injection of an immune adjuvant, alum, was recently found to induce fibrosis in three-spined stickleback (Hund et al. 2022). Applying this assay to 17 species of actinopterygian fishes from a wide diversity of sub-clades, Vrtílek and Bolnick (2021) showed that nearly all fishes retain the ability to initiate peritoneal fibrosis (Vrtílek and Bolnick 2021) with the exception of Otophysa species: the Mexican tetra (A. mexicanus) and the bleeding-heart tetra (H. erythrostigma). This result suggests that the fibrosis immune response (commonly associated with helminth infection) evolved in, or before, the most recent common ancestor of the Actinopterygii and hence most vertebrates. Fibrosis appears to be a highly conserved homologous trait that has been sustained across more than 350 million years of evolutionary divergence. Our results, however, show that the underlying gene expression program is far less homologous.

The conserved inducible fibrosis phenotype provides an opportunity to determine how the associated gene expression response evolves in a phylogenetically informed approach. Specifically, we characterize how associated differential expression is evolving across ray-finned fishes diversity (Hughes et al. 2018) in 14 of the experimental species from Vrtílek and Bolnick (2021). Here, we assigned transcripts to orthogroups which have a shared evolutionary history, revealing the predominant response to the immune stimulus is species-dependent with limited shared treatment effects. Most species do respond to alum injection and most initiate fibrosis (Vrtílek and Bolnick 2021), but each species uses a unique set of expression changes to achieve the apparently homologous phenotypic response. However, orthogroup differential expression was positively correlated across species. GO enrichment was also positively correlated, which together suggest each species is evolving independently while maintaining the ability to induce peritoneal fibrosis, yet there is a degree of functional constraint present with coordinated enriched transcript ontologies. While individual transcripts or orthogroups might be used differently from species to species, across the transcriptome as a whole, similar families of genes are nevertheless being engaged, for some (but not all) combinations of species being compared. In light of the omnigenic model (Boyle et al. 2017), for complex, phenotypes like fibrosis, large portions of the genome likely interact pleiotropically to generate the fibrosis phenotype, meaning the same set of genes that exist in these species are not all being induced or expressed consistently which could be due to pleiotropic effects that results in the same fibrosis phenotype.

As inferred from the differential transcriptomic response to alum, species maintain the ability to induce fibrosis through differing gene regulatory mechanisms, indicating a functional evolutionary constraint, yet multiple different expression pathways to arrive at the conserved phenotype. This result is contrary to a common assumption in comparative biology. For example, even with a strong selective force like hydrogen sulfide, with documented convergent evolution in morphological and physiological phenotypes in three sympatric fish lineages, the patterns of genomic divergence relative to sulfide tolerance are non-convergent with no consistent differentiated genomic regions (Greenway et al. 2024). The varied way by which fish arrive at the same fibrosis phenotype could be attributed to selection on specific components of the highly complex fibrosis gene regulatory network. Stabilizing selection on the ability to employ fibrosis can reflect its general importance as a physiological function, involved in many biological processes including development, wound healing, and immunity. But this pleiotropy may mean that different species, while maintaining fibrosis, evolve to deploy it in different contexts, prioritizing different subsets of the pleiotropic functions. For instance, the ability to induce fibrosis in response to flatworm infection is heterogeneously maintained in stickleback lake populations separated by only a few kilometers (Weber et al. 2022), facing similar parasite encounter rates. While in other species, fibrosis may be deployed more in wound-healing, or in embryonic development. These divergent applications of fibrosis may select for divergent gene regulatory pathways, producing the diversity of transcriptomic responses that we observe, even as fibrosis responses are retained. Therefore, even if the phenotypes under selection are conserved, the mechanism by which organisms arrive at the phenotype may not be conserved. As the gene regulatory architecture which gives rise to the fibrosis phenotype is generally inconsistent across the ray-finned fishes, this indicates that evolution is exploring neutral changes in pathways which retain the same function, evoking the concept of systems drift (or developmental systems drift) (True and Haag 2001; Schiffman and Ralph 2022).

While our study provides a snapshot into the gene regulatory divergence stimulated by a common immune challenge across macroevolutionary divergence in fishes, the immune response is temporally dynamic (Wuitchik et al. 2019) and that dynamism differs across tissues and even among cell types of the same tissue (Colquitt et al. 2021). Here, we focus on a single tissue, the head kidney sampled at a single time point for each fish following alum injection immune challenge. Previous work informed our temporal sampling design wherein Hund et al. (2022) found alum induced the fibrosis immune response within ten days following alum injection in three-spined stickleback. Therefore, we sampled fish 5 days following alum injection for all but cold-water O. mykiss which was sampled ten days post injection. In teleosts, the head kidney is similar to mammalian bone marrow because the head kidney produces immune cells like macrophages (Zapata et al. 2006). Macrophages play a multifaceted role in immune function from chemotaxis to antigen trapping (Agius and Roberts 2003); therefore, we used this immune representative tissue to estimate the transcriptional response to alum injection across fish diversity.

Overall, our experimental results represent the incongruence that can arise between the emergence of a conserved immune phenotype (here, fibrosis) and the gene regulatory mechanisms which are induced during the immune response. We demonstrate that vertebrate species responding to a common immune stimulus resulting in the same immune phenotype occurs through species-specific gene regulatory logic. Future cross-species transcriptomic integration may explore the assumption of constrained functional divergence among orthogroups. By examining expression connections among transcripts to determine if members of the same orthogroup are expressed with the same suite of transcripts will allow for inference about gene family evolution and functional diversification. Evaluating the consistency of these results across other vertebrates, particularly those employed in model immune research, highlights a need for phylogenetically diverse immune models systems. Moreover, the diverse gene regulatory structures suggest that no single species will be an effective genetic model for all human disease or physiological functions. A diverse portfolio of research organisms will better span the variety of immune regulatory systems in vertebrates.

Materials and methods

Fish selection and experimental methods

To gain an understanding of the evolution of gene expression associated with the formation of peritoneal fibrosis across ray-finned fishes (Actinopterygii), we intentionally selected species spread across the phylogenetic diversity of Actinopterygii which were small and commercially available. The 17 fish species were housed in laboratory aquaria and underwent experimental injections as described in (Vrtílek and Bolnick 2021). Because tapeworm proteins only limitedly induces fibrosis in ray-finned fishes (Vrtílek and Bolnick 2021), we focus on the gene expression response to alum, an immune adjuvant which activates the innate immune response (Kool et al. 2012) and broadly induces peritoneal fibrosis in Actinopterygii (Vrtílek and Bolnick 2021). Fish were peritoneally injected through their left flank with volume-mass corrected solutions at 20 µL per 1 g of average weight for a given species. We injected experimental fish with either 1× phosphate-buffered saline (PBS) or 1% alum (AlumVax phosphate) as our two treatments. These fish were euthanized five days following injection, but because stickleback and trout were raised at a lower temperature, they were euthanized ten days post injection. Following euthanasia, each sample was scored ordinarily for fibrosis which ranged between 0 and 3 where 0 indicates no fibrosis with freely moving internal organs and no fibrotic attachment to the peritoneal wall. Level 1 fibrosis is when internal organs adhere together and move as a unit. Level 2 fibrosis occurs when organs attach to the peritoneal wall but can be detached while fibrosis level 3 was the most extreme fibrosis phenotype where organ adhesion to the peritoneal wall could not be reattached without the peritoneal lining tearing apart and remaining on the internal organs when forcibly opened. This ordinal fibrosis score is described in Vrtlek and Bolnick (2021) and is similar to the ordinal classification of fibrosis presented in Hund et al. (2022). The ordinal scoring is highly repeatable among independent observers (Bolnick et al. 2024).

In addition to scoring fibrosis, the head kidney tissue was isolated and stored in RNAlater. We focused on characterizing the specific head kidney transcriptome response because in teleosts it is one of the tissues responsible for initializing the immune response, both innate and adaptive (Zapata et al. 2006). Of the 17 species used in the Vrtílek and Bolnick experiment, here we study 14 species whose head kidneys yielded sufficient quantities of high quality RNA and which had (or we could generate) reference transcriptomes from head kidneys.

Reference transcriptome sequencing

We extracted RNA from frozen head kidney samples from fish injected with tapeworm antigen homogenate (see Vrtílek and Bolnick (2021)) using MagMax96 Total RNA Isolation kit, following manufacturer's instructions. RNA yield was quantified using a Ribogreen fluorescent stain compared with a standard curve, with fluorescence measured on a Magellan 96-well plate fluorometer. While some species had sequences available from head kidney tissues (Table S1), there were species which lacked available head kidney sequences, so we generated de novo references using paired end 150 bp sequencing through Novogene on the Illumina NovaSeq platform (20 million reads per sample). We sequenced and assembled de novo transcriptomes for Polypterus senegalus, Hyphessobrycon erythrostigma, Tateurndina ocellicauda, Sphaeramia nematoptera, Chromis viridis, Salarias fasciatus, Nothobranchius furzeri, Lepomis macrochirus, Dichotomyctere nigroviridis, and Helostoma temminckii. We focused on head kidney transcriptomes because we wished to build references appropriate to the transcriptomes we later examine with 3′TagSeq, using the same tissues.

Transcriptome assembly and comparative transcriptomics

While some of the fish species had previously published transcriptomes, we wanted to ensure congruent assembly methods occurred for all species, so we independently de novo assembled transcriptomes for all species using previously published data where available and conducting full-length RNA-seq for a subset of species which had no head kidney RNA-seq data available (Table S2). De novo head kidney transcript assembly using previously published sequences included Danio rerio (Pasquier et al. 2016), Gasterosteus aculeatus (Huang et al. 2016), Oncorhynchus mykiss (Sun et al. 2023), and Oreochromis niloticus (Hou et al. 2023). For each species including those sequenced here and species with previously published data, we de novo assembled the head kidney-specific transcriptome using rnaSPAdes with paired end reads (Bushmanova et al. 2019). Then, to identify coding regions of the transcriptome we ran Transdecoder (Haas et al. 2013) which first identifies the longest open reading frame with a minimum length of 100 aa and then searches for protein domains with similarity to the Pfam-A Pfam database using HMMER3 (hmmer.org). Lastly, Transdecoder predicts candidate coding regions in the transcriptome.

The resultant Transdecoder (Haas et al. 2013) predicted peptide sequences for all species were used to identify orthologs, orthogroups, and infer a species tree based on duplication events through Orthofider2 (Emms and Kelly 2019) using DIAMOND (Buchfink et al. 2015) and the STRIDE (Emms and Kelly 2017) method to infer species tree. Previously, studies attempting to characterize the evolution of gene expression were limited by focusing on genes which were single-copy in all species. Moving beyond single-copy orthogroups to include orthogroup-level comparisons with species transcript content variability, we begin to understand how gene families with shared, potentially constrained function, evolve across phylogenetic diversity. For species comparisons, we retained orthogroups which had a minimum of one copy in all species.

3′TagSeq for differential expression of experimental fish

For samples from the injection experiment, we extracted RNA from head kidneys of saline injected and alum injected individuals of each species (sample sizes in Table S2). We used Ribogreen to quantify and normalize the extracted RNA to 20 ng/µL, and sent these to the Genome Sequencing and Analysis Facility at the University of Texas at Austin for 3′TagSeq (Table S3). 3′TagSeq a cost-effective approach that sequences just the 3′ ends, providing cheaper and more precise estimates of relative expression levels (albeit without information about splicing variants). 3′TagSeq libraries were constructed according to (Lohman et al. 2016, 2017) based on methods presented in (Meyer et al. 2011), aiming for 5 million reads on average per sample sequenced on the Illumina NovaSeq platform.

To generate read counts per transcript from the de novo head kidney transcriptomes assemblies, we used the coding sequence output from the Transdecoder predict algorithm. Using the pseudo-mapper Salmon (Patro et al. 2017) which maps and counts in a single step, we first indexed the transcriptome and subsequently quantified reads mapping to each transcript. After quantifying reads for each transcript, counts were imported into R (R Core Team 2024) and were normalized using the mean-of-medians approach described in DEseq2 (Love et al. 2014).

Orthogroup differential expression

For transcripts with a shared evolutionary history which were placed into orthogroups, the normalized transcript counts were summed for all transcripts which have membership in a given orthogroup for each sample individually. The resultant orthogroup count matrix was then analyzed using DEseq2 (Love et al. 2014) to estimate orthogroup differential expression and differential expression was plotted and grouped using heatmap R package (R Core Team 2024). To test for tree congruence between the heatmap hierarchical clustering and the inferred species tree, we implemented the method described in de Vienne et al. (2007).

Statistical analysis

Orthogroup expression was summed for all transcript members of a given orthogroup for each species, and then to determine the treatment effect for each orthogroup and species, we initially fit a linear model in R (lm). Then, we estimated the interaction effect for treatment and species for each orthogroup and then subsequently estimated the effect size (η2) for treatment and the interaction effect which allows us to determine whether species-specific treatment effects (treatment by species interaction) or shared treatment effects dominate variance in orthogroup expression.

To determine if differential orthogroup expression response to the alum treatment was consistent across species, for each orthogroup and species, we estimated the Wald test statistic using DeSeq2 (Love et al. 2014). Then to check whether the overall transcriptome-wide response was similar between species, we calculated correlations between all pairs of species, comparing the vectors of orthogroup Wald test statistics from each species.

Similarly, to determine if GO enrichment was consistent across all species, we first determined which transcripts were differentially expressed between treatments, within each species, using DeSeq2 (FDR adjusted P < 0.1) (Love et al. 2014). Then, for each de novo transcriptome assembly, we determine GO membership for each transcript using the web application TRAPID 2.0 (Bucchini et al. 2021). We then estimated the enrichment of differentially expressed transcripts for each GO term using a Fisher's exact test in R (fisher.test) estimated by an odds ratio test statistic. The resultant odds ratios for all GO terms were pairwise correlated for all species.

To determine the phylogenetic signal (Pagel's λ) for OG differential expression, we fit a mcmcglmm() with a random effect of the genetic distance variance-covariance matrix and a burn in of 10,000 chains iterated 600,000 times (Hadfield 2010). From the model, we estimated posterior variance, total variance, λ, λ's 95% CIs, and deviance information criterion (DIC). Simultaneously, we fit a reduced model without the random effect of phylogenetic relationship and estimated the DIC to compare fit of the reduced model with the full phylogenetic model.

Supplementary Material

msaf323_Supplementary_Data