We examined transcriptional gene regulation of blood by performing statistical overrepresentation analyses of Gene Ontology (GO) terms using PANTHERdb (Mi et al., 2019; Thomas et al., 2022). Accordingly, we scored annotated genes demonstrating significant two-fold (or higher) differences in gene expression in each cave population relative to surface fish. This resulted in six analyses, i.e., three pair-wise comparisons performed for both over- and under-expression. GO terms enrichments (FDR adjusted p-value <0.05) from each analysis were compared to identify convergent/divergent expression patterns between cave populations. A prior study using the same dataset provided expression validation through analysis of five genes subjected to quantitative real-time PCR (qPCR) and calculated delta Cq using ssr3 as the reference gene [see Boggs and Gross (2025)]. This study revealed an average correlation coefficient of 0.89, indicating a strong relationship (Cohen et al., 2009) and validation of our RNA-seq dataset.

Interestingly, we discovered substantial overlap between Tinaja and Pachón, and Pachón and Chica (Figure 1). These overlapping sets included every identified GO term for Pachón, however no overlap was observed for Chica and Tinaja cavefish. Many GO terms shared between Tinaja and Pachón cavefish were significant due to hemoglobin genes [see Boggs and Gross (2025)] including: oxygen transport, heme binding, oxygen binding, and hemoglobin complex (Figure 1). Overlapping terms between Pachón and Chica mostly reflected terms associated with immune system function, including: defense response to bacterium, innate immune response, phagocytosis recognition, and antigen binding.

These results were not entirely surprising given that prior GO enrichment studies in Astyanax identified convergent mechanisms of cave adaptation, including broad development processes (Riddle et al., 2020), metabolism (Krishnan et al., 2020), and immunity (Krishnan et al., 2022). Additionally, the number and diversity of genes within a test list affects the outcome of overrepresentation studies (Wijesooriya et al., 2022). Given that whole blood is a highly complex tissue (capable of predicting an estimated 60% of gene expression for dozens of tissues) (Basu et al., 2021), this likely impacted the statistical outcomes of the analysis.

We further aimed to investigate genes of potential interest that may not have been detected in these GO analyses. Accordingly, we created four lists representing genes that are biologically-relevant to hypoxia including: genes expressed higher or lower in all examined cave populations compared to surface fish (Table 1) and genes expressed higher and lower in Tinaja and Pachón compared to surface fish (Table 2), while excluding Chica, given their similarity in hemoglobin expression. A literature search for each of these genes was conducted to provide any potential relevance to adaptation to hypoxic caves.

Consistent with prior findings (Boggs and Gross, 2025), many hemoglobin genes were expressed higher in cavefish compared to surface fish (Tables 1 and 2) with the vast majority expressed higher only in Tinaja and Pachón (relative to Chica and Surface). We also identified numerous genes associated with canonical hypoxia response pathways. Interestingly, many genes normally activating these pathways were expressed lower in cavefish compared to surface fish, and genes typically suppressing these pathways were expressed higher in cavefish (Tables 1, 2). Notably, two hypoxia inducible factor (hif) genes (hif1al2 and hif1ab) and multiple genes contributing to HIF signaling, including cathepsin Ba (ctsba), lysosomal associated membrane protein 2 (lamp2), splicing factor 3b subunit 1 (sf3b1), eukaryotic translation initiation factor 5 (eif5), and p450 (cytochrome) oxidoreductase a (pora) were expressed lower in cavefish compared to surface fish. Additionally, two genes known to suppress HIF signaling (tp53inp1 and tcf20) were expressed higher in cavefish compared to surface fish.

