The Astyanax mexicanus complex includes two different morphs, a surface and a cave-adapted ecotype, found at three mountain ranges in Northeastern Mexico: Sierra de El Abra, Sierra de Guatemala, and Sierra de la Colmena (Micos). Since their discovery, multiple studies have attempted to characterize the timing and the number of events that gave rise to the evolution of these cave-adapted ecotypes. Here, using RAD-seq and genome-wide sequencing, we assessed the phylogenetic relationships, genetic structure, and gene flow events between the cave and surface Astyanax mexicanus populations, to estimate the time and mode of evolution of the cave-adapted ecotypes. We also evaluated the body shape evolution across different cave lineages using geometric morphometrics to examine the role of phylogenetic signal vs. environmental pressures. We found strong evidence of parallel evolution of cave-adapted ecotypes derived from two separate lineages of surface fish and hypothesize that there may be up to four independent invasions of caves from surface fish. Moreover, a strong congruence between the genetic structure and geographic distribution was observed across the cave populations, with the Sierra de Guatemala the region exhibiting most genetic drift among the cave populations analyzed. Interestingly, we found no evidence of phylogenetic signal in body shape evolution, but we found support for parallel evolution in body shape across independent cave lineages, with cavefish from the Sierra de El Abra reflecting the most divergent morphology relative to surface and other cavefish populations.

To characterize the morphological divergence between surface and cavefish, we performed landmark-based Geometric Morphometrics (GM) analyses on lateral left-hand view of the body shape, for a total of 146 individuals from surface and cave populations. All photographs were transformed to a TPS file using the TPSUtil program. We digitized 12 homologous landmarks and a curve with 15 pseudo-landmarks with the TPSDig2 v. 2.31 (Rohlf, 2015). A Generalized Procrustes Analysis (GPA) was performed with the “gpagen” in Geomorph v. 4.02 (R Core Team, 2019), and with the ‘curve’ argument the sliders were defined with the Procrustes distance criterion to optimize the position of the sliding reference point during the GPA. To avoid a size effect by allometry, the residuals of the regression of the shape on the centroid size (CS) were calculated with the function "procD.lm". Once the residuals were obtained, an allometry-free shape was created, using these residuals of the morphological data. The analyses were made on this allometry-free shape. To assess the degree of variation in body shape between groups, we used principal component analysis using the “gm.prcomp” function.

The archive can be visualized with any text editor; the data was worked in TPS software and analyzed in R, with the 'geomoprh' package.