
==== Front
Genome Biol Evol
Genome Biol Evol
gbe
Genome Biology and Evolution
1759-6653
Oxford University Press UK

39165136
10.1093/gbe/evae183
evae183
Article
AcademicSubjects/SCI01130
AcademicSubjects/SCI01140
Evolution of Key Oxygen-Sensing Genes Is Associated with Hypoxia Tolerance in Fishes
https://orcid.org/0000-0002-3870-6135
Babin Courtney H Department of Biological Sciences, University of New Orleans, New Orleans, LA 70148, USA

https://orcid.org/0000-0003-0249-9274
Leiva Félix P Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven 27570, Germany

https://orcid.org/0000-0002-0691-583X
Verberk Wilco C E P Department of Animal Ecology and Physiology, Radboud University Nijmegen, Nijmegen, The Netherlands

https://orcid.org/0000-0001-5636-1700
Rees Bernard B Department of Biological Sciences, University of New Orleans, New Orleans, LA 70148, USA

Bazykin Georgii Associate Editor
Corresponding author: E-mail: cchymel@uno.edu.
9 2024
21 8 2024
21 8 2024
16 9 evae18314 8 2024
03 9 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution.
2024
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Low dissolved oxygen (hypoxia) is recognized as a major threat to aquatic ecosystems worldwide. Because oxygen is paramount for the energy metabolism of animals, understanding the functional and genetic drivers of whole-animal hypoxia tolerance is critical to predicting the impacts of aquatic hypoxia. In this study, we investigate the molecular evolution of key genes involved in the detection of and response to hypoxia in ray-finned fishes: the prolyl hydroxylase domain (PHD)–hypoxia-inducible factor (HIF) oxygen-sensing system, also known as the EGLN (egg-laying nine)–HIF oxygen-sensing system. We searched fish genomes for HIFA and EGLN genes, discovered new paralogs from both gene families, and analyzed protein-coding sites under positive selection. The physicochemical properties of these positively selected amino acid sites were summarized using linear discriminants for each gene. We employed phylogenetic generalized least squares to assess the relationship between these linear discriminants for each HIFA and EGLN and hypoxia tolerance as reflected by the critical oxygen tension (Pcrit) of the corresponding species. Our results demonstrate that Pcrit in ray-finned fishes correlates with the physicochemical variation of positively selected sites in specific HIFA and EGLN genes. For HIF2A, two linear discriminants captured more than 90% of the physicochemical variation of these sites and explained between 20% and 39% of the variation in Pcrit. Thus, variation in HIF2A among fishes may contribute to their capacity to cope with aquatic hypoxia, similar to its proposed role in conferring tolerance to high-altitude hypoxia in certain lineages of terrestrial vertebrates.

hypoxia
critical oxygen tension
Actinopterygii
positive selection
hypoxia-inducible factor alpha
prolyl hydroxylase domain
==== Body
pmcSignificance

Decreased oxygen availability (hypoxia) poses a threat to animal life, especially in aquatic habitats where human activities have increased the duration, severity, and geographic extent of hypoxia. This surge in aquatic hypoxia is expected to unevenly impact the distribution of aquatic species, including ray-finned fishes, the most widespread and speciose group of vertebrates. This study investigates whether the hypoxia tolerance of ray-finned fishes, defined by the oxygen partial pressure threshold for sustaining aerobic metabolism, is linked to sequence variation in critical oxygen-sensing genes. The findings reveal that physicochemical variation in a key transcription factor is associated with variation in this measure of hypoxia tolerance of adult fish among the ray-finned fishes.

Introduction

Oxygen is essential for the development, growth, reproduction, and survival of contemporary life on Earth. Hence, a decrease in oxygen levels below fully air-saturated conditions (hypoxia) can be stressful, especially for water-breathing animals due to the much lower concentration and diffusion rates of oxygen in water compared to air (Nikinmaa and Rees 2005; Verberk et al. 2011). Although hypoxia occurs naturally in some aquatic habitats (Diaz and Breitburg 2009; Mandic and Regan 2018; Woods et al. 2022), the last several decades have seen a dramatic increase in the severity, duration, and geographic scope of aquatic hypoxia due to human activities (Breitburg et al. 2018; Sampaio et al. 2021; Xenopoulos et al. 2021). Higher water temperatures due to global warming decrease oxygen solubility, increase biological oxygen consumption rates, and intensify vertical stratification of the water column, thereby limiting the mixing of well-aerated surface waters with deeper oxygen-poor water. Nutrient runoff due to altered land use further contributes to oxygen depletion by stimulating the growth of nutrient-limited microorganisms. Thus, aquatic deoxygenation is recognized as an increasing threat to marine and freshwater ecosystems worldwide (Diaz and Rosenberg 2008; Breitburg et al. 2018; Jane et al. 2021).

The central molecular pathway regulating the cellular responses of animals to hypoxia includes the hypoxia-inducible factor (HIF) family of transcription factors and the prolyl hydroxylase domain (PHD) enzymes (Fong and Takeda 2008; Kaelin and Ratcliffe 2008; Semenza 2009; Ivan and Kaelin 2017). HIF is a family of heterodimeric transcription factors, with the functional transcription factor consisting of an oxygen-sensitive alpha subunit (HIFA) and a constitutively expressed beta subunit, also known as the aryl hydrocarbon receptor nuclear translocator (ARNT). The PHD enzymes catalyze the oxygen-dependent hydroxylation of specific proline residues of HIFA subunits, a modification that targets these subunits for ubiquitination and proteasomal degradation (Maxwell et al. 1999; Bruick and McKnight 2001; Epstein et al. 2001; Ivan et al. 2001; Jaakkola et al. 2001). The genes encoding PHD are also known as EGLN due to their homology with the egg-laying nine gene of the nematode Caenorhabditis elegans (Epstein et al. 2001). A drop in cellular levels of oxygen reduces rates of HIFA hydroxylation and degradation, allowing HIFA subunits to accumulate, bind to ARNT, and regulate the expression of hundreds of genes, many of which serve to improve oxygen transport to tissues or increase the capacity of cells to survive low oxygen levels (Wenger et al. 2005; Kaelin and Ratcliffe 2008; Semenza 2011, 2012).

Vertebrate animals have multiple copies of HIFA and EGLN genes (paralogs) arising from two rounds of genome duplication at the base of vertebrate evolution (Ohno 1970; Dehal and Boore 2005; Sacerdot et al. 2018). Most terrestrial vertebrates have three paralogs of HIFA (HIF1A, HIF2A, and HIF3A), which differ in their tissue distribution and target-gene specificity (Patel and Simon 2008; Keith et al. 2012; Duan 2016; Graham and Presnell 2017), and three paralogs of EGLN (EGLN1, EGLN2, and EGLN3), which vary in their affinity for HIFA subunits (Fong and Takeda 2008; Ivan and Kaelin 2017). The diversity of HIFA and EGLN is greater among ray-finned fishes (Actinopterygii) because of additional genome duplication events during their evolution, including the teleost-specific genome duplication (TGD, 350 to 320 Mya; Postlethwait et al. 2000; Hoegg et al. 2004; Volff 2005; Amores et al. 2011), and additional rounds of genome duplication in lineages leading to carp (19.5 Mya; Zhang et al. 2019; Wang et al. 2022), goldfish (12 to 10 Mya; Farhat et al. 2022; Wang et al. 2022), and salmonids (100 to 25 Mya; Berthelot et al. 2014; Macqueen and Johnston 2014). Phylogenetic analyses of HIFA in ray-finned fishes showed that many lineages retain the four HIFA genes predicted from two rounds of genome duplication in the ancestor to all vertebrates, with additional teleost-specific paralogs of HIF1A and HIF2A retained in the selected lineages (Rytkonen et al. 2013; Mandic et al. 2021; Townley et al. 2022). Moreover, evidence has been presented for additional HIFA paralogs arising from the genome duplication in lineages leading to carp (Zhang et al. 2019) and salmonids (Townley et al. 2022). Rytkonen et al. (2011) suggested that fishes have three EGLN genes that are homologous to those seen in terrestrial vertebrates, and more recently, Farhat et al. (2022) provided evidence of lineage-specific duplicates in goldfish.

Given this tremendous diversity in HIFA and EGLN among fishes, combined with the fact that Actinopterygii are the most speciose group of vertebrates and occur in aquatic habitats spanning a range of oxygen concentrations (Nelson et al. 2016), we hypothesized that sequence variation in HIFA and EGLN is related to variation in hypoxia tolerance among fishes. One measure of hypoxia tolerance is the critical oxygen tension (Pcrit), the lowest level of ambient oxygen at which the energetic costs of maintenance can be met by aerobic metabolism (Ultsch et al. 1981; Farrell and Richards 2009; Claireaux and Chabot 2016; Reemeyer and Rees 2019). At oxygen tensions below Pcrit, anaerobic processes are recruited or metabolism is suppressed. Because neither strategy can be maintained indefinitely, survival ultimately decreases at oxygen levels below Pcrit (Farrell and Richards 2009). It should be noted that the capacity to tolerate hypoxia is a complex phenotype, the utility of Pcrit has been questioned (Wood 2018), and alternative metrics of hypoxia tolerance of fishes exist (Alexander and McMahon 2004; Rees and Matute 2018; Seibel et al 2021). Nevertheless, Pcrit remains an ecologically and physiologically relevant index of the capacity for fish to extract oxygen from their surroundings (Regan et al. 2019), and it has been determined for a broader range of fishes than other measures of hypoxia tolerance (Rogers et al. 2016; Verberk et al. 2022a).

In the present study, therefore, we investigated whether genomic variation in the coding sequences (CDSs) of HIFA and EGLN is related to variation in the hypoxia tolerance of ray-finned fishes, as reflected by their Pcrit. Phylogenetic relationships were reconstructed for all HIFA and EGLN paralogs from currently available Actinopterygian genomes. Reasoning that positive selection could contribute to divergence within and among paralogs (Liao and Zhang 2006; Lian et al. 2020), we tested each gene for positive selection and grouped the HIFA and EGLN genes from different species according to the physicochemical properties of amino acid sites putatively under positive selection. While this analytical framework is similar to the recent analyses by Townley et al. (2022), the current study used more species, including 16 not represented in Townley et al. (2022), and a larger collection of paralogs, with over 200 new or updated sequences. The central and novel aspect of the current study, however, is that we took advantage of a comprehensive comparative analysis of Pcrit among fishes (Verberk et al. 2022a) and assessed the relationships between the physicochemical properties of HIFA and EGLN and the Pcrit values for the corresponding species. Prior to correcting for phylogeny, we found that physicochemical variation in HIFA and EGLN was significantly related to variation in Pcrit. After correcting for phylogeny, sequence variation in HIF1A and HIF2A was associated with variation in Pcrit, with variation in HIF2A being more strongly related to variation in Pcrit than other oxygen-sensing genes in this sample of ray-finned fishes.