In normoxic conditions, hif is continuously transcribed, but is controlled post-translationally by prolyl hydroxylase (PHD) and von Hippel-Lindau (VHL) proteins. During hypoxia, PHD activity is inhibited and Hif is not degraded. Thus, hif transcript abundance is not necessarily representative of Hif activity in mammals (Semenza and Wang, 1992; Maxwell et al., 1999; Ivan et al., 2001; Jaakkola et al., 2001; Bruick and McKnight, 2001; Epstein et al., 2001). In the Chinese sucker (Myxocyprinus asiaticus), a study revealed increased hif transcription is required to prevent degradation of Hif during hypoxia (Chen et al., 2012). Having said this, a recent study uncovered diverse reports of hifa mRNA abundance in fish exposed to hypoxia, as a likely function of varying methodologies (Murphy and Rees, 2024). Nevertheless, elasmobranch fish conditioned to hypoxia express hif higher than individuals that have not experienced hypoxia (Rytkönen et al., 2012). Additionally, certain hif family members are expressed higher in Pachón cavefish embryos (after normoxic rearing or exposure to hypoxia) than in surface fish (van der Weele and Jeffery, 2022). Thus, we were initially surprised to find that adult cavefish express two hif family members much lower than surface fish and express other known hypoxia response genes in similar, counterintuitive, patterns. In light of varying reports of hifa transcription in fish (Murphy and Reese, 2024), it will be essential to better characterize protein levels of hif1a in forthcoming studies through the use of Western blot analyses.

One explanation for these observed patterns may be the negative consequences associated with prolonged expression of hif. Hif is linked to many human pathologies including tumorigenesis, cardiovascular, metabolic, and reproductive diseases [reviewed in Chen et al. (2020)]. In mice, pharmacological knock-down of Hif protein relieved symptoms of rheumatoid arthritis (Hu et al., 2020). Hif pathways can also impair major histocompatibility complex function in culture, leading to an inability to recognize and eliminate cancerous and other harmful cells (Sethumadhavan et al., 2017). Additionally, Hif proteins influence ion fluctuations and homeostasis in fish, a well-characterized mechanism to conserve energy during hypoxia [reviewed in Pelster and Egg (2018)]. Thus, future work in Astyanax may determine if downregulation of hif and other known hypoxia response pathways are advantageous in cavefish to save energy, maintain proper immune function, and prevent disease and inflammation.

Prior work in Astyanax revealed that oxygen levels are a good deal lower in cave waters compared to surface waters (Boggs and Gross, 2025). An important consideration for this study is the fact that all experimental animals were reared in normoxic conditions. Indeed, our putative Chica cavefish were acquired from a commercial vendor, and therefore it is not possible for us to determine the extent to which transcriptomic changes are a function of assimilation to captivity. Interestingly, a number of cave populations maintain significantly elevated levels of hemoglobin despite rearing in normoxia for generations (Boggs and Gross, 2025). Given that the transcriptome can change markedly when comparing captive-bred versus wild-caught individuals (Krishnan et al., 2020), an essential future direction for this research includes examination of the blood transcriptome from individuals drawn from the natural population.

Astyanax cave and surface fish were reared in a satellite aquatic facility at the University of Cincinnati within a custom-designed reverse osmosis husbandry unit comprised of 5- and 10- gallon continuous flow tanks (Aquaneering, San Diego, CA). Animals were exposed to a 12:12 h light: dark cycle and fed a slurry of dry flake food (TetraMin Pro) and system water daily. Water in this system is processed through a series of filters including UV, 25-micron polypropylene felt, activated carbon, and dense particulate. Additionally, water conditions were adjusted using real-time dose monitoring of sodium bicarbonate and Instant Ocean sea salt to conductivity of 750 μS/cm (±50 μS/cm) and pH of 7.4 (±0.2). Water temperature was kept at 24°C (±2°C). Importantly, dissolved oxygen was not manipulated for this study meaning all fish were exposed to ample oxygen.

