A pangenome framework uncovers the role of deletions in repeated evolution of cave-derived traits

Overview of study design

Our sequencing and filtering strategy enabled us to examine SVs that likely arose after cave populations diverged from surface populations and compare SVs across independently derived cave lineages. Overall, using a total of 120 A. mexicanus genomes from four surface and seven cave populations, we identified structural insertions and deletions originating from two independent evolutionary events (Fig. 1A–D). These populations were selected to provide a robust sampling of cave and surface fish in each lineage (n = 56 for lineage 1, n = 64 for lineage 2), which allows for inter- and intralineage comparisons of structural variation.

Figure 1.

Mexican tetra population overview. (A) Representation of surface- and cave-morph Mexican tetra fish. (B) Phylogenetic relationships of populations used in the present study (subset from full phylogeny constructed by Moran et al. 2023). Populations belong to two distinct lineages, with cavefish evolving in each lineage. (C) Map of cave locations and surface population sampling points, San Luis Potosí and Tamaulipas, Mexico. (D) Sequencing data sets and sample sizes for each population.

Open in new tabDownload PowerPointLink to figure

To identify SVs, we first generated a pangenome graph from four reference genomes (Imarazene et al. 2021; Warren et al. 2024) and identified insertion and deletion variants within the graph (Fig. 2A). We refer to these variants as pangenome SVs and characterize them across the genome. To enable further analyses, we obtained genotype information for variants across populations by mapping short-read whole-genome sequencing (WGS) data to the pangenome graph (Fig. 2B). We refer to this data set as pangenome+short-read WGS (PGsr). We focused all further analyses on the study of deletions. To identify variants phenotypically relevant to cavefish, we isolated deletions that were statistically associated with cave-morph fish in each lineage. From these deletions, we identified those that were consistently called across different SV calling methods, not present in any sampled surface population (Río Choy, Mante, Jalpan, and Rascón), and high frequency (deletion allele frequency ≥0.8) in cave populations (Fig. 2C). We refer to this final data set of most-filtered deletions as “focal deletions.” In some tests below, we compare genes with focal deletions to “control genes,” which contain no SV insertions or deletions within the gene structure in any population in the pangenome. For full study details, please see the Methods section.

Figure 2.

Graphical overview of study design. (A) Graphical representation of pangenome graph construction using surface fish genome as a reference frame, as well as examples of insertion and deletions represented in the pangenome graph. (B) SV calling methods from linked-read sequencing (left) and mapping short-read sequencing to pangenome graph (right). (C) Integration of SV calls from different calling methods and filtering of SV deletion variants to generate a set of focal deletions used for downstream analyses.

Open in new tabDownload PowerPointLink to figure

Pangenome graph construction

Our Mexican tetra pangenome graph, generated from reference genomes of one surface- and three cave-morph individuals, consists of 52,571,860 nodes (sections of DNA sequence) and 72,140,066 edges (connections between nodes indicating that the given sequences are adjacent to one another in at least one genome). During graph construction, addition of the three cave-morph genomes to the pangenome graph contributed ∼400 Mbp of nonsurface sequence, attributable to 176 Mbp in Tinaja, 131 Mbp in Molino, and 94 Mbp in Pachón. We did not reach a plateau in new nonreference sequence nodes, indicating that inclusion of additional genomes would continue to grow the pangenome graph. Of the four reference genomes included, the pangenome graph used the Mexican tetra Río Choy surface fish reference (Warren et al. 2024) as a reference frame. Therefore, all results describe variation occurring in a population or individual relative to Río Choy, and genome positions are given in terms of the Río Choy coordinates.

Characteristics of pangenome SVs

From the pangenome graph, we identified variant sites: positions in the surface reference at which the cave reference genomes differ. For top-level indel variants (i.e., primary variants that are not nested within other variants, LV = 0) on chromosomes specifically, we identified a total of 1,100,450 sites with 1 bp indels, 1,359,534 sites with small insertions (<50 bp), 1,489,289 sites with small deletions (<50 bp), 236,786 sites with SV insertions (sequence gains 50 bp–100,000 bp), and 216,339 sites with SV deletions (sequence losses 50 bp–100,000 bp). We note that some sites affected by SVs are multiallelic, with unique SV alleles between caves. Across SV insertion sites, we found 338,631 unique insertion alleles, and across SV deletion sites, we found 255,331 unique alleles (Fig. 3A). SV alleles were more often shared within lineages than between lineages (lineage 1 = Molino, lineage 2 = Pachón and Tinaja). The majority (80.31%) of SV insertions and deletions we identified in the pangenome are ≤1 kb in length (Supplemental Fig. 1A). Alleles for structural insertions amount to an average of 138 Mbp per cave morph, whereas alleles for structural deletions amount to an average of 106 Mbp per cave morph (Fig. 3B). Our finding that insertions account for more base pairs compared with deletions is in line with findings from subterranean isopods, decapods, and mollusks, which have increased genome sizes relative to their surface counterparts (Lefébure et al. 2017). We find more total base pairs affected by SVs in the lineage 2 populations, Pachón and Tinaja, likely because the surface reference is from a Río Choy individual that belongs to lineage 1 and is therefore more closely related to the Molino cave population.

Figure 3.

Structural variation identified between Mexican tetra reference genome assemblies. (A) Comparison of pangenome SV deletion alleles (left) and SV insertion alleles (right) by presence in each cave population (lineage 1 = Molino, lineage 2 = Pachón and Tinaja). Overlapping sections represent shared variant alleles between populations within the pangenome. Variants called relative to the Río Choy reference genome. (B) Total number of base pairs affected by pangenome SV insertions (positive values) and SV deletions (negative values) in each cave population. (C) Pairwise comparison of duplications (lower triangle matrix) and translocations (upper triangle matrix) in each of the four A. mexicanus reference assemblies. Variants called from pairwise whole-genome alignments between the two assemblies indicated by each cell (row × column). Values in each cell represent the total number of structural differences detected between the specific pair of assemblies. Duplication counts include both tandem and interspersed duplications but do not include tandem repeats. Darker cells indicate fewer variants identified between two assemblies and, therefore, more structurally similar genomes.

Open in new tabDownload PowerPointLink to figure