Results and Discussion

Species Used for Analysis

This study included all species of ray-finned fishes for which sequenced genomes and critical oxygen tension (Pcrit) were available (Fig. 1). This collection of 28 species represents 14 orders of Actinopterygii spanning over 300 million years of evolution. This sample of ray-finned fishes included the spotted gar (Lepisosteus oculatus) that diverged approximately 320 Mya from other ray-finned fishes, i.e. prior to the TGD (Hoegg et al. 2004; Amores et al. 2011; Braasch et al. 2016; but see Davesne et al. 2021). Other species arose after the TGD and represent most major lineages of ray-finned fishes, including some that experienced additional lineage-specific genome or gene duplication events (Otocephala and Salmoniformes). These species exhibit a range of tolerance to reductions in ambient oxygen from very tolerant (low Pcrit) to very sensitive (high Pcrit).

Fig. 1. The phylogenetic relationships and hypoxia tolerance of the 28 species of Actinopterygii included in this study. The phylogeny was synthesized with the Open Tree of Life (opentreeoflife.org). Filled circles indicate genome duplication events at the base of vertebrate evolution (black) and during the evolution of teleost fishes (blue), carp (green), goldfish (yellow), and salmonids (orange). Species are color coded according to their phylogenetic placement: L. oculatus (spotted gar), representing a basal ray-finned fish that diverged prior to the TGD, black; the basal teleost A. anguilla (European eel), brown; Otocephala (Atlantic herring, tambaqui, zebrafish, grass carp, common carp, and goldfish), green; Salmonidae and their sister group (E. lucius, northern pike), orange; and Neoteleostei (more-derived teleosts), blue. On the right are values of standardized critical oxygen tensions (Pcrit in kPa), represented as dots (median) and bars (range) of values determined at three temperatures (15, 24, and 28 °C) for each species (see Materials and Methods). The vertical lines represent the lower (left) and upper (right) 20th percentile of the Pcrit values based on all ray-finned fish species in Verberk et al. (2022a) plus spotted gar, Atlantic herring, and northern pike. For comparison, the oxygen tension of air-saturated water is 21 kPa.

Phylogenetic Relationships of Actinopterygian HIFA and EGLN

Searching the genomes of these species resulted in 150 HIFA homologs and 109 EGLN homologs (supplementary table S1, Supplementary Material online). Phylogenetic analyses resolved four HIFA (Fig. 2) and three EGLN (Fig. 3) homology groups, largely congruent with earlier studies (Rytkonen et al. 2011, 2013; Li et al. 2022; Townley et al. 2022; Yu et al. 2023). Across all genes, the placement of species was essentially identical to the fossil-calibrated phylogeny of ray-finned fishes (Hughes et al. 2018). Importantly, our search did not recover HIF4A in Neoteleostei (Fig. 2), nor a fourth EGLN in any ray-finned fish (Fig. 3), suggesting that these were lost during (HIF4A) or prior to (EGLN4) the evolution of ray-finned fishes.

Fig. 2. ML phylogenetic inferences of Actinopterygian HIFA genes using CDSs. Evolutionary analyses were conducted in RAxML v 8.2.11 using the general time reversible substitution matrix (GAMMA+P-Invar model). The most likely tree inferred from 100 bootstrap replicates is shown for each gene with bootstrap values above the corresponding branches. Taxa are color coded as in Fig. 1. Where applicable, teleost-specific gene duplications are identified as “a” and “b” next to brackets. Genome duplications are indicated by filled circles: common carp-specific (green), goldfish-specific (yellow), and salmonid-specific (orange). Sequences are identified by the species names followed by the last 4 digits of the NCBI or Ensembl reference gene accession numbers (see supplementary table S1, Supplementary Material online, for a full list of genes).

Fig. 3. ML phylogenetic inferences of Actinopterygian EGLN genes using CDSs. Evolutionary analyses were conducted in RAxML v 8.2.11 using the general time reversible substitution matrix (GAMMA+P-Invar model). The most likely tree inferred from 100 bootstrap replicates is shown for each gene with bootstrap values above the corresponding branches. Taxa are color coded as in Fig. 1. Where applicable, teleost-specific gene duplications are identified as “a” and “b” next to brackets. Genome duplications are indicated by filled circles: common carp-specific (green), goldfish-specific (yellow), and salmonid-specific (orange). Sequences are identified by the species names followed by the last 4 digits of the NCBI or Ensembl reference gene accession numbers (see supplementary table S1, Supplementary Material online, for a full list of genes).

Spotted gar possesses only one paralog each of HIF1A-4A and EGLN1-3 as expected from its early divergence from the rest of ray-finned fishes. The remaining species investigated here occur in lineages that diverged after the TGD, resulting in teleost-specific paralogs, “a” and “b,” of each gene; however, two teleost-specific paralogs were only recovered for HIF1A, HIF2A, and EGLN1, and these were retained in only certain lineages. For HIF1A (Fig. 2), all searched species retained HIF1Aa, inferred to be more closely related to ancestral HIF1A by similarity of flanking genes (Gasanov et al. 2021; Townley et al. 2022), whereas only European eel (Anguilla anguilla) and Otocephala (including herring, tambaqui, zebrafish, carp, and goldfish) retained HIF1Ab. These results suggest that HIF1Ab was lost in the lineage leading to Neoteleostei after the divergence of Otocephala (228 Mya; Benton et al. 2015). On the other hand, most species have both teleost-specific paralogs of HIF2A. One form encodes a full-length protein in all species included here (supplementary table S1, Supplementary Material online). This paralog was the first HIF2A described in fishes (Powell and Hahn 2002) and shares more flanking genes with the ancestral HIF2A, represented by spotted gar (Townley et al. 2022). For these reasons, it is referred to as HIF2Aa (Rytkonen et al. 2013; Townley et al. 2022), although it has been designated as EPAS1b or HIF2Ab in many databases (supplementary table S1, Supplementary Material online). In contrast, the other paralog, herein referred to as HIF2Ab, encodes a full-length protein only in European eel and Otocephala and a truncated protein in Salmonidae (and their sister taxa Esox lucius) and Neoteleostei (supplementary table S1, Supplementary Material online; Rytkonen et al. 2013; Townley et al. 2022). EGLN1a was recovered in all taxa, but EGLN1b was present only in less-derived species (A. anguilla, Otocephala, E. lucius, Salmonidae) and absent in the Neotelestei (Fig. 3), consistent with the loss of this paralog after the divergence of Salmonidae (100 to 70 Mya; Shedko et al. 2013). The current analysis did not recover two paralogs of HIF3A, HIF4A, EGLN2, and EGLN3 in any species, suggesting that one form of each was nonfunctionalized shortly after the TGD, a common fate for gene duplicates (Lynch and Conery 2000).

Additional HIFA and EGLN paralogs were recovered in carp, goldfish, and salmonids, presumably resulting from lineage-specific genome or gene duplication events (Berthelot et al. 2014; Macqueen and Johnston 2014; Kuang et al. 2016; Zhang et al. 2019; Farhat et al. 2022; Wang et al. 2022). Carp-specific paralogs were recovered for HIF1Aa, HIF1Ab, HIF2Aa, HIF2Ab, and HIF3A, with additional goldfish-specific paralogs for HIF1Aa and HIF1Ab (green and yellow nodes, respectively, in Fig. 2), resulting in a minimum of 11 copies of various HIFA genes in common carp (Cyprinus carpio) and 14 copies in goldfish (Carassius auratus) (supplementary table S1, Supplementary Material online). Salmonid-specific duplicates of HIF1Aa, HIF2Aa, and HIF3A (orange nodes in Fig. 2) resulted in a total of eight copies of the various HIFA genes when including HIF2Ab and HIF4A (supplementary table S1, Supplementary Material online). Within the EGLN homology group, carp-specific duplicates were recovered for EGLN1a, EGLN1b, EGLN2, and EGLN3, and salmonid-specific duplicates were recovered for all EGLN genes except EGLN1a (green and orange nodes, respectively, in Fig. 3; supplementary table S1, Supplementary Material online). Many of these lineage-specific HIFA and EGLN paralogs have not been reported previously, and the elucidation of their expression and functions will require experiments capable of differentiating them (e.g. with paralog-specific probes).

Evidence of Positive Selection in the Evolution of Actinopterygian HIFA and EGLN

Reasoning that positive selection may have contributed to the diversity of HIFA and EGLN among fishes, we tested for positive selection using codon-based mixed effects models for episodic selection (mixed effects model of evolution [MEME]; Murrell et al. 2012) and gene-wide Bayesian tests of episodic diversification (branch-site unrestricted statistical test for episodic diversification [BUSTED]; Murrell et al. 2015). MEME detected episodic positive selection in all homology groups, with the number of sites putatively experiencing episodic positive selection ranging from eight for EGLN1 to 53 for HIF2A (Table 1). Most of these sites (>75%) also had BUSTED evidence ratios greater than 2, providing additional support for episodic positive selection occurring at these sites. Pervasive positive selection was evaluated with fixed effects likelihoods (FEL) for pervasive selection (Kosakovsky Pond and Frost 2005b) and found to occur much less frequently than episodic selection in all HIFA and not at all for EGLN genes (Table 1). This is consistent with the notion that episodic positive selection occurs at more sites than pervasive selection across a wide range of organisms (Murrell et al. 2012), potentially due to greater selective pressures in those lineages that experience more variable environments.

Table 1 Codons putatively under positive selection in HIFA and EGLN genes in ray-finned fishes

