For aerobic organisms, both the Hypoxia-Inducible Factor (HIF) pathway and the mitochondrial genomes are key players in regulating oxygen homeostasis. However, recent work has suggested that these mechanisms are not as highly conserved as previously thought, prompting more thorough surveys across animal higher taxonomic levels, which would in turn permit testing of hypotheses about the ecological conditions that may have facilitated evolutionary loss of such genes. The phylum Cnidaria is known to harbor wide variation in mitochondrial genome morphology, from typical single circular chromosomes to fragmented linear chromosomes. More recently, members of the cnidarian clade Myxozoa, comprising obligate endoparasites, were shown to have lost their mitochondrial genome, suggesting that variation in environmental oxygen availability may be a key determinant in the evolution of metabolic gene networks. Here, we surveyed genomes and transcriptomes across 42 cnidarian species for the presence of HIF pathway members (HIFa, EGLN, VHL), as well as for an assortment of hypoxia, mitochondrial, and stress-response toolkit genes. We find that presence of the HIF pathway, as well as number of genes associated with mitochondria, hypoxia, and stress response, do not vary based on mitochondrial genome morphology. More interestingly, we uncover evidence that myxozoans have lost the canonical HIF pathway repression machinery, potentially altering HIF pathway functionality to work under the specific conditions of their parasitic lifestyles. In addition, relative to other cnidarians, myxozoans show loss of large proportions of genes associated with the mitochondrion (~39%), and involved in response to hypoxia (~27.5%) and general stress (~32%). Our results provide additional evidence that the HIF regulatory machinery is evolutionarily labile and that variations in the canonical system have evolved in many animal groups.

Of the 35 genome assemblies publicly available for Cnidaria, 17 had annotated protein sequence files that were used in this study. These included classes Anthozoa (n = 9), Cubozoa (1), Hydrozoa (2), Myxozoa (3), and Scyphozoa (2) and were downloaded directly from their respective repositories (Supplemental Table 1; available online). To increase taxon sampling, transcriptomic data (Illumina paired-end mRNA reads) from an additional 25 species were downloaded from the NCBI Sequence Read Archive (SRA) and assembled de novo. Among these were species in the classes Cubozoa (4), Hydrozoa (9), Myxozoa (5), and Scyphozoa (4) and Staurozoa (3) (Supplemental Table 2; available online). The singular species of the group Polylpodium was also included.

For de novo transcriptome assemblies, downloaded reads were trimmed for base quality (>Q30) and adapters using BBDuk (Bushnell et al., 2017), and transcripts assembled de novo with Trinity (Haas et al., 2013) using default parameters. The assembled transcripts were then translated to protein sequences using Transdecoder, which included steps for predicting all open reading frames (ORFs), BLASTing of ORFs against the Uniprot/Swissprot database (release date March 2021), and retaining the best protein models for each transcript. For transcriptome assemblies from the endoparasitic myxozoan species, contamination from host tissue (fish or annelid) is likely to be present in the pool of reads, as observed for mRNA-based assemblies from parasitic copepods (Graham and Barreto, 2020). We hence filtered the myxozoan assemblies by identifying transcripts assembled from potential host RNA. We prepared a custom BLAST database containing protein sequences from three fish and one annelid obtained from Ensembl. These were Oncorhynchus mykiss (assembly Omyk_1.0, accession GCA_002163495.1), Cyprinus carpio (common_carp_genome, GCA_000951615.2), Sparus aurata (SpaAur1.1, GCA_900880675.1), and Capitella teleta (Capitella telata v1.0, GCA_000328365.1). We then performed BLASTP searches of our assembled myxozoan sequences against this database. Hits that had a match within these potential host sequences, with e-value ≤ 1e-75 and percent identity ≥ 85, were considered contaminant and removed from the assemblies, as used by Yahalomi et al. (2020).

Includes all transcriptomes used in this study - both de novo generated and from other sources.