The surface fish, Pachón cavefish and Tinaja cavefish used in this study were derived from breeding adults originally provided to our lab by Dr. Richard Borowsky (New York University). Specifically, the pedigrees used included Asty-152 and Asty-155 (surface fish), Asty-163 and Asty-138 (Pachón cavefish), and Asty-19 (Tinaja cavefish). Surface fish are descended from wild-caught individuals from the Río Sabinas and Río Valles drainages near Ciudad Valles in San Luis Potosí, Mexico. All Chica cavefish were acquired from the commercial pet trade. We extracted whole blood from (n = 4) surface, Pachón, Tinaja, and Chica populations (total n = 16) via the caudal vein using 31G syringes (BD Ultra-Fine™, BD Biosciences, San Jose, CA). In order to limit any potential effects outside the scope of this study, two male and two female fish were used from each population, fish were post-breeding age, and whole blood extractions were completed between 12:00 p.m - 1:00 p.m. All procedures were conducted in accordance with University of Cincinnati IACUC (Protocol# 22-01-06-01).

Immediately following whole blood extraction, whole RNA was isolated using an RNeasy Universal Mini Kit (Qiagen, Germantown, MD) according to the manufacturer’s directions. All RNA samples were subjected to quantification using a Nanodrop Lite spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA). The purity of samples was estimated based on the A260/A280 ratio, and only RNA samples of ∼2.0 were submitted for sequencing. Owing to the technical requirements of RNA-sequencing, samples had to be pooled by population. To mitigate potential effects of sequencing error, each pool was sub-aliquoted into three technical replicates with each replicate (n = 12) containing the same volume of RNA. Pools were submitted to the DNA Core at Cincinnati Childrens’ Hospital and Medical Center. There, additional RNA QC was conducted, polyA stranded libraries were generated, QC was conducted on the libraries, and they were subject to sequencing using an Illumina HiSeq 2,500 sequencer. This resulted in twenty-million 125bp paired end reads per sample. Raw reads were assessed for quality and length using FastQC (Wingett and Andrews, 2018) (version 0.11.8) and adapters were trimmed using Trimmomatic (Bolger et al., 2014) (version 0.39).

Analysis of gene expression was conducted by running a reference based analysis in CLC Genomics (Qiagen, Germantown, MD, version 12.0.1) using manufacturer recommended parameters. The latest Astyanax genome (AstMex3_surface, GCA_023375975.1) was used as the reference sequence. We used the latest annotation file for this reference from NCBI RefSeq (GCF_023375975.1, NCBI annotation release 103).

In order to increase efficiency of downstream transcriptome-wide analysis, we conducted a second RNA sequencing experiment in CLC Genomics using the “Astyanax-mexicanus-2.0” genome retrieved from Ensembl [GCA_000372685.2 (Warren et al., 2021)] as the reference and annotations from Ensembl release 106 were used to identify genes and determine expression. RNA-sequencing was validated using qPCR for five genes [see Boggs and Gross (2025)].

To investigate transcriptome-wide patterns of gene expression, we conducted a Gene Ontology (GO) Statistical Overrepresentation Test. Each cave population was assessed independently against surface fish. Thus, we created seven gene lists, one list representing genes expressed higher in a cave population versus surface fish, one list representing genes expressed lower in a cave population versus surface, and a list containing all genes detectable in this assay [noise threshold surpassed with TPM value of at least 2 (Wagner et al., 2013)]. Each list representing a comparison between a cave and surface population contained genes detectable for at least one of the two populations and with a fold change of at least 2x (any gene with a TPM value of 0 was substituted with the lowest TPM value in the entire dataset - 0.00192,433 in Tinaja fat1a–so that a fold change value could be calculated). We used PANTHERdb (Mi et al., 2019; Thomas et al., 2022) (version 17.0) to conduct a statistical overrepresentation test. Because Astyanax GO terms are not available in PANTHER, IDs in each list were converted to orthologous Danio rerio IDs by using BioMart (Smedley et al., 2009). We successfully converted 6,882 of 8,550 (∼80%) IDs from our and used these as our reference (Aleksander et al., 2023) for the statistical overrepresentation test. We used a Fisher’s Exact text to calculate p-values which were corrected using false discovery rate to determine statistical significance. Each of three categories of GO terms were assessed: biological process, molecular function, and cellular component. Results from each cave-to-surface analysis were then compared to determine convergence/divergence between cave populations.

In addition, we investigated genes of potential interest that may have been missed in the GO analysis. Thus, we compiled four additional lists of genes: two lists representing genes of putative biological relevance that are either expressed higher or lower in Chica, Tinaja, and Pachón cavefish compared to surface fish as well as two lists expressed higher or lower in Tinaja and Pachón compared to surface. Expression data derived from Chica cavefish was omitted from these lists owing to the difference in expression of hemoglobin compared to Tinaja and Pachón. Genes were ranked according to putative biological relevance. Rank was determined by subtracting the fold change (cavefish expression value divided by surface fish expression value) of a gene from each cavefish expression value and summing the absolute values from each cave population. Genes that have not been characterized were removed and the remaining genes were filtered for relevance to hypoxia using literature searches.