Gene	Episodic selection	Pervasive selection	
HIF1A	97, 140, 142, 183, 274, 295, 315, 354, 469, 483, 524, 629, 630, 642, 656, 660, 735, 805, 810, 872, 883, 894, 897, 905, 921, 923, 956, 976	543	
HIF2A	109, 156a, 199, 200, 201, 204, 268a, 281, 285, 288, 545, 547, 611, 613, 634, 651, 681a, 683, 684, 686, 707, 708, 712, 715, 717, 722, 724, 727, 731, 750, 761, 773, 786, 800, 816, 887, 888a, 912, 929, 953, 1009, 1079, 1085, 1088, 1090, 1107, 1118a, 1121, 1163, 1165, 1166, 1167, 1234	156a, 240, 268a, 681a, 888a, 1118a	
HIF3A	11 a , 25, 89, 176, 311, 456, 461, 462, 483a, 485, 491, 516a, 525, 577, 597, 599a, 602, 610, 643, 658, 670, 675, 688, 720, 832	11 a , 473, 483a, 516a, 599a	
HIF4A	71, 160, 303, 310, 394a, 546, 561, 715a, 779a, 807, 854, 867, 879, 892a	75, 394a, 715a, 779a, 892a	
EGLN1	13, 78, 79, 89, 90, 91, 167, 395	NAb	
EGLN2	21, 45, 107, 245, 289, 314, 316, 335, 379, 442, 489, 281, 593, 684, 687, 693, 695, 696, 697, 701	NAb	
EGLN3	29, 211, 223, 336, 376, 420, 461, 469, 476, 477	NAb	
Episodic positive selection was detected by MEME, and pervasive positive selection was detected by FEL. Codon numbers correspond to the positions in the multiple sequence alignment of all HIFA and EGLN genes used in each analysis. Codons with episodic positive selection supported by BUSTED evidence ratio > 2 are shown in bold type. Codons found to be under positive selection in Townley et al. (2022) are shown in italics. See supplementary tables S2 to S8, Supplementary Material online, for amino acid identities for each HIFA and EGLN gene.

aPositively selected codons detected by both MEME and FEL. bNA, not applicable: no sites detected by FEL.

Two observations can be made about the number of sites experiencing positive selection in HIFA and EGLN (Table 1). First, the number of sites is larger for the HIFA genes compared to EGLN genes, even after accounting for the longer CDSs for HIFA genes. Second, the number of codons under positive selection expressed as a proportion of the total number of codons was significantly greater in HIF2A than other HIFA genes (P < 0.05; χ2 test). Both observations could reflect a larger role for positive selection in the evolution in HIFA, HIF2A in particular, than for EGLN genes. However, we cannot exclude the possibility that those larger numbers of sites have experienced weaker positive selection, as compared to stronger positive selection acting on fewer sites in other genes.

Of note, many of the positively selected sites found here occurred in conserved structural or functional domains. Sites under positive selection in each HIFA homology group occurred in domains involved in DNA binding, protein dimerization, oxygen-dependent degradation, or activation of gene expression (supplementary tables S2 to S5, Supplementary Material online), while several sites under positive selection in the EGLN genes were mapped to the active site (supplementary tables S6 to S8, Supplementary Material online). In addition, for HIF2A, some of the positively selected sites aligned with sites known to be subject to posttranslational modification in mammals (supplementary table S3, Supplementary Material online; Albanese et al. 2020; Daly et al. 2021).

Using similar approaches, Townley et al. (2022) recently measured the prevalence of positive selection on a smaller set of fish HIFA genes, including some paralogs from different species. The current results agree with Townley et al. (2022) in two important ways: HIF2A was found to have more positively selected sites than other HIFA homologs and positive selection in all HIFA genes was found at sites predicted to be critical for protein structure and function. The specific codons identified in these two studies, however, were largely different. For HIF1A, HIF2A, and HIF3A, only five codons from each gene were shared between the current study and Townley et al. (2022), and for HIF4A, only one site was common (Table 1, italicized codons). Tests of positive selection are known to be sensitive to the number, quality, and alignment of sequences (Wong et al. 2008; Murrell et al. 2012) indicating that caution should be exercised when ascribing functional importance to any specific site identified by these tests.

Nevertheless, the current results and those of Townley et al. (2022) complement studies reporting signatures of positive selection in the EGLN–HIF pathway during adaptation to hypoxia associated with high altitude. Targets of positive selection include HIF2A and EGLN1 in Tibetan human populations (Beall et al. 2010; Eichstaedt et al. 2017; Pamenter et al. 2020); HIF2A in the Tibetan mastiff (Li et al. 2014); HIF2A in high-altitude deer mice (Schweizer et al. 2019); and HIF1A and HIF2A in fishes from the Tibetan plateau (Guan et al. 2014; Wang et al. 2015a, b). Thus, it appears that positive selection on HIFA and EGLN has occurred among aquatic and terrestrial vertebrate lineages living in high-altitude habitats characterized by reductions in oxygen availability, potentially contributing to the hypoxia tolerance of certain species.

HIFA and EGLN Homologs Group According to Properties of Positively Selected Sites

To better understand the potential physicochemical consequences of this positive selection, we distinguished groups of HIFA and EGLN homologs by the characteristics of amino acids putatively under positive selection using discriminant analysis of principal components (DAPC). This is a novel application of a method often employed to evaluate genetic diversity and population structure based on single nucleotide polymorphisms (SNPs) and microsatellites (Jombart et al. 2010; Dotti do Prado et al. 2017; Kajungiro et al. 2019; Sunde et al. 2020). In brief, variation in the physicochemical properties among amino acids occurring at each positively selected site for a given HIFA or EGLN homology group was summarized by principal components, followed by clustering groups of homologs (DAPC groups) along orthogonal axes of variation, or linear discriminants (LDs). Depending upon the gene, two to four DAPC groups were distinguished by their scores along one to three LDs (Table 2).

Table 2 Number of gene groups identified by DAPC of the physicochemical properties of positively selected sites in HIFA and EGLN in ray-finned fishes

Homolog	DAPC groups	LD1	LD2	LD3	
HIF1A	4	0.62	0.24	0.14	
HIF2A	4	0.80	0.14	0.06	
HIF3A	3	0.57	0.43	NA	
HIF4A	2	1.00	NA	NA	
EGLN1	4	0.56	0.42	0.02	
EGLN2	2	0.77	0.23	NA	
EGLN3	2	1.00	NA	NA	
DAPC groups were distinguished by their scores on one to three LDs, and the proportion of physicochemical variation explained by each LD is shown. See supplementary tables S2 to S8, Supplementary Material online, for specific genes in each DAPC group and codons that heavily weighted on each LD.

The species and homologs making up the DAPC groups differed for each gene (supplementary figs. S1 to S7, Supplementary Material online). As expected, phylogeny was a main contributor to the homolog groupings, e.g. segregating homologs from Otocephala, Salmonidae, and Neoteleostei for HIF1A (supplementary fig. S1, Supplementary Material online), HIF3A (supplementary fig. S3, Supplementary Material online), EGLN1 (supplementary fig. S5, Supplementary Material online), and EGLN2 (supplementary fig. S6, Supplementary Material online). Among Otocephala, the cyprinids (zebrafish, Danio rerio; grass carp, Ctenopharyngodon idella; common carp, Cyprinus carpio; goldfish, Carassius auratus) frequently grouped together and apart from other Otocephala (Atlantic herring, Clupea harengus; tambaqui, Colossoma macropomum), which grouped with various taxa depending upon the gene analyzed. This analysis also discriminated between teleost-specific paralogs of HIF2A and EGLN1 (supplementary figs. S2 and S5, Supplementary Material online, respectively). The codons that loaded most heavily on each LD for each gene were determined (supplementary figs. S1 to S7, Supplementary Material online), along with the amino acid residues at these sites in all genes (supplementary tables S2 to S8, Supplementary Material online).

Relating Physicochemical Variation of HIFA and EGLN to Hypoxia Tolerance

While Townley et al. (2022) used similar approaches to distinguish among a different set of HIFA paralogs in fishes, the current study is novel in determining if physicochemical variation of either HIFA or EGLN is related to the hypoxia tolerance of fishes. We used phylogenetic generalized least squares (PGLS) to assess the relationship between the LD scores of each HIFA and EGLN homolog and the standardized Pcrit values from the corresponding species. These analyses were performed for standardized Pcrit values determined at three temperatures, 15, 24, and 28 °C. Moreover, these analyses were done excluding the influence of phylogeny (λ = 0.001), as well as after accounting for the relationships among the paralogs in a given analysis (using model-selected values of λ).

Tables 3 and 4 present the best PGLS models describing variation in Pcrit at 24 °C excluding and including the influence of phylogeny, respectively. Without accounting for phylogeny, variation in Pcrit was found to be associated with physicochemical variation for all genes except HIF4A and EGLN3 as shown by the retention of one or more LD in the best models (Table 3). Although the goodness of fit (R2) and the significance of individual LDs (P values) differed somewhat for Pcrit at 15 °C (supplementary table S9, Supplementary Material online) and Pcrit at 28 °C (supplementary table S10, Supplementary Material online), the general results were largely congruent across temperatures. When phylogeny was taken into account, Pcrit was found to be associated with physicochemical variation in HIF1A and HIF2A (Table 4), which, again, was consistent across temperatures (supplementary table S11, Supplementary Material online, for Pcrit at 15 °C; supplementary table S12, Supplementary Material online, for Pcrit at 28 °C). After accounting for phylogeny, PGLS models relating variation in HIF3A, EGLN1, or EGLN2 to Pcrit were not significant at any temperature.

Table 3 Relationships between the physicochemical properties of positively selected amino acid sites in Actinopterygian HIFA and EGLN and critical oxygen tension (Pcrit) at 24 °C, without accounting for phylogeny

Effect	Estimate	SE	t	P	
HIF1A (λ = 0.001; R2 = 0.480)					
 (Intercept)	5.581	0.407	13.727	<0.001	
 LD2	0.084	0.091	0.917	0.364	
 LD3	0.635	0.113	5.632	<0.001	
HIF2A (λ = 0.001; R2 = 0.507)					
 (Intercept)	5.952	0.462	7.603	<0.001	
 LD1	−0.158	0.033	−4.002	<0.001	
 LD2	−0.411	0.077	−4.836	<0.001	
HIF3A (λ = 0.001; R2 = 0.238)					
 (Intercept)	6.496	0.619	10.486	<0.001	
 LD1	−0.260	0.175	−1.484	0.149	
 LD2	0.513	0.210	2.444	0.021	