Additionally, we identified duplications and translocations between all pairwise combinations of the four A. mexicanus reference assemblies in the pangenome (Río Choy, Molino, Pachón, and Tinaja) (Fig. 3C). Each pairwise comparison resulted in about 1000–2300 variant calls for both duplications and translocations, representing far fewer variants compared with insertions and deletions. This finding is consistent with other within-species studies across taxa showing that indels outnumber translocations and duplications by more than two orders of magnitude (Long et al. 2018; Dhakal et al. 2024; Kaur et al. 2024). We found fewer duplications and translocations in the comparison of cavefish from the same lineage (Pachón and Tinaja) than in comparisons between cavefish lineages (Pachón or Tinaja compared with Molino). This finding further indicates that the Pachón and Tinaja genomes are more structurally similar to one another than to Molino.

Across the 26,736 Mexican tetra protein-coding genes, 18,016 (67%) have SV deletions and 17,715 (66%) have SV insertions in Molino, Pachón, or Tinaja within intron, exon, or UTR regions. Longer genes offer more mutational targets and are therefore more likely to contain variation (Lopes et al. 2021), a trend that is well supported in Mexican tetras (Moran et al. 2023; Roback et al. 2025). Likewise, we found that genes that contain no SV insertions or deletions are shortest in total base pairs (median = 4847 bp), followed by genes with SV insertions (median = 6538 bp), then genes with SV deletions (median = 7177 bp), and finally genes with both SV insertions and deletions (median = 23,880 bp; all groups different from one another through pairwise comparisons using Wilcoxon rank-sum test with continuity correction, all P < 0.01, average Cohen's d across all comparisons = 0.33) (Supplemental Table S1A). For genes that contain SV deletions, we find a bimodal distribution in which most genes have ≤25% of their total base pairs deleted, as well as a notable group of genes that are completely deleted (Supplemental Fig. 1B).

Overall, 50.85% of all SV deletions and 41.55% of all SV insertions affect introns. In contrast, relatively few SV insertions and deletions occur in exons. Across the three cave morphs in the pangenome, SV deletions affect exons of 2397 protein-coding genes and SV insertions affect exons of 1332 protein-coding genes. We find that SV insertions and deletions overlap CDSs less frequently than random (permutation test P < 0.005 for both insertions and deletions) (Fig 4A; Supplemental Table S1B) and overlap introns more frequently than random (P = 0.04 for deletions and P < 0.005 for insertions) (Fig. 4B; Supplemental Table S1B). Although both insertions and deletions are overrepresented in introns, the degree of overrepresentation is much greater for insertions. In intergenic regions (i.e., regions of the genome not within the gene structure), there are more insertions and fewer deletions than expected through random accumulation (P < 0.005 for both insertions and deletions) (Supplemental Table S1B). These results suggest that negative selection is likely removing SVs from coding regions, in line with findings in other animals (Liang et al. 2024).

Figure 4.

Insertions and deletions are found nonrandomly within the genome. Results of a permutation test that places identified SVs into the genome at random and determines how frequently the variant overlaps regions of interest. The histograms illustrate the number of intersections between the SV and the given region in 300 runs of the random permutation test, and the dotted line indicates the observed number of intersections. (A) Intersections of deletions (left) and insertions (right) with coding sequence (CDS) regions. (B) Intersections of deletions (left) and insertions (right) with intron regions.

Open in new tabDownload PowerPointLink to figure

Because transposable elements can amplify and insert multiple times, they are inherently repetitive and abundant among insertion events. True to these expectations, SV insertions found in cavefish relative to the Río Choy surface reference are 78.6% repetitive on average across Molino, Pachón, and Tinaja (Supplemental Table S1C). Of this repetitive content, the largest contributors are retroelements (30.7% on average, including 21.3% from LTR elements) and DNA transposons (19.3%). The fact that insertions are largely composed of transposable elements may explain our finding that insertions are overrepresented within introns and intergenic regions (and more so than deletions), as transposable elements show location bias for noncoding regions, especially introns and regions flanking genes (Wei et al. 2016).

Finally, we predicted functional impacts of variants in the pangenome as characterized by SnpEff (Cingolani et al. 2012a). SV deletions represent 30.92% of all high-impact variants, and SV insertions account for 10.65% (Supplemental Fig. 1C; Supplemental Table S1D).

Identification of SV deletions from linked-read sequencing

To identify deletions using a different approach unrelated to the pangenome graph, we used linked-read sequencing from each of two surface populations and seven caves and called deletions relative to the linear Río Choy surface reference genome (Fig. 5; Supplemental Table S2). We obtained fewer deletion calls in Molino, Pachón, and Tinaja through the linked reads compared with the pangenome graph (∼1.7× more deletions were identified in the pangenome graph in Molino and 1.5× more in Pachón and Tinaja). This highlights the sensitivity of SV detection provided by the pangenome graph method compared with use of a linear reference, as also demonstrated by Rice et al. (2023).

Figure 5.

Shared deletions between and within cavefish lineages. Heatmap depicting results of SV deletion calling from linked-read sequencing. Cells on the diagonal represent the number of deletions (50 bp–100 kb) called in each individual relative to the Río Choy reference after filtering for variants passing all Long Ranger quality filters. Darker cells on the diagonal represent a greater number of deletion calls. Cells within the triangle matrix represent the number of shared deletions between individuals from each pair of populations. Deletions were considered shared if they had 90% reciprocal overlap between populations. This cutoff increases the confidence that the variants originate from a shared SV-generating event, rather than repeated origins in a similar location. Darker cells in the triangle matrix represent a greater number of shared deletions between two populations. More deletions are shared within lineages than between lineages, in accordance with the best supported phylogenetic history of these populations.

Open in new tabDownload PowerPointLink to figure

From the linked reads, we identified the fewest deletions in the Jalpan surface individual (75,975 deletions), which is the most closely related population to the Río Choy reference (both lineage 1 surface populations). It is logical that we found the fewest deletions in the population most closely related to the reference used; however, the presence of nearly 76,000 deletions identified between individuals of the same lineage and ecotype suggests that there is substantial variation even between closely related populations with shared environmental pressures. We found the most deletions in the Rascón surface individual (94,619 deletions), which represents the divergence between lineage 1 and lineage 2, the deepest divergence time in A. mexicanus (Herman et al. 2018; Garduño-Sánchez et al. 2023; Moran et al. 2023). This finding is in line with the documented genetic diversity existing between surface populations belonging to different lineages. For example, the observed value of D_XY (a measure of absolute nucleotide divergence) between Río Choy and Rascón is greater than any comparison of either surface population to a cave population (Herman et al. 2018). Although most population genetic measures of diversity such as D_XY are derived from SNP data, our results highlight that there is also a high degree of structural diversity between surface populations from different lineages.

As with our finding for pangenome SVs, we found that linked-read deletions are most often shared between individuals of the same cave lineage rather than between cave lineages (Fig. 5), supporting current evidence that at least two cave lineages diverged from at least two separate surface lineages (Ornelas-García et al. 2008; Bradic et al. 2012; Herman et al. 2018; Moran et al. 2023). Additionally, the pattern of variation being shared within a lineage more commonly than between lineages has also been found in examination of other types of variants in this system (Roback et al. 2025). Overall, the number of deletions detected in each population relative to Río Choy and the comparative number of deletions shared between versus within lineages follow our expectations based on the documented evolutionary history of A. mexicanus.

Identification of high-confidence cave-associated focal deletions

To further explore deletions segregating in Mexican tetra populations, we used the pangenome as a reference for the alignment and genotyping of 107 short-read WGS samples from nine total populations (Fig. 1D; Supplemental Table S2). This approach allows SVs present in the pangenome to be confidently and accurately genotyped in many populations and individuals using short reads, which would not be possible with traditional linear reference alignments (Hickey et al. 2020; Sirén et al. 2021). We note, however, that this method does not identify additional variants that may exist in the WGS data, and therefore, there are likely many false negatives. In other words, our data set of genotyped variants is not an exhaustive list of all variants that exist in A. mexicanus populations.

Our genotyping yielded 279,995 PGsr deletions. To identify deletions statistically associated with cave-morph fish, we identified PGsr deletions present in each lineage and conducted lineage-specific runs of SnpSift CaseControl (Cingolani et al. 2012b) comparing lineage 1 caves to lineage 1 surface populations (Río Choy and Mante) and lineage 2 caves to the lineage 2 surface population (Rascón). We used the allelic model and Fisher's exact test with a Bonferroni–Holmes correction to assess differences in deletion allele frequencies between cave populations (i.e., cases) and surface populations (i.e., controls). We identified 23,032 deletions in lineage 1 and 6606 deletions in lineage 2 that were statistically associated with cave morphs (adjusted P < 0.05). The analysis in lineage 2 likely returned fewer deletions owing to artifacts of SnpSift CaseControl not because lineage 2 truly contains fewer cave-associated deletions (for full commentary, see Methods).

In line with best practices, we used multiple SV calling methods and retained SVs with consensus between linked-read and pangenome graph mapping approaches to increase confidence in our PGsr cave-associated deletion calls (Ho et al. 2020). Taking a consensus approach that combines both pangenome graph and linear reference-based methods significantly improves the precision of rare SV calls, which are more prone to false positives (Groza et al. 2024). Of the deletions associated with cave morphs by SnpSift CaseControl, 14,050 deletions in lineage 1 and 5246 deletions in lineage 2 were supported by linked-read deletion calls. Repeated sourcing of alleles from standing genetic variation is a prominent and important source of genetic convergence (Barrett and Schluter 2008) and is well documented in systems such as three-spined stickleback fishes, in which alleles associated with freshwater adaptation are frequently present in marine populations at low frequency (Colosimo et al. 2005; Jones et al. 2012). We identified whether our linked-read-supported, cave-associated deletions were present as a minor allele in related surface populations by identifying those with a surface allele frequency ≤0.2 and greater than zero. In lineage 1, 3437 cave-associated deletions (24%) are present as a minor allele in Río Choy and/or Mante based on our PGsr genotyping (mean deletion allele frequency; Río Choy = 0.11, Mante = 0.10). In lineage 2, 3222 cave-associated deletions (61%) are present as a minor allele in Rascón (mean deletion allele frequency = 0.09). The smaller proportion of lineage 1 deletions present as a minor allele in surface populations is likely an artifact of deletions being called against the Río Choy surface reference; if a deletion is present in the Río Choy population and represented in the Río Choy reference genome, the deletion will not be called in our data set because deletions are identified when sequence is missing relative to the Río Choy reference.

Although deletions can be repeatedly sourced from standing genetic variation, we specifically sought to explore deletions that have arisen because cave colonization and compare putatively independent deletion events between cave lineages. Therefore, we combined our cave-associated deletions from each lineage and isolated cave-specific deletions by removing deletions present at any frequency in any surface population (Río Choy, Mante, Jalpan, and Rascón; including all individuals from the linked-read and PGsr data sets). Finally, to identify deletions that have risen in frequency toward fixation, we required a deletion allele frequency of ≥0.8 in at least one cave population. As a result, we obtained a final set of 2431 unique deletions hereafter referred to as “focal deletions” (Supplemental Table S3).

Deletions identified as significantly associated with cave-morph fish in one lineage are not necessarily absent from the other; rather, they may exist at lower frequencies or lack sufficient genotype coverage to achieve statistical significance. Among the focal deletions in lineage 1, 73% are unique to that lineage (i.e., not present at any frequency in lineage 2), and 79% of focal deletions in lineage 2 are unique to lineage 2 (Supplemental Table S4A). Deletions unique to one lineage are almost always present in multiple caves within the lineage (99% of unique lineage 1 deletions are present in two or more lineage 1 caves and 97% of unique lineage 2 deletions are present in two or more lineage 2 caves). Approximately 26% of focal deletions are shared between lineages. These shared deletions may represent introgression between cave populations, reflect standing genetic variation that was undetected in our sampling of surface fish, or have arisen through similar de novo mutations.

In sum, we identify a set of 2431 focal deletions that are strong candidates to contribute to cave-derived phenotypes. These deletions are significantly associated with cave-morph fish according to SnpSift, are supported by both the PGsr data set and linked-read data set (i.e., are consistently identified across different SV calling methods), are absent in all sampled surface fish, and are high frequency (deletion allele frequency is ≥0.8) in one or more cave populations.

Characteristics of focal deletions

Like the set of all deletions identified in the pangenome, the majority (96%) of focal deletions (n = 2431) are ≤1 kb (mean = 202 bp, median = 80 bp) (Supplemental Table S4B) and occur in CDS regions significantly less often than would be expected through random accumulation (only 14 focal deletions are within exons, permutation test P < 0.003) (for more detail, see Supplemental Table S4C). Again, this suggests that overarching purifying selection removes deletions that occur in coding regions. Focal deletions intersect introns and intergenic regions no more often than expected under random accumulation (permutation test, P = 0.303 and P = 0.376, respectively) (Supplemental Table S4C). For focal deletions that occur anywhere within the structure of protein-coding genes (n = 1530), we do not find evidence that the distribution of the starting position of the deletion within the gene length differs from a uniform distribution (chi-squared comparison, X = 9.098, df = 9, P = 0.4283) (Supplemental Table S4D). Thus, focal deletions are not concentrated at certain positions along the gene structure.

The position of deletions in the genome may be biased by the sequence context in that region. For example, repeated deletions within the pelvic enhancer region of stickleback fish are thought to be stimulated by stretches of TG sequence repeats that are highly mutagenic (Xie et al. 2019). Because TG repeats can promote deletion mutations, we explored whether any of our focal deletions occurred at TG repeat sites. We found six deletions in lineage 1 and four deletions in lineage 2 that contained a TG repeat region of ≥20 bp (TG₁₀₊) (Supplemental Table S4E); however, none of the focal deletions at TG repeat sites are shared between lineages or occur within an overlapping region. We note one deletion of interest that occurs at a 54 bp TG repeat (TG₂₇), a 1410 bp deletion fixed in the Molino and Escondido caves that begins 5413 bp upstream of ptprdb, a gene associated with diabetes risk and insulin resistance in humans (Chang et al. 2012).

We identified only 13 shared focal deletions that are associated with the cave morph in both lineages (i.e., only 13 focal deletions are significant by SnpSift in both lineages) (Supplemental Table S5A) but find 67 instances of focal deletions in the same protein-coding gene but at different positions in different lineages (Fig. 6A; Supplemental Table S5B). Previous research has indicated that genes with large mutational target size are more likely to undergo repeated evolution (Gompel and Prud'homme 2009; Moran et al. 2023); correspondingly, we find that genes that repeatedly incurred unique deletions in each lineage are three times longer (in terms of total base pairs) than the rest of genes that bear focal deletions (median length for genes with repeated deletions = 231,341 bp, median length for all other genes impacted by focal deletions = 74,387 bp, pairwise comparisons using Wilcoxon rank-sum test with continuity correction P < 0.0001). This suggests that genomic biases such as length influence which genes are affected repeatedly in independent bouts of evolution.

Figure 6.

Repeated deletions and gene reuse between independently evolved cavefish lineages. (A) Details of the different deletions found in genes pcsk2 and ksr2 within each cavefish lineage. Introns have been scaled to 100 bp to visualize the entire gene structure in one image. Positions of gene start, deletion variant, and gene end are given in the Río Choy reference genome coordinates. Tube maps illustrate the deletion call within the pangenome graph, and pie charts show the genotype calls and population-level deletion allele frequency found through mapping short-read sequencing to the pangenome graph. (B) Putative relationships between genes containing deletions in both lineages and cave-derived phenotypes. Functional relevance for each gene was identified by exploring documented functional impact of the gene ortholog in zebrafish, mice, and humans, leveraging the databases Zfin, MGI, and OMIM and reviews of current literature (Supplemental Table S5C). Each bubble represents a suite of cavefish phenotypes, which are illustrated and compared with surface fish phenotypes in the corresponding boxes. The genes listed in each bubble have evidence of functional impact related to the cavefish phenotype indicated by the bubble. Examples of enriched pathways identified in ingenuity pathway analysis of focal deletions pertinent to the phenotypes are given in dashed boxes (for full list, see Supplemental Table S7A).

Open in new tabDownload PowerPointLink to figure

Putative functional impacts of focal deletions and association with cave-derived traits

We explored the potential functional impacts for genes containing focal deletions in both lineages through inspection of documented phenotypes in gene orthologs from zebrafish, humans, and mice (Supplemental Table S5C). Many genes containing deletions in both lineages have documented phenotypic impacts related to cave-derived phenotypes (Fig. 6B).

Additionally, we identified biological pathways overrepresented among the set of genes affected by focal deletions in each lineage through ingenuity pathway analysis (IPA) of human orthologs of Mexican tetra genes. We found 224 enriched pathways in lineage 1 and six enriched pathways in lineage 2 (Benjamini–Hochberg adjusted P-values < 0.05) (Supplemental Table S6A). The smaller number of enriched pathways in lineage 2 reflects the lower number of cave-associated deletions identified in this lineage, resulting in fewer genes input for this analysis. All pathways enriched in lineage 2 were also enriched in lineage 1: RHO GTPase cycle, autism signaling pathway, synaptic adhesion-like molecules, ROBO SLIT signaling pathway, GABAergic receptor signaling pathway (enhanced), and neurexins and neuroligins. Despite this overall concordance in enrichment, mostly different genes in the same pathway incurred deletions between the two lineages (Supplemental Table S6B), indicating a pattern of convergence at the level of the pathway rather than of convergence at the gene level. To consider cases of deletions affecting a gene's regulatory sequence, we also ran pathway enrichment for genes with focal deletions within 5 kb upstream (four enriched pathways in lineage 1, one enriched pathway in lineage 2, no shared pathways) (Supplemental Table S6C) of and 5 kb downstream (no enriched pathways) from their coding sequence. Many enriched pathways correspond to the production of cavefish phenotypes including metabolism, eyes, heart, and behavior (when considering pathways significant in either lineage) (Fig. 6B).

Because most available RNA-seq data are derived from Pachón cavefish and Río Choy surface fish, we investigated the impact of focal deletions on Pachón gene expression. WGS genotype data and RNA-seq are not available for the same individuals, so we used only fixed focal deletions found in Pachón (obtained by filtering to only include deletions that have F_ST = 1 between Pachón and Río Choy surface, n = 155 deletions). We then identified genes containing a fixed deletion in cumulative exons and introns for the gene or within 5000 bp upstream of and downstream from the gene (n = 134 genes). Correspondingly, for this test we restricted our set of control genes (i.e., genes that contain no SV insertions or deletions in the gene structure) to also exclude genes with SVs within 5000 bp upstream of and downstream from the gene. We found that having a fixed focal deletion does not affect the likelihood that a gene is differentially expressed between Pachón cavefish and Río Choy surface fish in the eye at 54 h postfertilization (Gore et al. 2018) or in liver tissues (Krishnan et al. 2020) compared with control genes (Fisher's exact test, P > 0.05) (Supplemental Table S7A). Because surface Mexican tetras can regenerate heart muscle after injury whereas cavefish retain permanent scarring, we also explored differential expression in heart tissues at 3, 7, and 14 days after injury (Stockdale et al. 2018). We found that genes with fixed focal deletions are more likely to be differentially expressed at 7 days after injury (Fisher's exact test, P < 0.05) but not at 3 or 14 days after injury (Supplemental Table S7A). Additionally, we used a hypergeometric test that assesses whether a disproportionate number of differentially expressed genes contain fixed structural deletions in Pachón and obtain the same pattern of significance (P < 0.05 for hearts 7 days after injury, P > 0.05 for all other tissues) (Supplemental Table S7A). We might have identified differential expression in deletion genes in only the heart samples because examining RNA-seq at multiple time points can better capture dynamic expression patterns that are missed in single time point analyses like the liver and eye samples.

Of the 134 genes that have a fixed focal deletion in Pachón, we identified 20 genes with fixed deletions that are differentially expressed in the eye, 30 in the liver, and 33 in the heart across all time points, as well as 14 genes with fixed deletions that are differentially expressed in multiple tissues (Supplemental Table S7A). In accordance with the eye loss phenotype seen in cavefish, 85% of differentially expressed deletion genes in the eye are downregulated in the cave (Supplemental Table S7A). These results suggest that a subset of deletions likely have impacts on expression and therefore cavefish phenotypes, and these genes merit further examination.

Most QTLs identified in Mexican tetras are derived from crosses between Pachón cavefish and Río Choy surface fish (Wiese et al. 2024). Therefore, we again used fixed focal deletions in Pachón (n = 155 deletions) to identify whether deletions were associated with regions of the genome thought to contribute to cave-derived traits. Using previously established QTLs in Pachón mapped to the Río Choy surface reference by Wiese et al. (2024), we found that of the 155 fixed focal deletions in Pachón, 98 are found within a QTL region. We then ran a permutation test and identified that fixed deletions are found within QTLs more frequently than expected by chance and significantly more frequently than expected based on the frequency of Pachón deletions in the pangenome within QTLs (permutation test P < 0.005; Fisher's exact test P < 0.0005) (Fig. 7A; Supplemental Table S7B). These results show that focal deletions are enriched within genomic regions associated with cave-derived phenotypes such as reduction of pigmentation and expansion of taste buds (Protas et al. 2007), supporting their potential role in the genetic basis of these traits.

Figure 7.

Deletions are enriched within QTL for cavefish traits and genes with evidence of selection. (A) Overlap between deletions and QTLs from previous studies of the Pachón cavefish population (see Supplemental Table S8B). The histogram illustrates the number of deletions fixed in Pachón (F_ST = 1 between Pachón and Río Choy, n = 155 deletions) that were found within Pachón QTL regions in each run of a permutation test that produces a simulated VCF with deletions placed into the genome at random. The dotted line indicates the true number of deletions within QTL regions (98 deletions). (B) Proportion of genes that contain deletions and selective sweeps in three unique gene sets: (PG) Genes containing deletions that are present in the pangenome but did not have a statistical association with cave-morph fish (n = 16,617), (FOCAL) genes containing focal deletions that have a significant association with the cave-morph and are not present in surface fish (n = 1260), and (REUSE) genes that contain a unique focal deletion in each cavefish lineage (n = 67). Genes must have a selective sweep and a deletion in the same cave population. Genes are placed into a category based on their highest level of association (e.g., a gene containing both a PG deletion and a FOCAL deletion is put in the FOCAL category). (C) Proportion of control genes (i.e., genes that contain no insertion or deletion SVs, n = 5265) that contain one or more selective sweeps in each cave population.

Open in new tabDownload PowerPointLink to figure

Genes with cave-associated deletions are more likely to be under selection

To assess whether deletions associated with cave phenotypes are more likely to be under selection, we identified whether genes containing deletions also showed evidence of selective sweeps in the specific cave population(s) in which the deletion is present. Evidence of selection was identified using diploS/HIC (Kern and Schrider 2018), a deep convolutional neural network approach that identifies sweeps from patterns of genetic variation. We identified the presence of sweeps in 5 kb windows across the genome within the Molino, Pachón, and Tinaja cave populations. In line with previous analyses (Moran et al. 2023; Roback et al. 2025), genes that contained one or more windows with a hard or soft sweep call were considered to have evidence of selection in the population(s) with the sweep call (i.e., a gene may have evidence of selection in no cave populations, one population, or multiple populations). Within these genes, we required that the deletion and selective sweep were present in the same population.

After identifying the presence of sweep windows within all genes, we compared three deletion gene sets: PG genes (genes with deletions in the pangenome that are not statistically associated with cave-morph fish and may also be present in surface fish), focal genes (genes with deletions that are associated with cave morphs and are not present in any sampled surface fish, consistent with above usage of this term), and reuse genes (genes that contain focal deletions at different positions in each cavefish lineage). We found that the proportion of genes with evidence of selection increased significantly with the level of cave-association: PG genes had the lowest proportion under selection (53%) followed by focal genes (62%), and reuse genes had the highest proportion under selection (99%) (pairwise comparison of each gene set using chi-squared test, all P < 0.0001) (Fig. 7B; Supplemental Table S5D). Therefore, genes with cave-associated deletions are significantly more likely to have population-specific selective sweeps than are genes with deletions not associated with the cave-morph phenotype. Further, genes that harbor unique deletions in each cave lineage are more likely to have evidence of population-specific selection, supporting their potential role in the repeated evolution of cave-adapted traits. We note however, that our analysis does not definitively show that there has been selection for the deletion variant specifically, only that a deletion and selective sweep co-occur within a gene in a population.

To provide an additional comparison, we explored how often control genes with no insertion or deletion SVs had evidence of selection in each population. Only 8%, 9%, and 15% of control genes had evidence of selection in Molino, Pachón, and Tinaja, respectively (Fig. 7C; Supplemental Table S5D), substantially less than the proportion of genes with deletions under selection. This result further supports the hypothesis that structural variation, particularly deletions, contributes to adaptation in cavefish.

Pangenome construction and processing

We built a pangenome from autosomal sequences of the surface-morph and three cave-morph A. mexicanus assemblies AstMex3_surface (Río Choy), AstMex3_Molino (Molino cave), AstMex3_Tinaja (Tinaja cave), and AMEX_1.1 (Pachon cave) (Imarazene et al. 2021; Warren et al. 2024). We ran the Minigraph-Cactus pipeline (Hickey et al. 2024) using the cactus v2.4.2 Docker image as described previously (Rice et al. 2023). We specified the Río Choy surface assembly as the reference to aid in the identification of variation existing in cave-adapted fish relative to those on the surface. The graph was output in multiple formats, including GFAv1.1, hal, vg, and VCF.

As default in cactus 2.4.2, the VCF was postprocessed with vcfbub, which exposes the snarl tree to annotate top-level (LV = 0) and nested variants (LVs > 0). We subset the graph to only include top-level variants (“vcfbub -l 0”) for all subsequent analyses. Although an updated version of Minigraph-Cactus left-aligns this VCF with BCFtools norm (Danecek et al. 2021), the version used in this study (v2.4.2) lacked normalization as a default. We initially postprocessed the top-level VCF with BCFtools norm but, in doing so, observed negative impacts of the left alignment, wherein variants were reassigned new coordinates that made them both incompatible with other formats of the graph (e.g., HAL, GFA, and VG) and incompatible with one another (e.g., an individual cannot have two different alternate genotypes at the same site) (https://github.com/ComparativeGenomicsToolkit/cactus/issues/1557). We chose to use the top-level VCF without left alignment for all subsequent analyses. As a final postprocessing step, we subset the top-level, unnormalized VCF to lines in which at least one individual (Molino, Pachón, or Tinaja) had an alternate genotype. Given that our input assemblies are haploid compressed genome assemblies, each “genotype” in the pangenome is a single allele (e.g., “1” instead of “1/1”).

Pangenome genotyping

In the postprocessed VCF, variant types were categorized based on the length of the alternate (ALT) allele and the length of the reference allele (REF) after splitting multiallelic sites. SVs were defined as those ≥50 bp and up to 100 kbp, in accordance with other pangenome approaches (e.g., Fang and Edwards 2024). To set equal thresholds across insertions and deletions, we categorized the structural insertions >100 kbp as “large structural insertions.” Similarly, 1 bp indels, small insertions (2 bp–49 bp), and small deletions (2 bp–49 bp), were defined through inspection and subtraction of the ALT and REF alleles. We defined 1 bp indels as 1 bp sequence gains or sequence losses that do not, otherwise, meet the criteria of being a multinucleotide polymorphism.

Pangenome summary statistics

Variants on unplaced scaffolds were excluded from pangenome summary statistics reported in the paper. Throughout the paper, a gene is defined in accordance with the NCBI definition used for GTF files, in which “gene” describes coordinates encompassing all isoforms of a gene (the union of all introns and exons). For this study, we only analyzed protein-coding genes. Statistics regarding the relationship of SVs to different genomic features, other than repetitive elements, were calculated based on the annotation of the pangenome VCF (see section “Annotating genomic features of SVs”). Calculations of gene length, proportion of gene deleted, and position of deletions within a gene were all calculated based on the gene coordinates specified in the GTF.

Identification of repetitive regions within insertions

We used RepeatModeler 2.0.4 (Flynn et al. 2020) with the -LTRStruct option to create libraries of transposable elements for the surface and each of the cave-morph genome assemblies. For each assembled cave morph, we subset the VCF to contain only insertions present in that genome assembly. We converted column 5 of the resulting VCF to a FASTA file and used this, along with the repeat-families.fa output of RepeatModeler as described above, as input to RepeatMasker v4.1.7-p1 (Tarailo-Graovac and Chen 2009) with default parameters.

Identification of duplications and translocations

We used the Synteny and Rearrangement Identifier SyRi v1.6 (Goel et al. 2019) to identify duplications and translocations between all pairwise combinations of the four A. mexicanus reference assemblies in the pangenome (Río Choy, Molino, Pachón, and Tinaja). SyRI was run with default parameters, and DUP and TRANS parent-level records were extracted from the output and counted for each pairwise comparison.

Annotating genomic features of SVs

To identify the impacts of SVs in A. mexicanus, we annotated the pangenome VCF using SnpEff 5.2a (Cingolani et al. 2012a) with informed breakpoints. The SnpEff database was built from the AstMex3_surface reference using the NCBI GTF. Custom annotations were provided for population-specific sweeps, as well as QTLs from Wiese et al. (2024) using the –interval option. Additionally, we added the tag “CDSID:” using coordinates from AstMex3_surface NCBI GTF with the –interval option.

Short-read WGS genotyping of surface and cave populations

We mapped and genotyped short-read surface- and cave-morph population data (for full information and sequence accessions, see Supplemental Table S2) to the pangenome graph described in this study using vg giraffe with default options (Sirén et al. 2021). GAM format files were surjected to BAM format with surface serving as the reference genome using the command “vg surject” with default options. Only SVs within 50 bp–100 kbp were retained for downstream analyses.

Linked reads generation and SV genotyping

To identify deletions using a different approach unrelated to the pangenome graph, we used 10x Genomics linked-read sequencing of nine individuals from two surface populations (Jalpan and Rascón) and seven caves (Escondido, Molino, Japones, Pachón, Tigre, Tinaja, and Vásquez). All processing of the linked-read data and SV calling was completed using Long Ranger (v2.2.2) (Zheng et al. 2016). We generated FASTQ files from the Illumina binary base call sequence output by running longranger mkfastq and created a 10x compatible reference from the A. mexicanus surface fish reference genome (GCF_023375975.1) using the longranger mkref command. We then ran Long Ranger in whole-genome mode (wgs). The longranger wgs command performs alignment of sequences to the reference genome, deduplication, and filtering and uses the Chromium molecular barcodes to call and phase SNPs, indels, and SVs. From the phased SVs called by Long Ranger, we retained only variants 50 bp–100 kb in length that passed all filters tested in the pipeline. These filters use barcode and phasing information to improve the quality of variant calls and reduce false positives (for a description of the filters in the Long Ranger pipeline, see https://support.10xgenomics.com/genome-exome/software/pipelines/latest/output/vcf). To determine the number of shared deletions between each population pair, we used BEDTools intersect (Quinlan 2014) to identify and count deletion SV calls with ≥90% reciprocal overlap. This level of overlap increases the likelihood that deletions being counted as shared are indeed comparable variants rather than unique deletions in similar positions.

Identifying a high-confidence set of cave candidate deletions

We integrated short-read genotyping with linked-read genotyping, employing stringent filtering criteria and a statistical test to identify SVs associated with cave-morph fish. First, we identified deletions in the PGsr data set statistically associated with cave-morph fish using SnpSift (Cingolani et al. 2012b). The allelic model implemented in SnpSift's CaseControl analysis (v5.2) utilizes the Fisher's exact test to assess statistically significant associations between genotype and phenotype. This test is effective in identifying variants exhibiting a high frequency of ALT alleles in case samples with a high frequency of REFs in control samples. We applied this test in a phylogenetically informed framework, treating cave samples as the “case” and surface samples as the “control.” To account for the independent origins of the cave phenotype in lineage 1 and lineage 2, we performed separate analyses for each lineage using the postprocessed PGsr data. Only SV deletions were included as input, with multiallelic sites joined. We adjusted P-values with Holm–Bonferroni correction using a custom code written with pysam (v. 0.22.1) (https://github.com/pysam-developers/pysam) and the multipletests() function in statsmodels 0.14.0 (Seabold and Perktold 2010), setting the significance threshold to adjusted P = 0.05.

We likely identified fewer significant deletions in the lineage 2 SnpSift CaseControl analysis owing to (1) reduced statistical power resulting from fewer surface individuals in lineage 2 (lineage 1 = 19 controls [i.e., surface], 31 cases [i.e., cave]; lineage 2 = 12 controls, 45 cases) and (2) an artifact of using the Río Choy surface reference wherein a multiallelic site with an alternate Rascón allele (the lineage 2 surface population/control) would be excluded by SnpSift owing to the absence of reference genotype in the control samples. We also note that variants deemed significant in the SnpSift CaseControl analysis are only those for which there was high genotype coverage resulting from mapping the short-read WGS to the pangenome. Therefore, we can confidently say these deletions are associated with cave-morph fish, however, there are likely additional deletions that are potentially biologically relevant but do not meet our threshold for significance. For instance, deletions in oca2 are well documented in cavefish (first identified by Protas et al. 2006) and present in our PGsr deletion data set but are not identified as significant through SnpSift CaseControl owing to reduced genotype coverage in our samples.

In investigations of structural variation, the confidence of SV calls is increased using multiple SV calling methods and retaining calls with consensus between programs (Ho et al. 2020). Generally, a 50% reciprocal overlap between callers is used to determine high-confidence SV calls (Sudmant et al. 2015; Audano et al. 2019; Kosugi et al. 2019; Ho et al. 2020). We generated a set of high-confidence SV deletions from PGsr genotyping by requiring a ≥50% reciprocal overlap with a deletion call from Long Ranger in a linked-read cave sample from the same lineage.

Next, we isolated cave-specific deletions in a two-step filtering process. We first removed deletions with a ≥50% reciprocal overlap with deletions called in the linked-read surface samples (Jalpan for lineage 1 and Rascón for lineage 2). Second, we calculated the population-specific frequency of our high-confidence deletions from PGsr genotyping using a custom script and required an ALT allele frequency of zero in all sampled surface populations (lineage 1: Río Choy and Mante; lineage 2: Rascón).

Detection of genome-wide selective sweeps

To find genomic regions possibly associated with adaptation, we identified regions of the genome with evidence of selective sweeps using diploS/HIC (Kern and Schrider 2018). This method combines a convoluted neural network approach and calculation of population genetic statistics from sequence data to infer selection. diploS/HIC requires empirical population genetic information in VCF format. We generated this input using the same srWGS data used to genotype pangenome SVs, focusing on the populations with the most individuals sequenced: two surface populations (Río Choy and Rascón) and three cave populations (Molino, Pachón, and Tinaja) (for sequence accessions, see Supplemental Table S2). We first performed genotyping against the Río Choy surface reference (Warren et al. 2024) using the Genome Analysis Toolkit (GATK) (Van der Auwera and O'Connor 2020) v4.4.0 GenotypeGVCFs tool to produce a VCF file for each population. We then subset these VCFs to only include SNPs and removed any SNPs that occurred within indels or a 3 bp buffer region on either side. We followed GATK best practices (https://gatk.broadinstitute.org/hc/en-us/articles/360035890471-Hard-filtering-germline-short-variants) to hard-filter SNPs and removed SNPs that occurred within repetitive regions using the NCBI WindowMasker files.

Next, we simulated population genetic data under neutral, soft, and hard sweeps in discoal (Kern and Schrider 2016), using the same demographic parameters as in a past run of diploS/HIC in Mexican tetra (Moran et al. 2023). We created separate simulated data sets for surface and cave to account for past population size changes and reduce the potential for false positives resulting from population bottlenecks that can cause neutrally evolving genomic regions to exhibit low genetic variation mimicking a selective sweep. For example, the cave demographic model includes a harsh bottleneck in line with the timing of cave colonization that dramatically reduces the effective population size. We then generated summary statistics and feature vectors for the simulated surface and cave data sets, which we used to train and test the convoluted neural network and produce a model for cave and surface. Finally, we supplied our VCFs containing real population genetic data to the trained diploS/HIC model, which generated another set of summary statistics and feature vectors that were used to predict selection in 5 kb windows across the genome. Each window was classified as being neutral, as containing hard sweeps or soft sweeps, or as being linked to a region that experienced a hard or soft sweep. Windows with no call were removed. Although diploS/HIC is a powerful approach to detect regions of the genome containing selective sweeps versus those that do not, it is not as reliable at distinguishing hard sweeps from soft sweeps (Kern and Schrider 2018). For this reason, we grouped hard and soft sweeps as evidence of selection and do not distinguish between the two. To further increase the confidence of our selection calls, we then swapped the underlying demographic models (running the cave population genetic data set with the surface demographic model and vice versa) and retained only windows with consistent calls between models.

To determine whether a gene had evidence of selection, we identified diploS/HIC windows within genes using the NCBI AstMex3_surface GTF gene coordinates. Following the procedure from previous analyses (Moran et al. 2023; Roback et al. 2025), genes that contained one or more windows with a hard or soft sweep call were considered to have evidence of selection in the population(s) with the call. Given that SV deletions span many base pairs and result in removed sequence context, it is difficult to accurately identify whether deletions are found exactly within a selective sweep. However, the population-specific nature of this analysis enabled us to identify genes that harbor a deletion in a certain population and simultaneously have evidence of selection in the same population. We then compared three sets of genes containing deletions: PG genes (genes with deletions in the pangenome that are not statistically associated with cave-morph fish; these genes may be present in surface fish), focal genes (filtered set of genes with deletions that are associated with cave morphs and are not present in any sampled surface fish), and reuse genes (genes that contain focal deletions at different positions in each cavefish lineage and are not present in our surface fish samples). The proportion of genes with evidence of population-specific selection (i.e., the deletion variant and selective sweep are both present in the gene in the same population) was then compared between each gene set using a pairwise chi-squared test (Supplemental Table S5D). Control genes (those containing no insertion or deletion SVs in any population) were considered to have evidence of selection if they contained one or more windows with a hard or soft sweep call.

Permutation testing of SV prevalence

To assess how the distribution of structural deletions across different genomic features (CDSs, introns, and intergenic regions) differ from chance, we ran a permutation test using a custom script. This test randomizes positions in the VCF, restricted to the lengths of the chromosomes, and generates a user-specified number of permuted VCFs. We generated 300 permutations of pangenome structural deletions. We annotated these with SnpEff following the methods previously described (see section “Annotating genomic features of SVs”). We also ran permutation tests on lineage 1 and lineage 2 focal deletions to test relationships of cave-associated deletions to selective sweeps in cave populations, QTLs, CDSs, introns, and intergenic regions. We used a Fisher's exact test to ascertain how the distribution of focal deletions differs from overall population-specific deletions in the pangenome.

Orthology inference

We used OrthoFinder v3.0.1b1 (Emms and Kelly 2019) to identify human orthologs of genes with focal deletions. Briefly, we used the same methods and species sampling as our prior study (Warren et al. 2024), but whereas the prior study used Ensembl gene predictions as input, here we used NCBI gene predictions. The metadata associated with input species/annotations and a summary of the orthogroup assignments are reported in Supplemental Table S8. Assignment of gene symbols to corresponding stable identifiers was accomplished with scripts reported in our prior study (Warren et al. 2024).

Pathway analysis

We performed signaling pathway overrepresentation tests for statistical significance using IPA (Krämer et al. 2014). For each lineage, we identified three sets of genes: (1) genes with deletions within an intron or exon, (2) genes with deletions within 5 kb upstream of the start codon, and (3) genes with deletions within 5 kb downstream from the stop codon. For each Mexican tetra gene, we identified the human ortholog from our OrthoFinder analysis. Genes with single-copy human orthologs were used as input for IPA because A. mexicanus is unavailable in this software. Each of the three gene sets was run independently for each lineage (six runs in total). Significant pathways were those that had a Benjamini–Hochberg adjusted P-value < 0.05. To not inflate enrichment, we restricted pathway enrichment analysis to single-copy human orthologs (including 1:1 and multi:1 orthologs between A. mexicanus and H. sapiens). This is consistent with methods in functional annotations overall, which rely on orthology relationships beyond 1:1 orthology (see, e.g., Hernández-Plaza et al. 2022). We note that because multiple Mexican tetra genes may map to the same human ortholog, the number of genes within a pathway that appear to be impacted in both cavefish lineages is overrepresented (i.e., different cavefish genes with deletions may map to the same human ortholog, which can give the false impression that the same gene is affected in both lineages).

Quantification of cave versus surface gene expression

Previously published RNA-seq data sets were utilized for differential gene expression: wild-type livers (Krishnan et al. 2020), postinjury hearts across three time points (Stockdale et al. 2018), and 54 hours postfertilization (hpf) eyes (Gore et al. 2018). Reads were mapped using the Amex 3.0 NCBI annotation with STAR v2.7.11b (Dobin and Gingeras 2015). We enabled –quantmode to count reads during mapping, which generates identical results to HTseq (Anders et al. 2014). We did not trim reads given that trimming is redundant when mapping involves soft-clipping (Liao and Shi 2020), a default in the STAR aligner. We used DESeq2 to normalize raw counts based on size factors, remove outliers based on Cooke's distance, perform differential expression, and run PCAs on variance stabilized counts (Love et al. 2014). All comparisons were performed with Pachón versus Río Choy, as both populations were available in all data sets. For the livers, eyes, and hearts, the significance threshold was set to an adjusted P-value < 0.05.

Associating candidate deletions with gene expression

Genotyping of the individuals used for RNA-seq is not available. Therefore, to ensure the deletion was present in the Pachón individual used for RNA-seq, we calculated Weir and Cockerham's F_ST using VCFtools on the lineage 2 focal deletions (Danecek et al. 2011). We filtered the data set to only include variants in Pachón that have F_ST = 1 compared with Río Choy. We then identified all genes with a fixed focal deletion that intersects the gene structure, 5000 bp upstream of or downstream from the start and stop site. In addition, we identified a set of control genes that did not have any evidence of structural variation (including both insertions and deletions) within the gene structure or 5000 bp upstream of or downstream from the start and stop site.

We first tested whether the likelihood that a gene is differentially expressed differs between fixed focal deletions and control deletions in each tissue using a Fisher's exact test. Additionally, we tested whether fixed focal genes are differentially expressed more frequently than expected by chance in each tissue using a hypergeometric test with the phyper() function in R (R Core Team 2025), setting lower.tail=FALSE. This calculates whether candidate genes are overrepresented among the genes that are differentially expressed between cave and surface. Differentially expressed genes that are focal genes are treated as “successes” in the test. The significance of the rate of success is based on N, the number of protein-coding genes in the genome; K, the number of differentially expressed genes; M, the number of candidate genes; and x, the number of differentially expressed genes that are focal genes (i.e., “successes”).