EGLN1 (λ = 0.001; R2 = 0.340)					
 (Intercept)	6.004	0.491	12.235	<0.001	
 LD1	0.479	0.118	4.057	<0.001	
 LD2	0.318	0.163	1.945	0.059	
EGLN2 (λ = 0.001; R2 = 0.121)					
 (Intercept)	6.691	0.619	10.810	<0.001	
 LD2	−0.419	0.203	−2.065	0.047	
The effects of LD scores from DAPC analyses of each gene on Pcrit were tested with PGLS with lambda = 0.001 (no phylogenetic influence) after model reduction by ANOVA (see Materials and Methods). Lambda values (λ) and R2 values for each model are given in parentheses.

Table 4 Relationships between the physicochemical properties of positively selected amino acid sites in Actinopterygian HIFA and EGLN and critical oxygen tension (Pcrit) at 24 °C, after accounting for phylogeny

Effect	Estimate	SE	t	P	
HIF1A (λ = 0.000; R2 = 0.481)					
 (Intercept)	5.577	0.403	13.850	<0.001	
 LD2	0.084	0.091	0.919	0.364	
 LD3	0.635	0.113	5.642	<0.001	
HIF2A (λ = 0.534; R2 = 0.346)					
 (Intercept)	8.130	1.083	7.505	<0.001	
 LD1	−0.138	0.073	−1.903	0.064	
 LD2	−0.460	0.122	−3.776	<0.001	
The effects of LD scores from DAPC analyses of each gene on Pcrit were tested with PGLS using the model-selected values of lambda (to account for phylogeny) after model reduction by ANOVA (see Materials and Methods). Lambda values (λ) and R2 values for each model are given in parentheses.

We expected that including phylogeny would reduce the power of sequence variation to explain variation in Pcrit because both the gene sequences and Pcrit (Verberk et al. 2022a) have strong phylogenetic signals. This is not to say that phylogenetically related variation in the physicochemical properties of HIF3A, EGLN1, and EGLN2 is unimportant. Rather, it is simply not possible to distinguish the role of adaptation versus shared evolutionary history in shaping this variation. The more striking result was that, even after including the influence of phylogeny, physicochemical variation in positively selected sites of HIF1A and HIF2A was significantly associated with variation in Pcrit across temperatures.

The relationships between the standardized Pcrit at 24 °C and each LD retained in the phylogenetically informed PGLS models (P < 0.1) are shown in Fig. 4 (supplementary fig. S8, Supplementary Material online, for Pcrit at 15 °C; supplementary fig. S9, Supplementary Material online, for Pcrit at 28 °C). Variation in HIF1A LD3 tended to discriminate between Otocephala, characterized by low Pcrit (high hypoxia tolerance; green symbols), and Salmonidae, characterized by high Pcrit (low hypoxia tolerance; orange symbols) (Fig. 4a; supplementary figs. S8a and S9a, Supplementary Material online). Variation in HIF2A LD1 separated Otocephala HIF2Ab, specifically those from cyprinids, from HIF2Aa in all other fishes (Fig. 4b; supplementary figs. S8b and S9b, Supplementary Material online). Variation in HIF2A LD2 tended to segregate salmonid HIF2Aa from all others (Fig. 4c; supplementary figs. S8c and S9c, Supplementary Material online). While all of these relationships are significant, or nearly so (e.g. P = 0.064 for HIF2A LD1 at 24 °C), it is important to note that these LDs explain vastly different amounts of the overall physicochemical variation in the respective genes (Table 2). HIF2A LD1 explains ∼80% of the variation in HIF2A, whereas HIF2A LD2 and HIF1A LD3 each explain 14% of the physicochemical variation in those genes.

Fig. 4. Relationships between physicochemical variation in Actinopterygian HIF1A and HIF2A and critical oxygen tension (Pcrit) at 24 °C. Standardized Pcrit (see Materials and Methods) are plotted against LD scores from DAPC analysis of a) HIF1A or b) and c) HIF2A. Only LDs having P < 0.10 in the final PGLS models are plotted (see Table 4). The regression lines are from PGLS, including the effects of phylogeny. Symbol colors are basal ray-finned fish (spotted gar), black; basal teleost (European eel), brown; Otocephala, green; Salmonidae and northern pike, orange; and Neoteleostei, blue. Symbol shapes are ancestor of teleost-specific duplicates, cross; teleost-specific “a” paralog, circles; and teleost-specific “b” paralog, triangles.

We carried out a separate analysis using randomly selected amino acid sites to test whether the patterns observed above were due to chance. For this analysis, we selected 52 amino acid sites from HIF2A that were not shown to be under positive selection (compared to 53 amino acids found to be under positive selection for this gene, Table 1). DAPC returned four DAPC groups distinguished from one another by three LDs (the same number of groups and LDs as for positively selected sites from HIF2A). PGLS without considering phylogeny (λ = 0.001) showed that LD2 was related to standardized Pcrit at all three temperatures (P < 0.1; supplementary table S13, Supplementary Material online), although the PGLS models themselves had low R2 and were not significant. When phylogeny was included in the analysis (model-selected λ), no individual LDs were significantly related to Pcrit at any temperature (supplementary table S14, Supplementary Material online), and correspondingly, no models were significant. These results suggest that the physicochemical properties of randomly selected amino acid sites can be weakly related to variation in Pcrit, but that this variation is completely explained by phylogeny. Therefore, the relationships between physicochemical variation at positively selected sites in HIF1A and HIF2A and Pcrit that were found with PGLS after accounting for phylogeny (Table 4) cannot be explained by random chance.

Variation in Pcrit Is Best Explained by Physicochemical Variation in HIF2A

PGLS models including HIF2A LD1 and LD2 explained between 20% and 39% of the variation in Pcrit, after accounting for phylogeny (model R2 from Table 4 and supplementary tables S11 and S12, Supplementary Material online). In addition, LD1 and LD2 together accounted for 94% of the physicochemical variation in amino acid sites putatively under positive selection in this sample of ray-finned fish HIF2A (Table 2). These results strongly suggest that sequence variation in HIF2A is associated with hypoxia tolerance of fishes. As mentioned above, positive selection in HIF2A has been proposed to contribute to the adaptation of various vertebrate groups to high-altitude hypoxia (Beall et al. 2010; Guan et al. 2014; Li et al. 2014; Wang et al. 2015a, b; Eichstaedt et al. 2017; Schweizer et al. 2019; Pamenter et al. 2020). In Tibetan human populations, sequence variation in HIF2A is associated with a diminished increase in hematocrit at altitude and lower incidence of chronic mountain sickness (Beall et al. 2010; Basang et al. 2015). In North American deer mice (Peromyscus maniculatus), a high-altitude variant of HIF2A (Epas1) is associated with higher heart rate under hypoxia and potentially blunted adrenal catecholamine release (Schweizer et al. 2019).

In the current study, the first major axis of variation (LD1) distinguished HIF2Ab from the Cyprinidae, having high hypoxia tolerance, from all other HIF2A, and the second major axis (LD2) separated HIF2Aa from Salmonidae, having low hypoxia tolerance, from other HIF2A. Recently, Zhao et al. (2023) showed that SNP variation in HIF2Ab is correlated with hypoxia tolerance in the blunt nose bream (Megalobrama amblycephala), a member of the Cyprinidae. In ray-finned fishes, HIF2A is primarily expressed in gill tissue (Townley et al. 2022), where it may be involved in oxygen sensing (Pan et al. 2022). Furthermore, carp and goldfish exhibit gill remodeling during hypoxia or elevated temperature conditions, resulting in greater surface area available for gas exchange (Sollid et al. 2003; Sollid et al. 2005). Whether the variation we document in HIF2A contributes to hypoxia tolerance by altering the development, physiology, or morphology of gills in ray-finned fishes remains to be explored.

In comparison to HIF2A, PGLS models relating variation in HIF1A to Pcrit, while explaining a similar amount of variation in Pcrit (Table 4; supplementary tables S11 and S12, Supplementary Material online), retained LD3 which only explained 14% of the physicochemical variation of the positively selected sites in HIF1A (Table 2). This result suggests that physicochemical variation in HIF1A is a poor predictor of hypoxia tolerance of fishes. This result is consistent with Rytkonen et al. (2007), who failed to demonstrate clear relationships between HIF1A amino acid variation and the oxygen sensitivity of a smaller sample of ray-finned fishes. In addition, Mandic et al. (2020) demonstrated that deletion of both paralogs of HIF1A (HIF1Aa/Ab) in zebrafish had no effect on Pcrit of adult fish. Other metrics of hypoxia tolerance, however, were affected by this double knockout: as larvae, knockout fish displayed higher Pcrit (less hypoxia tolerance) than wild-type larvae, and adult knockout fish became impaired (lost equilibrium) sooner than wild-type adults when exposed to water with oxygen tensions lower than Pcrit (Mandic et al. 2020). Recently, Elbassiouny et al. (2024) documented variation in HIF1Ab among knifefish species (order Gymnotiformes, also included in the Otocephala) that was correlated with oxygen levels of their natural habitats. Based upon in vitro assays, these authors proposed that knifefishes from habitats more prone to hypoxia have a variant of HIF1Ab capable of greater transactivation of gene expression during conditions of low oxygen. Thus, although our data provide limited support for a relationship between sequence variation in HIF1A and Pcrit of adults, there is little doubt that it contributes to other aspects of hypoxia tolerance in fishes.

Conclusions

The hypoxia tolerance of fishes is a complex organismal phenotype that is almost certainly controlled by several genes (Healy et al. 2018). Despite this complexity, we demonstrate that physicochemical variation in amino acid sites putatively under positive selection in a single gene, the transcription factor HIF2A, is strongly associated with variation in the Pcrit across a broad taxonomic range of fishes. The Pcrit of adult fishes is only one measure of how well fishes tolerate reductions in ambient oxygen, albeit the one most widely documented. We cannot exclude the possibility that the variation of other HIFA and EGLN homologs described here may be related to other metrics of hypoxia tolerance of fishes or linked to variation in other aspects of their life histories. Future studies should explore the functional links between sequence variation in HIFA and how fishes, and other aquatic organisms, deal with naturally occurring and human-induced episodes of low dissolved oxygen.

Materials and Methods

Gene Identification and Alignment

We used several search terms (including “hif”, “hypoxia”, “hypoxia inducible factor”, “epas”, “endothelial”, “endothelial PAS”, “egl”, “egln”, “prolyl”, and “prolyl hydroxylase”) to identify four HIFA and three EGLN genes from 28 Actinopterygian species with published genomes at NCBI (https://www.ncbi.nlm.nih.gov) or Ensembl (http://www.ensembl.org/) through March 2023 (supplementary table S1, Supplementary Material online). Full gene sequences were downloaded, and we used BLASTn and exploratory multiple sequence alignments to identify and group sequences with questionable gene descriptions (e.g. “hypoxia-inducible factor 1-like”). We excluded exact duplicate, low quality, or partial sequences from the curated data sets for all genes except for four instances, where partial sequences that had similar gene names or were on the same chromosome encoded sequential exons of a given gene. We suspected these were called partial sequences or were on unplaced scaffolds due to errors in sequencing or annotating. Thus, these sequences were concatenated and used for all analyses (see Gene ID in supplementary table S1, Supplementary Material online, for details). The longest corresponding CDSs were extracted from the full gene sequences, and nucleotide alignments for each gene were conducted in MACSE v2.06 (Ranwez et al. 2011, 2018) with corresponding amino acid alignments being produced with the BLOSUM62 score matrix (Henikoff and Henikoff 1992).

Phylogenetic Inference

Maximum likelihood (ML) gene trees were inferred for each gene (e.g. HIF1A or EGLN1) using rapid bootstrapping and subsequent ML search in RAxML v8.2.11 (Stamatakis 2014) as implemented through Geneious Prime v2023.1.2 (www.geneious.com). This was accomplished by drawing bootstrap support values on the best-scoring ML tree from 100 bootstrap inferences using the general time reversible substitution matrix (GAMMA+P-Invar model). All trees were visualized and edited in Geneious Prime v2023.1.2 (www.geneious.com).

Positive Selection Analyses

We created combined fasta files of the CDS nucleotide alignments and ML trees for each gene to perform selection analyses in the HyPhy package (Kosakovsky Pond et al. 2005; Kosakovsky Pond et al. 2020) through the Datamonkey webserver (Kosakovsky Pond and Frost 2005a; Delport et al. 2010; Weaver et al. 2018). Frameshifts and stop codons that may result from multiple sequence alignments are not allowed in the data used for selection analyses. Therefore, these characters were replaced with the exportAlignment program in MACSE v2.06 (Ranwez et al. 2011, 2018). For each gene and corresponding phylogeny, we used MEME (Murrell et al. 2012) to assess whether individual sites were subject to episodic selection on a proportion of branches and FEL (Kosakovsky Pond and Frost 2005b) to identify sites of pervasive selection. We then implemented BUSTED (Murrell et al. 2015), which evaluates whether a gene has experienced positive selection at any site on at least one branch given a phylogeny. Evidence ratios from BUSTED for sites identified to be under positive selection by MEME and/or FEL were used as additional support for positive selection. The evidence ratio is a log-likelihood ratio that, when >2, provides support for positive selection compared to the null model (Murrell et al. 2015).

MEME provides a robust quantitative improvement compared to traditional models that are unable to detect instances of episodic positive selection due to variable levels of purifying selection pressure across different lineages, resulting in qualitatively different conclusions (Murrell et al. 2012). While FEL evaluates whether some sites evolved primarily under significant purifying selection, MEME can detect the signature of positive selection on certain branches. Unlike traditional models, MEME does not assume constant selective forces across all lineages, which allows the strength and direction of natural selection to vary both from site to site and from branch to branch at a site. This flexibility enables MEME to capture widespread episodic selection, which may have been previously underestimated. However, MEME is limited in that it assumes independent selective pressures between branches, and this assumption might be violated if the ratio of nonsynonymous to synonymous substitutions (ω) changes slowly across a phylogeny (Murrell et al. 2012). Because MEME uses a mixed effects model that allows for different branches and sites to have different ω rates, it is sensitive to updated alignments when sequence data have changed, additional paralogs are included, or species are added to a data set, resulting in a change in ω estimates for some sites or their significance levels, depending on the model fit and phylogenetic signal. It should be noted that even traditional tests of positive selection are sensitive to alignment methods which can lead to different conclusions in detecting positively selected sites (Wong et al. 2008).

Physicochemical Similarity Assessment

For each HIFA and EGLN homology group, codons identified as sites under positive selection were scored by five z-descriptors of amino acid physicochemical properties as described by Sandberg et al. (1998). The descriptors included hydrophobicity (z1), steric bulk (z2), polarity (z3), and electronic effects (z4 and z5). When codons were missing due to gaps in the multiple sequence alignments, the column means for z-descriptors were used. Salmonidae and Neoteleostei retain truncated forms of HIF2Ab with CDSs less than half the length of “full-length” HIF2Ab (Rytkonen et al. 2013; Townley et al. 2022). Preliminary analyses showed that they formed their own physicochemical group (Townley et al. 2022) and they were removed from the analysis of other HIF2A because of the large number of missing sites compared to the longer genes. The adegenet package (Jombart 2008; Jombart and Ahmed 2011) in R v4.2.1 (R Development Core Team 2021) was used to summarize all descriptors from all sites by principal components analysis (PCA). DAPC was performed on the minimum number of retained PCs that explained approximately 90% of the variation. DAPC is a robust multivariate method that serves as a bridge between PCA and discriminant analysis (DA). By transforming data using PCA before applying DA, DAPC ensures that the variables under consideration are uncorrelated, which reduces the dimensionality of the data set and ensures that the number of variables is less than the number of analyzed individuals (a necessary condition for DA). DAPC uses the k-means algorithm to infer the minimum number of gene groups by clustering PCs to construct linear combinations of the input variables with the greatest variation between groups and the smallest variation within groups. This results in orthogonal LDs that are axes of variation in retained PCs. The amino acids whose physicochemical properties were heavily loaded on each LD (90th percentile) were mapped to associated codons in the multiple sequence alignments for each gene.

Critical Oxygen Tensions for Target Species

We employed a database of Pcrit values from marine, freshwater, and brackish fishes (Verberk et al. 2022a, b). This database provides multiple Pcrit values for each species and includes metrics describing methodological details associated with each Pcrit estimation, such as the test temperature and test salinity. Additionally, the database encompasses biological factors (e.g. body size, genome size, and metabolic rates) that collectively contribute to the observed variability in Pcrit (Verberk et al. 2022a). A data set of 171 species with a complete set of variables was supplemented with three species for which HIFA and EGLN sequences were available, but lacking Pcrit values: spotted gar (L. oculatus), Atlantic herring (C. harengus), and northern pike (E. lucius). The first of these represents a lineage of ray-finned fishes lacking the TGD (Braasch et al. 2016), while the latter two represent sister groups of Cyprinidae and Salmonidae, respectively. Pcrit data for these three species were obtained by searching the Web of Science Core Collection using similar keyword combinations as Verberk et al. (2022a). We refit the model of Verberk et al. (2022a) incorporating these three species to estimate standardized values for Pcrit at 15, 24, and 28 °C (corresponding to the 25th percentile, median, and 75th percentile of measurement temperatures included in the database) for the 28 species of ray-finned fishes for which genomes were available. The standardized Pcrit values accounted for the effects of test salinity, test temperature, metabolic rate, genome size, and body mass, as well as the interactions between temperature and genome size and temperature and body mass (Verberk et al. 2022a).

Relating Hypoxia Tolerance to HIFA and EGLN Physicochemical Variation

The effects of physicochemical variation at positively selected amino acids in HIFA and EGLN on standardized Pcrit values were assessed by PGLS models. The scores for each LD were extracted from the DAPC analysis for each member of a gene family (e.g. HIF1A) and regressed against standardized Pcrit values at a given temperature (15, 24, or 28 °C). Two PGLS models were conducted for each of the three temperatures using the best-scoring ML gene tree inferred for each gene: a model that does not account for phylogeny (λ = 0.001) and a model corrected for phylogeny by optimizing lambda using ML (λ = ML). Model reduction was performed with analysis of variance (ANOVA) to retain LDs that had P < 0.1, and then, PGLS models of the reduced ANOVAs were summarized. A post hoc analysis of 52 random amino acid sites from HIF2A was also conducted to compare with the performance and results of the positively selected sites used in the DAPC and PGLS analyses. All methods for this post hoc analysis were the same as those used for the experimental data sets with one exception: PGLS models were not reduced by ANOVA, which would have resulted in null models. Thus, we present results for the full models both without and with phylogenetic correction.

All analyses were conducted in R v4.1.0 (R Development Core Team 2021) by using the following packages: “ape v5.7-1” (Paradis and Schliep 2019), “geiger v2.0.11” (Pennell et al. 2014), “caper v1.0.1” (Orme et al. 2018), “tidyverse v2.0.0” (Wickham et al. 2019), and “phytools v1.5-1” (Revell 2012).

Supplementary Material

evae183_Supplementary_Data

Acknowledgments

We gratefully acknowledge Iris L. E. van de Pol for her contributions to this work. Funding provided by the Greater New Orleans Foundation was granted to B.B.R., the Alexander von Humboldt Fellowship for postdoctoral researchers was granted to F.P.L., and W.C.E.P.V. gratefully acknowledges funding from the Dutch Research Council (NWO-VIDI 016.161.321).

Supplementary Material

Supplementary material is available at Genome Biology and Evolution online.

Data Availability

Data files and code supporting the analyses, figures, and tables of this study are publicly available on GitHub (https://github.com/felixpleiva/Genetic_basis_Pcrit) and Zenodo (Babin et al. 2023).
==== Refs
Literature Cited

Albanese  A, Daly  LA, Mennerich  D, Kietzmann  T, See  V. The role of hypoxia-inducible factor post-translational modifications in regulating its localization, stability, and activity. Int J Mol Sci. 2020:22 (1 ):268. 10.3390/ijms22010268.33383924
Alexander  JE  Jr, McMahon  RF. Respiratory response to temperature and hypoxia in the zebra mussel Dreissena polymorpha. Comp Biochem Physiol Part A. 2004:137 (2 ):425–434. 10.1016/j.cbpb.2003.11.003.
Amores  A, Catchen  J, Ferrara  A, Fontenot  Q, Postlethwait  JH. Genome evolution and meiotic maps by massively parallel DNA sequencing: spotted gar, and outgroup for the teleost genome duplication. Genetics. 2011:188 (4 ):799–808. 10.1534/genetics.111.127324.21828280
Basang  Z, Wang  B, Li  L, Yang  L, Liu  L, Cui  C, Lanzi  G, Yuzhen  N, Duo  J, Zheng  H, et al  HIF2A variants were associated with different levels of high-altitude hypoxia among native Tibetans. PLoS One. 2015:10 (9 ):e0137956. 10.1371/journal.pone.0137956.26368009
Beall  CM, Cavalleri  GL, Deng  L, Zheng  YT. Natural selection on EPAS1 (HIF2α) associated with low hemoglobin concentration in Tibetan highlanders. Proc Natl Acad Sci U S A.  2010:107 (25 ):11456–11464. 10.1073/pnas.1002443107.
Benton  MJ, Donoghue  PCJ, Asher  RA, Friedman  M, Near  TJ, Vinther  J. Constraints on the timescale of animal evolutionary history. Palaeont Electr. 2015:18 (1 ):1–16. 10.26879/424.
Berthelot  C, Brunet  F, Chalopin  D, Juanchich  A, Bernard  M, Noel  B, Bento  P, Da Silva  C, Labadie  K, Alberti  A, et al  The rainbow trout genome provides novel insights into evolution after whole-genome duplication in vertebrates. Nat Commun. 2014:5 (1 ):3657. 10.1038/ncomms4657.24755649
Braasch  I, Gehrke  AR, Smith  JJ, Kawasaki  K, Manousaki  T, Pasquier  J, Amores  A, Desvignes  T, Batzel  P, Catchen  J, et al  The spotted gar genome illuminates vertebrate evolution and facilitates human-teleost comparisons. Nat Genet. 2016:48 (4 ):427–437. 10.1038/ng.3526.26950095
Breitburg  D, Levin  LA, Oschlies  A, Gregoire  M, Chavez  FP, Conley  DJ, Garcon  V, Gilbert  D, Gutierrez  D, Isensee  K, et al  Declining oxygen in the global ocean and coastal waters. Science. 2018:359 (6371 ):eaam7240. 10.1126/science.aam7240.29301986
Bruick  RK, McKnight  SL. A conserved family of prolyl-4-hydroxylases that modify HIF. Science. 2001:294 (5545 ):1337–1340. 10.1126/science.1066373.11598268
Claireaux  G, Chabot  D. Responses by fishes to environmental hypoxia: integration through Fry's concept of aerobic metabolic scope. J Fish Biol. 2016:88 (1 ):232–251. 10.1111/jfb.12833.26768976
Daly  LA, Brownridge  PJ, Batie  M, Rocha  S, See  V, Eyers  CE. Oxygen-dependent changes in binding partners and post-translational modifications regulate the abundance and activity of HIF-1α/2α. Sci Signal. 2021:14 (692 ):eabf6685. 10.1126/scisignal.abf6685.34285132
Davesne  D, Friedman  M, Schmitt  AD, Fernandez  V, Carnevale  G, Ahlberg  PE, Sanchez  S, Benson  RBJ. Fossilized cell structures identify an ancient origin for the teleost whole-genome duplication. Proc Natl Acad Sci U S A. 2021:118 (30 ):e2101780118. 10.1073/pnas.2101780118.34301898
Dehal  P, Boore  JL. Two rounds of whole genome duplication in the ancestral vertebrate. PLoS Biol. 2005:3 (10 ):e314. 10.1371/journal.pbio.0030314.16128622
Delport  W, Poon  AFY, Frost  SDW, Kosakovsky Pond  SL. Datamonkey 2010: a suite of phylogenetic analysis tools for evolutionary biology. Bioinform. 2010:26 (19 ):2455–2457. 10.1093/bioinformatics/btq429.
Diaz  RJ, Breitburg  DL. The hypoxic environment. In: Richards  JG, Farrell  AP, Brauner  CJ, editors. Fish physiology. Vol. 27 . (NY): Academic Press; 2009. p. 1–23.
Diaz  RJ, Rosenberg  R. Spreading dead zones and consequences for marine ecosystems. Science. 2008:321 (5891 ):926–929. 10.1126/science.1156401.18703733
Dotti do Prado  F, Fernandez-Cebrian  R, Hashimoto  DT, Senhorini  JA, Foresti  F, Martinez  P, Porto-Foresti  F. Hybridization and genetic introgression patterns between two South American catfish along their sympatric distribution range. Hydrobiologia. 2017:788 (1 ):319–343. 10.1007/s10750-016-3010-5.
Duan  C . Hypoxia-inducible factor 3 biology: complexities and emerging themes. Am J Physiol Cell Physiol. 2016:310 (4 ):C260–C269. 10.1152/ajpcell.00315.2015.26561641
Eichstaedt  CA, Pagani  L, Antao  T, Inchley  CE, Cardona  A, Morseburg  A, Clemente  FJ, Sluckin  TJ, Metspalu  E, Mitt  M, et al  Evidence of early-stage selection on EPAS1 and GPR126 genes in Andean high altitude populations. Sci Rep. 2017:7 (1 ):13042. 10.1038/s41598-017-13382-4.29026132
Elbassiouny  AA, Buck  LT, Abatti  LE, Mitchell  JA, Crampton  WGR, Lovejoy  NR, Chang  BSW. Evolution of a novel regulatory mechanism of hypoxia inducible factor in hypoxia-tolerant electric fishes. J Biol Chem. 2024:300 (3 ):105727. 10.1016/j.jbc.2024.105727.38325739
Epstein  AC, Gleade  JM, McNeill  LA, Hewitson  KS, O’Rourke  J, Mole  DR, Mukherji  M, Metzen  E, Wilson  MI, Dhanda  A, et al  C. elegans EGL-9 and mammalian homologs define a family of dioxygenases that regulate HIF by prolyl hydroxylation. Cell. 2001:107 (1 ):43–54. 10.1016/S0092-8674(01)00507-4.11595184
Farhat  E, Talarico  GGM, Gregoire  M, Weber  J-M, Mennigen  JA. Epigenetic and post-transcriptional repression support metabolic suppression in chronically hypoxic goldfish. Sci Rep. 2022:12 (1 ):5576. 10.1038/s41598-022-09374-8.35368037
Farrell  AP, Richards  JG. Defining hypoxia: an integrative synthesis of the responses of fish to hypoxia. In: Richards  JG, Farrell  AP, Brauner  CJ, editors. Fish physiology. Vol. 27 . (NY): Academic Press; 2009. p. 487–503.
Fong  G-H, Takeda  K. Role and regulation of prolyl hydroxylase domain proteins. Cell Death Differ. 2008:15 (4 ):635–641. 10.1038/cdd.2008.10.18259202
Gasanov  EV, Jedrychowska  J, Kuznicki  J, Korzh  V. Evolutionary context can clarify gene names: teleosts as a case study. BioEssays. 2021:43 (6 ):2000258. 10.1002/bies.202000258.
Graham  AM, Presnell  JS. Hypoxia inducible factor (HIF) transcription factor family expansion, diversification, divergence and selection in eukaryotes. PLoS One. 2017:12 (6 ):e0179545. 10.1371/journal.pone.0179545.28614393
Guan  L, Chi  W, Xiao  W, Chen  L, He  S. Analysis of hypoxia-inducible factor alpha polyploidization reveals adaptation to Tibetan plateau in the evolution of schizothoracine fish. BMC Evol Biol. 2014:14 (1 ):192. 10.1186/s12862-014-0192-1.25205386
Healy  TM, Brennan  RS, Whitehead  A, Schulte  PM. Tolerance traits related to climate change resilience are independent and polygenic. Glob Change Biol. 2018:24 (11 ):5348–5360. 10.1111/gcb.14386.
Henikoff  S, Henikoff  JG. Amino acid substitution matrices from protein blocks. Proc Natl Acad Sci U S A. 1992:89 (22 ):10915–10919. 10.1073/pnas.89.22.10915.1438297
Hoegg  S, Brinkmann  H, Taylor  JS, Meyer  A. Phylogenetic timing of the fish-specific genome duplication correlates with the diversification of teleost fish. J Mol Evol. 2004:59 (2 ):190–203. 10.1007/s00239-004-2613-z.15486693
Babin  CH, Leiva  FP, Verberk  WCEP, Rees  BB. Paper data and code of the manuscript: evolution of key oxygen-sensing genes is associated with hypoxia tolerance in fishes. Zenodo; 2023. 10.5281/zenodo.10026473.
Verberk  WCEP, Sandkler  JF, van de Pol  ILE, Urbina  MA, Wilson  RW, McKenzie  DJ, Leiva  FP. Data and code for: body mass and cell size shape the tolerance of fishes to low oxygen in a temperature-dependent manner. Zenodo; 2022b. 10.5281/zenodo.6123770.
Hughes  LC, Orti  G, Huang  Y, Sun  Y, Baldwin  CC, Thompson  AW, Arcila  D, Betancur-R  R, Li  C, Becker  L, et al  Comprehensive phylogeny of ray-finned fishes (Actinopterygii) based on transcriptomic and genomic data. Proc Natl Acad Sci U S A. 2018:115 (24 ):6249–6254. 10.1073/pnas.1719358115.29760103
Ivan  M, Kaelin  WG  Jr. The EGLN-HIF O2-sensing system: multiple inputs and feedbacks. Mol Cell. 2017:66 (6 ):772–779. 10.1016/j.molcel.2017.06.002.28622522
Ivan  M, Kondo  K, Yang  H, Kim  W, Valiando  J, Ohh  M, Salic  A, Asara  JM, Lane  WS, Kaelin  WG  Jr. HIFalpha targeted for VHL-mediated destruction by proline hydroxylation: implications for O2 sensing. Science. 2001:292 (5516 ):464–468. 10.1126/science.1059817.11292862
Jaakkola  P, Mole  DR, Tian  YM, Wilson  MI, Gielbert  J, Gaskell  SJ, von Kriegsheim  A, Hebestreit  HF, Mukherji  M, Schofield  CJ, et al  Targeting of HIF-alpha to the von Hippel-Lindau ubiquitylation complex by O2-regulated prolyl hydroxylation. Science. 2001:292 (5516 ):468–472. 10.1126/science.1059796.11292861
Jane  SF, Hansen  GJA, Kraemer  BM, Leavitt  PR, Mincer  JL, North  RL, Pilla  RM, Stetler  JT, Williamson  CE, Woolway  RI, et al  Widespread deoxygenation of temperate lakes. Nature. 2021:594 (7861 ):66–70. 10.1038/s41586-021-03550-y.34079137
Jombart  T . Adegenet: a R package for the multivariate analysis of genetic markers. Bioinform. 2008:24 (11 ):1403–1405. 10.1093/bioinformatics/btn129.
Jombart  T, Ahmed  I. Adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinform. 2011:27 (21 ):3070–3071. 10.1093/bioinformatics/btr521.
Jombart  T, Devillard  S, Balloux  F. Disciminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010:11 (1 ):94. 10.1186/1471-2156-11-94.20950446
Kaelin  WG  Jr, Ratcliffe  PJ. Oxygen sensing by metazoans: the central role of the HIF hydroxylase pathway. Mol Cell. 2008:30 (4 ):393–402. 10.1016/j.molcel.2008.04.009.18498744
Kajungiro  RA, Palaiokostas  C, Lopes Pinto  FA, Mmochi  AJ, Mtolera  M, Houston  RD, de Koning  DJ. Population structure and genetic diversity of Nile tilapia (Oreochromis niloticus) strains cultured in Tanzania. Front Genet. 2019:10 :1269. 10.3389/fgene.2019.01269.31921307
Keith  B, Johnson  R, Simon  M. HIF1α and HIF2α: sibling rivalry in hypoxic tumour growth and progression. Nat Rev Cancer. 2012:12 (1 ):9–22. 10.1038/nrc3183.
Kosakovsky Pond  SL, Frost  SDW. Datamonkey: rapid detection of selective pressure on individual sites of codon alignments. Bioinform. 2005a:21 (10 ):2531–2533. 10.1093/bioinformatics/bti320.
Kosakovsky Pond  SL, Frost  SDW. Not so different after all: a comparison of methods for detecting amino acid sites under selection. Mol Biol Evol. 2005b:22 (5 ):1208–1222. 10.1093/molbev/msi105.15703242
Kosakovsky Pond  SL, Frost  SDW, Muse  SV. Hyphy: hypothesis testing using phylogenies. Bioinform. 2005:21 (5 ):676–679. 10.1093/bioinformatics/bti079.
Kosakovsky Pond  SL, Poon  AFY, Velazquez  R, Weaver  S, Hepler  NL, Murrell  B, Shank  SD, Magalis  BR, Bouvier  D, Nekrutenko  A, et al  Hyphy 2.5—a customizable platform for evolutionary hypothesis testing using phylogenies. Mol Biol Evol. 2020:37 (1 ):295–299. 10.1093/molbev/msz197.31504749
Kuang  Y-Y, Zheng  X-H, Li  C-Y, Li  X-M, Cao  D-C, Tong  G-X, Lv  W-H, Xu  W, Zhou  Y, Zhang  X-F, et al  The genetic map of goldfish (Carassius auratus) provided insights to the divergent genome evolutions in the Cyprinidae family. Sci Rep. 2016:6 (1 ):34849. 10.1038/srep34849.27708388
Li  Y, Wu  D-D, Boyko  AR, Wang  G-D, Wu  S-F, Irwin  DM, Zhang  Y-P. Population variation revealed high-altitude adaptation of Tibetan Mastiffs. Mol Biol Evol. 2014:31 (5 ):1200–1205. 10.1093/molbev/msu070.24520091
Li  X, Zhang  M, Ling  C, Sha  H, Zou  G, Liang  H. Molecular characterization and response of prolyl hydroxylase domain (PHD) genes to hypoxia stress in Hypophthalmichthys molitrix. Animals. 2022:12 (2 ):131. 10.3390/ani12020131.35049755
Lian  S, Zhou  Y, Liu  Z, Gong  A, Cheng  L. The differential expression patterns of paralogs in response to stresses indicate expression and sequence divergences. BMC Plant Biol. 2020:20 (1 ):277. 10.1186/s12870-020-02460-x.32546126
Liao  B-Y, Zhang  J. Low rates of expression profile divergence in highly expressed genes and tissue-specific genes during mammalian evolution. Mol Biol Evol. 2006:23 (6 ):1119–1128. 10.1093/molbev/msj119.16520335
Lynch  M, Conery  JS. The evolutionary fate and consequences of duplicate genes. Science. 2000:290 (5494 ):1151–1155. 10.1126/science.290.5494.1151.11073452
Macqueen  DJ, Johnston  IA. A well-constrained estimate for the timing of the salmonid whole genome duplication reveals major decoupling from species diversification. Proc Biol Sci.  2014:281 (1778 ):20132881. 10.1098/rspb.2013.2881.24452024
Mandic  M, Best  C, Perry  SF. Loss of hypoxia-inducible factor 1α affects hypoxia tolerance in larval and adult zebrafish (Danio rerio). Proc Biol Sci.  2020:287 (1927 ):20200798. 10.1098/rspb.2020.0798.32453991
Mandic  M, Joyce  W, Perry  SF. The evolutionary and physiological significance of the Hif pathway in teleost fishes. J Exp Biol. 2021:224 (18 ):jeb231936. 10.1242/jeb.231936.34533194
Mandic  M, Regan  MD. Can variation among hypoxic environments explain why different fish species use different hypoxic survival strategies?  J Exp Biol. 2018:221 (21 ):jeb161349. 10.1242/jeb.161349.30381477
Maxwell  PH, Wiesener  MS, Chang  GW, Clifford  SC, Vaux  EC, Cockman  ME, Wykoff  CC, Pugh  CW, Maher  ER, Ratcliffe  PJ. The tumour suppressor protein VHL targets hypoxia-inducible factors for oxygen-dependent proteolysis. Nature. 1999:399 (6733 ):271–275. 10.1038/20459.10353251
Murrell  B, Weaver  S, Smith  MD, Wertheim  JO, Murrell  S, Aylward  A, Eren  T, Pollner  T, Martin  DP, Smith  DM, et al  Gene-wide identification of episodic selection. Mol Biol Evol. 2015:32 (5 ):1365–1371. 10.1093/molbev/msv035.25701167
Murrell  B, Wertheim  JO, Moola  S, Weighill  T, Scheffler  K, Kosakovsky Pond  SL. Detecting individual sites subject to episodic diversifying selection. PLoS Genet. 2012:8 (7 ):e1002764. 10.1371/journal.pgen.1002764.22807683
Nelson  JS, Grande  TC, Wilson  MVH. Fishes of the world. Hoboken (NJ), USA: John Wiley & Sons, Inc.; 2016.
Nikinmaa  M, Rees  BB. Oxygen-dependent gene expression in fishes. Am J Physiol Regul Integr Comp Physiol. 2005:288 (5 ):R1079–R1090. 10.1152/ajpregu.00626.2004.15821280
Ohno  S . Evolution by gene duplication. Heidelberg, Germany: Springer Berlin; 1970.
Orme  D, Freckleton  R, Thomas  G, Petzoldt  T, Fritz  S, Isaac  N, Pearse  W. The Caper Package: Comparative Analyses of Phylogenetics and Evolution in R. Version 1.0.1.  https://cran.r-project.org/package=caper.
Pamenter  ME, Hall  JE, Tanabe  Y, Simonson  TS. Cross-species insights into genomic adaptations to hypoxia. Front Genet. 2020:22 (11 ) :743. 10.3389/fgene.2020.00743.
Pan  W, Godoy  RS, Cook  DP, Scott  AL, Nurse  CA, Jonz  MG. Single-cell transcriptomic analysis of neuroepithelial cells and other cell types of the gills of zebrafish (Danio rerio) exposed to hypoxia. Sci Rep. 2022:12 (1 ):10144. 10.1038/s41598-022-13693-1.35710785
Paradis  E, Schliep  K. Ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinform. 2019:35 (3 ):526–528. 10.1093/bioinformatics/bty633.
Patel  SA, Simon  MC. Biology of hypoxia-inducible factor-2alpha in development and disease. Cell Death Differ. 2008:15 (4 ):628–634. 10.1038/cdd.2008.17.18259197
Pennell  M, Eastman  J, Slater  G, Brown  J, Uyeda  J, Fitzjohn  R, Alfaro  M, Harmon  L. Geiger v2.0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees. Bioinform. 2014:30 (15 ):2216–2218. 10.1093/bioinformatics/btu181.
Postlethwait  JH, Woods  IG, Ngo-Hazelett  P, Yan  YL, Kelly  PD, Chu  F, Huang  H, Hill-Force  A, Talbot  WS. Zebrafish comparative genomics and the origins of vertebrate chromosomes. Genome Res. 2000:10 (12 ):1890–1902. 10.1101/gr.164800.11116085
Powell  WH, Hahn  ME. Identification and functional characterization of hypoxia-inducible factor 2α from the estuarine teleost, Fundulus heteroclitus: interaction of HIF-2α with two ARNT2 splice variants. J Exp Zool. 2002:294 (1 ):17–29. 10.1002/jez.10074.11932946
Ranwez  V, Douzery  EJP, Cambon  C, Chantret  N, Delsuc  F. MACSE v2: toolkit for the alignment of coding sequences accounting for frameshifts and stop codons. Mol Biol Evol. 2018:35 (10 ):2582–2584. 10.1093/molbev/msy159.30165589
Ranwez  V, Harispe  S, Delsuc  F, Douzery  EJP. MACSE: multiple alignment of coding sequences accounting for frameshifts and stop codons. PLoS One. 2011:6 (9 ):e22594. 10.1371/journal.pone.0022594.21949676
R Development Core Team . R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2021.
Reemeyer  JE, Rees  BB. Standardizing the determination and interpretation of Pcrit in fishes. J Exp Biol. 2019:222 :jeb210633. 10.1242/jeb.210633.31511343
Rees  BB, Matute  LA. Repeatable interindividual variation in hypoxia tolerance in the gulf killifish, Fundulus grandis. Physiol Biochem Zool. 2018:91 (5 ):1046–1056. 10.1086/699596.30125141
Regan  MD, Mandic  M, Dhillon  RS, Lau  GY, Farrell  AP, Schulte  PM, Seibel  BA, Speers-Roesch  B, Ultsch  GR, Richards  JG. Don’t throw the fish out with the respirometry water. J Exp Biol. 2019:222 (6 ):jeb200253. 10.1242/jeb.200253.30923073
Revell  L . Phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol Evol. 2012:3 (2 ):217–223. 10.1111/j.2041-210X.2011.00169.x.
Rogers  NJ, Urbina  MA, Reardon  EE, McKenzie  DJ, Wilson  RW. A new analysis of hypoxia tolerance in fishes using a database of critical oxygen level (Pcrit). Conserv Physiol. 2016:4 (1 ):cow012. 10.1093/conphys/cow012.27293760
Rytkonen  KT, Akbarzadeh  A, Miandare  HK, Kamei  H, Duan  C, Leder  EH, Williams  TA, Nikinmaa  M. Subfunctionalization of cyprinid hypoxia-inducible factors for roles in development and oxygen sensing. Evolution. 2013:67 (3 ):873–882. 10.1111/j.1558-5646.2012.01820.x.23461336
Rytkonen  KT, Vuori  KAM, Primmer  CR, Nikinmaa  M. Comparison of hypoxia-inducible factor-1 alpha in hypoxia-sensitive and hypoxia-tolerant fish species. Comp Biochem Physiol Part D Genomics Proteomics. 2007:2 :177–186. 10.1016/j.cbd.2007.03.001.20483291
Rytkonen  KT, Williams  TA, Renshaw  GM, Primmer  CR, Nikinmaa  M. Molecular evolution of the metazoan PHD-HIF oxygen-sensing system. Mol Biol Evol. 2011:28 (6 ):1913–1926. 10.1093/molbev/msr012.21228399
Sacerdot  C, Louis  A, Bon  C, Berthelot  C, Roest Crollius  H. Chromosome evolution at the origin of the ancestral vertebrate genome. Genome Biol. 2018:19 (1 ):166. 10.1186/s13059-018-1559-1.30333059
Sampaio  E, Santos  C, Rosa  IC, Ferreira  V, Portner  H-O, Duarte  CM, Levin  LA, Rosa  R. Impacts of hypoxic events surpass those of future ocean warming and acidification. Nat Ecol Evol. 2021:5 (3 ):311–321. 10.1038/s41559-020-01370-3.33432134
Sandberg  M, Eriksson  L, Jonsson  J, Sjostrom  M, Wold  S. New chemical descriptors relevant for the design of biologically active peptides: a multivariate characterization of 87 amino acids. J Med Chem. 1998:41 (14 ):2481–2491. 10.1021/jm9700575.9651153
Schweizer  RM, Velotta  JP, Ivy  CM, Jones  MR, Muir  SM, Bradburd  GS, Storz  JF, Scott  GR, Cheviron  ZA. Physiological and genomic evidence that selection on the transcription factor Epas1 has altered cardiovascular function in high-altitude deer mice. PLoS Genet. 2019:15 (11 ):e1008420. 10.1371/journal.pgen.1008420.31697676
Seibel  BA, Andres  A, Birk  MA, Burns  AL, Shaw  CT, Timpe  AW, Welsh  CJ. Oxygen supply capacity breathes new life into critical oxygen partial pressure (Pcrit). J Exp Biol. 2021:224 (8 ):jeb242210. 10.1242/jeb.242210.33692079
Semenza  GL . Regulation of oxygen homeostasis by hypoxia-inducible factor 1. Physiology. 2009:24 (2 ):97–106. 10.1152/physiol.00045.2008.19364912
Semenza  GL . Oxygen sensing, homeostastis, and disease. N Engl J Med. 2011:365 (6 ):537–547. 10.1056/NEJMra1011165.21830968
Semenza  GL . Hypoxia-inducible factors in physiology and medicine. Cell. 2012:148 (3 ):399–408. 10.1016/j.cell.2012.01.021.22304911
Shedko  SV, Miroshnichenko  IL, Nemkova  GA. Phylogeny of salmonids (salmoniformes: Salmonidae) and its molecular dating: analysis of mtDNA data. Russ J Genetics. 2013:49 (6 ):623–637. 10.1134/S1022795413060112.
Sollid  J, De Angelis  P, Gundersen  K, Nilsson  GE. Hypoxia induces adaptive and reversible gross morphological changes in crucian carp gills. J Exp Biol. 2003:206 (20 ):3667–3673. 10.1242/jeb.00594.12966058
Sollid  J, Weber  RE, Nilsson  GE. Temperature alters the respiratory surface area of crucian carp Carassius carassius and goldfish Carassius auratus. J Exp Biol. 2005:208 (6 ):1109–1116. 10.1242/jeb.01505.15767311
Stamatakis  A . RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinform. 2014:30 (9 ):1312–1313. 10.1093/bioinformatics/btu033.
Sunde  J, Yildirim  Y, Tibblin  P, Forsman  A. Comparing the performance of microsatellites and RADseq in populations genetic studies: analysis of data for pike (Esox lucius) and a synthesis of previous studies. Front Genet. 2020:11 :218. 10.3389/fgene.2020.00218.32231687
Townley  IK, Babin  CH, Murphy  TE, Summa  CM, Rees  BB. Genomic analysis of hypoxia inducible factor alpha in ray-finned fishes reveals missing Ohnologs and evidence of widespread positive selection. Sci Rep. 2022:12 (1 ):22312. 10.1038/s41598-022-26876-7.36566251
Ultsch  GR, Jackson  DC, Moalli  R. Metabolic oxygen conformity among lower vertebrates: the toadfish revisited. J Comp Physiol. 1981:142 (4 ):439–443. 10.1007/BF00688973.
Verberk  WCEP, Bilton  DT, Calosi  P, Spicer  JI. Oxygen supply in aquatic ectotherms: partial pressure and solubility together explain biodiversity and size patterns. Ecology. 2011:92 (8 ):1565–1572. 10.1890/10-2369.1.21905423
Verberk  WCEP, Sandker  JF, van de Pol  ILE, Urbina  MA, Wilson  RW, McKenzie  DJ, Leiva  FP. Body mass and cell size shape the tolerance of fishes to low oxygen in a temperature-dependent manner. Glob Change Biol. 2022a:28 (19 ):5695–5707. 10.1111/gcb.16319.
Volff  J-N . Genome evolution and biodiversity in teleost fish. Heredity. 2005:94 (3 ):280–294. 10.1038/sj.hdy.6800635.15674378
Wang  Y, Li  X-Y, Xu  W-J, Wang  K, Wu  B, Xu  M, Chen  Y, Miao  L-J, Wang  Z-W, Li  Z, et al  Comparative genome anatomy reveals evolutionary insights into a unique amphitriploid fish. Nat Ecol Evol. 2022:6 (9 ):1354–1366. 10.1038/s41559-022-01813-z.35817827
Wang  Y, Yang  L, Wu  B, Song  Z, He  S. Transcriptome analysis of the plateau fish (Triplophysa dalaica): implications for adaptation to hypoxia in fishes. Gene. 2015a:565 (2 ):211–220. 10.1016/j.gene.2015.04.023.25869933
Wang  Y, Yang  L, Zhou  K, Zhang  Y, Song  Z, He  S. Evidence for adaptation to the Tibetan Plateau inferred from Tibetan loach transcriptomes. Genome Biol Evol. 2015b:7 (11 ):2970–2982. 10.1093/gbe/evv192.26454018
Weaver  S, Shank  SD, Spielman  SJ, Li  M, Muse  SV, Kosakovsky Pond  SL. Datamonkey 2.0: a modern web application for characterizing selective and other evolutionary processes. Mol Biol Evol. 2018:35 (3 ):773–777. 10.1093/molbev/msx335.29301006
Wenger  RH, Stiehl  DP, Camenisch  G. Integration of oxygen signaling at the consensus HRE. Sci STKE. 2005:2005 (306 ):re12. 10.1126/stke.3062005re12.16234508
Wickham  H, Averick  M, Bryan  J, Chang  W, McGowan  LD, Francois  R, Grolemund  G, Hayes  A, Henry  L, Hester  J, et al  Welcom to the tidyverse. J Open Source Softw. 2019:4 (43 ):1686. 10.21105/joss.01686.
Wong  KM, Suchard  MA, Huelsenbeck  JP. Alignment uncertainty and genomic analysis. Science. 2008:319 (5862 ):473–476. 10.1126/science.1151532.18218900
Wood  CM . The fallacy of the Pcrit—are there more useful alternatives?  J Exp Biol. 2018:221 (22 ):jeb163717. 10.1242/jeb.163717.30420494
Woods  HA, Moran  AL, Atkinson  D, Audxijonyte  A, Berenbrink  M, Borges  FO, Burnett  KG, Burnett  LE, Coates  CJ, Collin  R, et al  Integrative approaches to understanding organismal responses to aquatic deoxygenation. Biol Bull.  2022:243 (2 ):85–103. 10.1086/722899.36548975
Xenopoulos  MA, Barnes  RT, Boodoo  KS, Butman  D, Catalan  N, D’Amario  SC, Fasching  C, Kothawala  DN, Pisani  O, Solomon  CT, et al  How humans alter dissolved organic matter composition in freshwater: relevance for the Earth's biogeochemistry. Biogeochemistry. 2021:154 (2 ):323–348. 10.1007/s10533-021-00753-3.
Yu  Y, He  J, Liu  W, Li  Z, Weng  S, He  J, Guo  C. Molecular characterization and functional analysis of hypoxia-responsive factor prolyl hydroxylase domain 2 in mandarin fish (Siniperca chuatsi). Animals. 2023:13 (9 ):1556. 10.3390/ani13091556.37174593
Zhang  J, Dong  C, Feng  J, Li  J, Li  S, Feng  J, Duan  X, Sun  G, Xu  P, Li  X. Effects of dietary supplementation of three strains of Lactococcus lactis on HIFs genes family expression of the common carp following Aeromonas hydrophila infection. Fish Shellfish Immunol. 2019:92 :590–599. 10.1016/j.fsi.2019.06.040.31252044
Zhao  S-S, Su  X-L, Yang  H-Q, Zheng  G-D, Zou  S-M. Functional exploration of SNP mutations in HIF2αb gene correlated with hypoxia tolerance in blunt snout bream (Megalobrama amblycephala). Fish Physiol Biochem. 2023:49 (2 ):239–251. 10.1007/s10695-023-01173-w.36859574
