The genetic regulators of regressive craniofacial morphologies are poorly understood. To shed light on this problem, we examined the freshwater fish Astyanax mexicanus, a species with surface-dwelling and multiple independent eyeless cave-dwelling forms. Changes affecting the skull in cavefish include morphological alterations to the intramembranous circumorbital bones encircling the eye. Many of these modifications, however, have evolved separately from eye loss, such as fragmentation of the third suborbital bone. To understand the genetic architecture of these eye-independent craniofacial alterations, we developed and scored 33 phenotypes in the context of an F2 hybrid mapping pedigree bred from Pachón cavefish and surface fish. We discovered several individuals exhibiting dramatic left–right differences in bone formation, such as extensive fragmentation on the right side only. This observation, along with well-known eye size asymmetry in natural cave-dwelling animals, led us to further evaluate left–right genetic differences for the craniofacial complex. We discovered three phenotypes, inclusive of bone fragmentation and fusion, which demonstrated a directional heritable basis only on one side. Interestingly, the overall areas of affected bones were genetically symmetric. Phenotypic effect plots of these novel craniofacial QTL revealed that cave alleles are associated with abnormal conditions such as bony fusion and fragmentation. Moreover, many linked loci overlapped with other cave-associated traits, suggesting regressive craniofacial changes may evolve through linkage or as antagonistic pleiotropic consequences of cave-associated adaptations. These novel findings illuminate significant craniofacial changes accompanying evolution in complete darkness and reveal complex changes to the skull differentially influenced by genetic changes affecting the left and right sides.
- quantitative trait locus analysis
- regressive phenotypic evolution
- circumorbital bone series
NATURAL model systems, such as Darwin’s finches (Abzhanov et al. 2004, 2006; Mallarino et al. 2011) and African cichlid fishes (Albertson et al. 2003; Albertson and Kocher 2006; Streelman and Albertson 2006), provide living proof that the skull is remarkably labile over evolutionary time. Dramatic morphological changes accompany, or perhaps enable, expansion and exploitation of niches through rapid evolution of adaptive feeding modes (Anyonge and Baker 2006; Cooper et al. 2010). While a great deal is known regarding the developmental origin of cranial bones mediating these constructive changes (Smith 2006; Streelman et al. 2007; Jheon and Schneider 2009), little is known of the genetic basis for craniofacial changes that evolve in the absence of obvious selective pressures.
To explore this phenomenon, we examined alterations to the craniofacial skeleton of the freshwater fish Astyanax mexicanus. This species consists of an extant surface-dwelling form and multiple independently derived cave-dwelling forms (Strecker et al. 2003). As a consequence of invading the subterranean environment millions of years ago (Bradic et al. 2012), cave-dwelling morphs have evolved a series of regressive (e.g., eye loss) and constructive (e.g., increased lateral line sensitivity) phenotypes (Wilkens 1971; Montgomery et al. 2001; Jeffery 2009; Yoshizawa et al. 2010, 2012). Phenotypic loss is believed to arise through genetic drift (Wilkens 1988), direct selection (Klaus et al. 2013), or indirect selection via linkage or pleiotropy (Yamamoto et al. 2009); however, the evolutionary mechanism that drives regressive loss remains unclear (Gross 2012). Our natural model system enables us to investigate the extent to which evolutionary modifications of the craniofacial complex evolve as an indirect consequence of pleiotropy or close physical linkage between causative gene(s) mediating craniofacial traits and other regressive traits, such as eye loss.
Among the most significant changes affecting the skull in cavefish are alterations to circumorbital bone morphology (Alvaréz 1946, 1947). These modifications demonstrate a spectrum of severity, comprising both bone fragmentation and fusion, much of which is endemic to each independent cavefish population (Mitchell et al. 1977). In surface-dwelling forms, each member of the circumorbital series is present as a single, intact bone (Mitchell et al. 1977). In contrast, fragmentation of the suborbital bones [particularly the first suborbital (SO1) and third suborbital (SO3) bones] has been reported for eight wild cavefish populations, including the Pachón cave (Mitchell et al. 1977). Classic studies, which were essentially descriptive accounts, assumed craniofacial changes in cavefish evolved as a secondary consequence of the loss of the eye.
Yamamoto et al. (2003) tested this hypothesis, using lentectomy and intermorphotype grafting to determine the extent to which craniofacial phenotypes were influenced by experimental removal of the eye. Certain traits were affected, including the distance between the nasal and antorbital bones, the inner sectors of the SO3 and supraorbital bones, and position of the SO3 bone relative to the orbit of the eye (Yamamoto et al. 2003). However, other craniofacial traits were not affected by eye loss, such as number of SO3 bony elements, positioning of SO4–6 bones relative to the opercular bone, and opercular bone shape (Yamamoto et al. 2003). Protas et al. (2008) first investigated the genetic basis for craniofacial defects in Astyanax by evaluating variation in SO3 width on the right side of the face only.
In this study, we searched for additional craniofacial traits demonstrating a genetic basis. We scored 33 phenotypes in a hybrid mapping pedigree derived from blind Pachón cave-dwelling and eyed surface-dwelling forms of A. mexicanus (Supporting Information, Table S1). Prior studies using the same pedigree revealed the genetic basis for a number of cave-associated phenotypes, including eye size reduction, pigmentation loss, and chemical sensitivity (Protas et al. 2008). Evaluation of additional phenotypes in the context of this previous work enabled us to determine whether pleiotropic or tightly linked loci are shared between craniofacial and other cave-associated traits. Certain individuals within our pedigree exhibited dramatic left–right differences in bone formation, such as extensive fragmentation on one side only (Figure 1, A–F). This observation, along with the well-documented asymmetry in eye size in natural cave-dwelling populations (Wilkens 2001, 2010; Pouilly and Miranda 2003), motivated us to extend our studies to evaluate left–right genetic symmetry for craniofacial traits. Indeed, we discovered numerous phenotypes with a genetic basis that was detectable on only one side of the cranium. Significant craniofacial changes have accompanied adaptation to the subterranean habitat, and the work herein reveals that complex evolutionary changes to the skull can be influenced by genetic changes affecting the left and right sides differently.
Materials and Methods
We reanalyzed an F2 pedigree of Pachón cavefish × surface fish hybrid individuals (n = 539) with previously collected genotypic data for 164 microsatellite markers (Protas et al. 2007) (Table S2). These individuals represent a full-sibling F2 pedigree from a sibling cross of F1 individuals derived from a mating of a surface-dwelling male and a Pachón cave-dwelling female. We based our linkage map calculations on supplementary information published in Protas et al. (2007), to which we added genotypic information for 11 markers (see Table S2). Each specimen of this pedigree was cleared and stained using Alizarin red to visualize ossified bone (Protas et al. 2008) and stored in sterile glycerol at 4°.
Craniofacial imaging, phenotypic analyses, and scoring
The left and right sides of each individual were imaged using high-resolution microscopy under identical lighting conditions. For each specimen, a high-resolution “montage” image was consolidated from photographs taken at numerous focal positions, which were aligned and flattened into a single, uniformly focused image, as described in Gross and Wilkens (2013). All micrographs were collected using a Leica (Wetzlar, Germany) M205FA stereomicroscope equipped with a DFC310FX camera. Images were collected utilizing the MultiFocus module within the Leica Application Suite (LAS) v3.8 software package. All subsequent measurements were collected from these montage tiff-formatted images.
We collected lateral images at several magnifications (7.81×, 11.0×, 12.5×, 18.0×, 20.0×, and 30.0×) to evaluate 33 craniofacial phenotypes for both the left and right sides of the cranium, including total lateral head area, supraorbital bone area, suborbital bone area (SO2–6), suborbital bone fusion (SO1+2, SO4+5, and SO5+6), suborbital bone fragment number (for bones SO3, SO4, and SO5), maxillary area, the number of maxillary teeth, opercular area, area of the nares, and antorbital bone area (see Table S1).
Bony area measurements were compiled using the freehand tool in ImageJ (National Institutes of Health, Bethesda, MD). Raw values were obtained as the total number of pixels and converted to square millimeters based on the resolution of each raw image. For bones demonstrating fragmentation (SO3–5), we recorded total number of fragments, as well as individual area measurements for each fragment. In cases of fragmentation, we summed the individual fragment areas to recover the total area of the bone. Finally, we scored certain traits as binary, e.g., the presence (1) or absence (0) of fusion for bone pairs vulnerable to fusion, SO1+2, SO4+5, and SO5+6 (Table S1). We also performed additional analyses of bone area, accounting for size of the fish, in which bony area measurements were regressed against standard length of each specimen. QTL analyses were then performed based on the residual values for each bony area measurement (Table 3).
Linkage map construction
We recalculated linkage groups with the R/qtl package, using genotypic information for 175 markers (Protas et al. 2007) (Table S2). We identified 29 linkage groups for a total map length of 2545 cM and an average intermarker distance of 17.4 cM. The karyotypic number of A. mexicanus is 25 (Kavalco and Almeida-Toledo 2007). We attribute the higher number of linkage groups we identified to the reduced number of meiotic events analyzed, as well as potential genotyping errors. Recombination frequencies were collected for 539 individuals; however, craniofacial phenotypes were investigated in 237 cleared-and-stained individuals.
QTL and effect plot analyses
All QTL analyses were carried out using the software program R/qtl (Broman et al. 2003). Each trait was analyzed using marker regression (MR) (Kearsey and Hyne 1994), expectation maximization (EM) (Xu 2010), and Haley–Knott (HK) (Haley and Knott 1992) mapping methods. Most traits were evaluated as parametric or binary phenotypes (above). The number of meioses and genetic markers evaluated in our study was limited and may have reduced our power to detect QTL (Beavis 1998). Therefore, all phenotypes were analyzed using each of the three methods above. To identify prospective QTL from this analysis, we first used a low-stringency significance threshold (LOD = 3) to detect QTL that may otherwise be missed due to the number of genomic markers used and/or the number of hybrid individuals evaluated in our study.
We balanced this approach by discarding certain pseudomarkers, calculated using interval mapping, which harbored inflated LOD values and incongruous genomic positions compared to other mapping methods. These putative “ghost QTL” (Broman and Speed 1999) are artifacts that can arise using interval mapping procedures (Haley and Knott 1992; Martínez and Curnow 1992) and were discarded from further analysis. Effect plots were created for the closest linked markers, using phenotypic information for the left and right sides of each member of the pedigree. Genome-wide LOD significance thresholds (P < 0.05) were calculated in R/qtl for each phenotypic trait (Table 1), using a permutation test (of 1000 permutations) to identify statistically significant associations.
Following one-scan analyses, we performed multiple-QTL mapping (MQM), using R/qtl software (as described in Arends et al. 2010). MQM mapping is a powerful analytical tool for further understanding the genetic architecture of complex traits, such as craniofacial phenotypes (see Parnell et al. 2012). The automated MQM mapping procedure we utilized involves first “augmenting” missing genotypic data and then automatically selecting cofactors, using multiple-regression and iterative backward elimination procedures (Arends et al. 2010). Then, QTL were interval mapped using maximum likelihood (Table 2) and represented in circular genome interaction plots for each craniofacial phenotype (see examples in Figure 1, K and L) to identify prospective epistatic interactions between cofactors.
We also examined the potential effects of sex on our craniofacial phenotypes following procedures established by Broman and Sen (2009). We first performed an analysis of variance for 19 craniofacial trait values, using R/qtl. Accordingly, group means of males and females for each trait were compared (Table S3), followed by a QTL analysis with and without sex as a covariate. This technique utilizes standard interval mapping (Broman and Sen 2009), and therefore traits for which a QTL was detected using only marker regression could not be evaluated. The interval mapping method used in each analysis (EM or HK; Table S3) was selected based on the significant QTL reported for each trait (Table 1 and Table 3). All additional statistical analyses were carried out using Microsoft Excel (v.12.3.4 for Mac) or JMP (v.10).
Syntenic analyses and anchoring to the Danio rerio genome
We identified or estimated positions of linked microsatellite markers for the seven phenotypic traits that yielded significant QTL in this study. Marker position was identified in Danio rerio based on a previously generated “integrated” Astyanax linkage map (Gross et al. 2008). Eleven markers from the present study were not included in the prior analysis, so their homologous positions were estimated, as previously described (Gross et al. 2008). For each marker, we identified or estimated marker positions in the latest draft of the D. rerio genome (Zv9; release 72). Relative marker positions and synteny between the A. mexicanus linkage map and the D. rerio genome were developed as a circular representation, using the Circos software program (Krzywinski et al. 2009).
Seven novel craniofacial phenotypes demonstrate a genetic basis in Astyanax
We scored 33 craniofacial phenotypes on both the left and right sides of the head to determine the extent to which complex craniofacial alterations in a subterranean-dwelling fish are controlled by heritable factors. Seven of these traits yielded a genetic basis, four of which were present on both the left and right sides of the head, including (1) fusion of the first and second suborbital bones (SO1+2), (2) total bony area of the second suborbital bone (SO2), (3) total bony area of the third suborbital bone (SO3), (4) number of SO3 bony fragments, (5) total bony area of the fourth suborbital bone (SO4), (6) fusion of the fourth and fifth suborbital bones (SO4+5), and (7) total bony area of the fifth suborbital bone (SO5). We discovered a previously unappreciated asymmetry in the genetic architecture of three traits: SO1+2 fusion, SO3 fragmentation, and SO4+5 fusion, when comparing these traits on the right and left sides of the cranium (below).
Sex-specific body shape differences are apparent in A. mexicanus (Hinaux et al. 2011); however, dimorphic differences in cranial morphology have not been described. To evaluate the potential influence of sex on these craniofacial phenotypes, we performed a two-part analysis (see Materials and Methods) (Broman and Sen 2009). First, the results of an analysis of variance revealed only three traits differed significantly (P < 0.05) with respect to group means for phenotypic scores between males and females (Table S3). Interestingly, these traits included the residual value QTL identified for the SO2 bone (both right and left sides) and the SO2 area for the left side only (Table S3).
We then performed multiple covariate analyses, in which sex was the assigned covariate, in R/qtl for each trait (1000 permutations; see Broman and Sen 2009). We observed a difference in the LOD value for the QTL detected at genomic marker 227A for both right (R)SO2 area residuals and left (L)SO2 area residuals (Table S3), suggesting that this trait differs between males and females in our hybrid pedigree. However, this analysis did not reveal a difference in the strength of the LOD for the two markers (NYU27 and 229B) associated with LSO2 area. Thus, there appears to be a gender-specific difference in SO2 bone area, after correction for body length. For the majority of craniofacial traits analyzed, however, prior QTL associations did not change (P > 0.05), suggesting most of these traits do not differ between sexes (Table S3). We did detect a change in the LOD values for two other traits, LSO3 area and RSO3 number (Table S3); however, it is unclear whether these are biologically relevant since the ANOVA results for both traits were insignificant.
The majority of the craniofacial traits that we analyzed did not yield significant QTL. These included the total area of the opercular bone, the lateral cranial area of the head, the area of the nasal opening (nares), the area of the antorbital bone, and the area of the supraorbital bone. It was surprising that some of these traits did not yield a genetic effect. For instance, opercular bone shape varies considerably between different Astyanax cave populations (Jeffery 2009) and has undergone marked morphological changes within ray-finned fish (Kimmel et al. 2005). However, unlike that observed for bone area for four of the suborbital bones, we did not observe a genetic effect associated with opercular bone area.
In stickleback fish, this bone demonstrates significant variation between freshwater and oceanic forms (Kimmel et al. 2005, 2012a,b; Arif et al. 2009). Recent morphological analyses of the opercle bone revealed a major-effect locus associated with opercle bone shape in an F2 mapping pedigree of anadromous and lake stickleback fish (Kimmel et al. 2005). Perhaps future morphometric analyses will similarly reveal a genetic effect associated with shape differences in cavefish populations compared to surface morphs. Nonetheless, our collective results indicate that, while some regions of the skull are controlled by heritable genetic changes that evolved following colonization of the subterranean environment, other aspects of the skull remained unchanged.
Fragmentation and fusion phenotypes affecting the suborbital bone series demonstrate genetic asymmetry
Three phenotypes, including the number of SO3 bony fragments and fusion of the SO1+2 and SO4+5 bones, demonstrated a directional genetic basis in our hybrid pedigree. We discovered this pattern by observing that certain individuals with severe fragmentation phenotypes (e.g., five SO3 fragments) on one side of their cranium had less severe fragmentation on the contralateral side (Figure 1). Therefore, we scored both left and right sides of the cranium in our analyses.
The number of SO3 elements in our pedigree ranged from one to five on the right side and from one to three on the left side (Figure S1D). The most common phenotype was presence of a single bony element for both sides of the head. A QTL analysis identified two significantly associated markers, NYU53 (LODMR = 7.31) and 206A (LODMR = 5.43), for the right side SO3 number (Figure 1G). Effect plots revealed a significant effect of genotype for both loci on the right side only (Figure 1H). The heterozygous genotype was highly similar to the surface fish genotype, implying a dominant effect at this locus (Figure 1H). An insignificant phenotypic effect was noted for the same two loci on the left side (Figure 1J). Interestingly, bone size varied minimally between individuals, suggesting that the total area of the bone is constrained on right and left sides, despite uneven fragmentation patterns. We tested whether this result was determined by a few individuals harboring extreme fragmentation (e.g., four and five fragments) by reanalyzing SO3 fragmentation as a binary trait. This analysis identified the same loci, NYU53 and 206A, with lower but significant LODMR values (4.10 and 3.88, respectively). MQM also identified marker 206A (LODMQM = 4.06) as significantly associated with right SO3 fragmentation; however, no cofactors were observed (Figure 1K). For left-sided SO3 fragmentation, MQM analysis identified one marker and cofactor (233D) with an insignificant LOD score (2.15; Figure 1L).
An analysis of SO1+2 bony fusion similarly revealed a right-sided directional genetic effect (Figure 2). Fusion was scored as a binary trait and displayed a similar frequency on both sides of the head (Figure S1A). However, analysis of this trait identified two different loci that exceeded a LOD value of 3.0, 119C (LODMR = 3.28) and 229B (LODHK = 4.41), on the right side only (Figure 2, G and I). Effect plots revealed a dominant effect of the surface allele at marker 119C and an intermediate phenotypic effect at 229B (Figure 2H). MQM analyses supported the results of our one-scan mapping methods, identifying either the neighboring (NYU27; LG6) or the same markers (229B; LG27) for right-sided SO1+2 fusion (Figure 2K). Interestingly, two cofactors were identified on the same linkage groups (223C on LG27; 233D on LG6), which demonstrate a positive epistatic interaction for this trait (blue line, Figure 2K). Consistent with one-scan analyses, no significant markers or cofactors were detected using MQM for left-sided SO1+2 fusion (Figure 2L). As with the SO3 bone, the total area occupied by these bones did not differ between left and right sides.
SO4+5 fusion was also scored as a binary trait and revealed an association with marker 112A (LODMR = 3.03) only on the left side of the head (Figure 3). The effect plot for 112A indicates fusion is associated with the homozygous cave condition and that the heterozygous genotype produces an intermediate phenotypic effect (Figure 3, I and J). MQM analyses revealed the same genetic marker (112A; LODMQM = 2.79) as one-scan mapping methods and identified a cofactor (NYU31) that also resides on LG11 (Figure 3L). Consistent with one-scan results, no significant genetic markers were identified for right-sided SO4+5 fusion (Figure 3K). Despite the asymmetric association with fusion of these two elements, the overall area collectively occupied by the SO4 and SO5 bones did not differ across individuals.
Genetic symmetry is preserved in several cave-associated and craniofacial phenotypes
In a prior analysis using the pedigree evaluated here, eye size reduction was scored only on the right side of individual fish. We reanalyzed this phenotype and confirmed the identity of loci associated with quantitative changes in eye size for both left and right eye measurements (Table 1). We further validated the same QTL position for albinism, a Mendelian trait (Şadoğlu 1957a,b; Protas et al. 2006), using the previously published data set (Table 1).
Of the craniofacial traits analyzed, only area measurements for the SO2 and SO3 bones produced QTL that were identical when scored on the left and right sides (i.e., genetically symmetric). Analysis of SO4 and SO5 area revealed partially symmetric results. For the SO4 bone, one QTL was shared for both the right and left sides, while two additional markers were present on the left side only (Figure 6). For the SO5 bone, two nearby markers were detected on the right side (216C and NYU53) and one QTL was detected on the left side (2B); however, all three markers reside on the same linkage group (26; Figure 7). Area measurements for the SO2, SO3, and SO5 bones were distributed normally (Figure S1, B, C, and G); however, SO4 area measurements were positively skewed. We transformed SO4 area measurements to log10 values to create a normal distribution for QTL association studies (Figure S1E).
SO2 area is associated with three loci on both the left and right sides: 119C/NYU27 on linkage group 6, 209A on linkage group 7, and 229B on linkage group 27 (Figure 4). Despite alternative positions for 119C (right side; LODMR = 4.43) and NYU27 (left side; LODEM = 3.88), these loci are merely 3.2 cM apart and effect plots for both markers yield nearly identical patterns (Figure 4, G–J). The markers 209A (left side LODMR = 3.74; right side LODMR = 3.42) and 229B (left side LODEM = 4.28; right side LODHK = 3.94) were associated with increased bony area in cavefish (Figure 4, H and J). MQM analyses largely supported the results of one-scan methods, identifying three markers on the right side (NYU27, 209A, and 229B; Figure 4K), but only two markers on the left side (209A and 229B; Figure 4L). MQM analyses revealed four cofactors (Figure 4K) for right-sided SO2 area, but only two for the left side (Figure 4L). Interestingly, a positive epistatic interaction was observed between two cofactors (NYU25 and 110B) on the right side, and a negative epistatic interaction was found between two markers (223C and 110B) on the left side (Figure 4, K and L). This result demonstrates a complex interaction between the same marker (110B) and two different cofactors (NYU25 and 223C) on different linkage groups, which exert contrasting epistatic effects on the right and left sides, respectively.
SO3 area is associated with two loci on both right and left sides: 55B on linkage group 18 (left side LODHK = 4.02; right side LODEM+HK = 3.24) and 229B (left side LODEM = 3.28; right side LODEM+HK = 3.66) on linkage group 27 (Figure 5). Interestingly, both SO2 and SO3 areas shared a genetic association with the same genetic marker (229B). Effect plots for both bone area phenotypes reveal the same polarity, wherein two copies of the cave allele are always associated with larger SO2 (Figure 4, H and J) and SO3 bones (Figure 5, H and J). The heterozygous genotype for marker 55B produces an intermediate phenotype between the cave- and surface-dwelling phenotypes (Figure 5, H and J); however, the heterozygous genotype for marker 229B is similar to the cave-dwelling phenotype, possibly indicating dominance of the cave allele at the linked locus (Figure 5, H and J). MQM analyses identified the same QTL as one-scan methods for the both the right (55B, LODMQM = 3.38; 229B, LODMQM = 2.9) and left (55B, LODMQM = 4.16; 229B, LODMQM = 3.02) SO3 bones (Figure 5, K and L). Moreover, the same cofactors were selected for both sides (131C, 223C) with the exception of marker 26A, which was identified on the right side only. No epistatic interactions between cofactors were observed on either side (Figure 5, K and L).
SO4 area demonstrated partial genetic symmetry (Figure 6). The same genetic marker, 214F, is associated with area changes on both the left and right sides (Figure 6, G and I), although the LOD value is higher on the left side. Two additional genetic markers, 222E and 229B, are associated with SO4 area differences on the left side only (Figure 6I). Interestingly, the effect plot for the symmetric marker 214F demonstrates intermediate dominance on the right side (Figure 6H), but dominance for the cave allele on the left side (Figure 6J). Marker 222E demonstrates a weak dominant effect for the surface allele, while marker 229B demonstrates intermediate dominance (Figure 6J). Interestingly, the phenotypic polarity is different for the three markers: cave alleles for 222E and 214F are associated with smaller SO4 area, while cave alleles for 229B are associated with a larger SO4 area (Figure 6, H and J). MQM analyses supported the notion of partial genetic symmetry for SO4 area (Figure 6, K and L), identifying the same genetic marker (229B) on both right (LODMQM = 3.12) and left sides (LODMQM = 3.81). However, only one cofactor was selected for the right side (223C; Figure 6K) compared to three cofactors that were selected for the left side (202E, NYU53, and 223C; Figure 6L).
SO5 area is associated with the same linkage group (26) on both the left and right sides of the cranium (Figure 7). The effect plots for identified markers on the right (216C, NYU53) and left sides (2B) are nearly identical, demonstrating dominance at the cave locus (Figure 7, I and J). In both cases, the phenotypic polarity is the same—cave alleles are associated with a smaller SO5 area. One marker, 222E, harbored a LOD value just below our 3.0 threshold on the left side (Figure 7H). The same marker on the right side demonstrated an insignificant LOD value (Figure 7G). MQM analyses identified different results for this trait compared to one-scan analyses. On the right side, marker 214F on LG26 demonstrated an association with SO5 area (LODMQM = 2.29; Figure 7K). However, on the left side, a different marker on LG27 (229B) demonstrated an association with SO5 area (LODMQM = 2.88; Figure 7L). Two cofactors were identified on the right side (NYU53 and 233D); however, no cofactors were selected on the left side.
We performed additional analyses of area for the SO2, SO3, and SO5 bones based on residual values following normalization for standard length of each specimen. Interestingly, these results identified both symmetric and asymmetric QTL. For instance, a symmetric genetic signal was detected for right- and left-sided SO2 area at marker 227A (LODMR = 5.56 and 5.69 for the right and left sides, respectively). However, one marginally significant QTL (106C; LODMR = 4.0) was detected only on the right side (Table 3). SO2 area effect plots again revealed a larger bone associated with the homozygous cave genotype (“CC”).
Residual values for the SO3 bony area yielded similarly mixed results. Both 26A and 136B were identified for both the left and right sides; however, two additional markers (203F and 6A) were detected on the left side only. QTL effect plots from this analysis yielded contrasting results. Markers 136B, 203F, and 6A were all associated with larger values in the cave homozygous condition; however, marker 26A was associated with a smaller area. These results may explain, in part, the discrepancy between our results and those of Yamamoto et al. (2003), particularly if the causative gene(s) associated with marker 26A exerts a stronger effect at the earlier developmental stages evaluated in their study. At present, however, the gene(s) mediating this potential effect remains unknown. Finally, residual values for the SO5 bone also demonstrated asymmetry since significant QTL were strictly detected on the left side (marker 2B; LODMR = 3.61). Interestingly, this marker demonstrated a smaller value associated with the homozygous cave condition.
QTL associated with craniofacial phenotypes map to syntenic regions of the D. rerio genome
At present, there is not a fully annotated genome sequence available for A. mexicanus. Therefore, as a step toward nominating candidate genes linked to the QTL we detected, we analyzed the positions of microsatellite loci on an integrated linkage map that is anchored to the D. rerio genome (Figure 8). This integrated map represents a comprehensive analysis of multiple cave × surface fish crosses and therefore is a proxy for the genomic positions of several markers (Gross et al. 2008). Three analyses have identified numerous syntenic blocks shared between Astyanax and D. rerio (Gross et al. 2008, 2013; O’Quin et al. 2013). These blocks were identified through sequence similarity (BLAST hit) queries of noncoding sequences flanking polymorphic microsatellites or SNP markers used in RAD-seq studies.
The microsatellite marker 119C on linkage group 1 (Figure 8A), linked to right-sided SO1+2 fusion and SO2 area, lies at approximate position 24.2 Mb on chromosome 24, ∼8.8 cM away from the gene musculin (msc). The microsatellite marker NYU27, which is linked to left-sided SO2 area, is positioned on chromosome 24 in D. rerio (Figure 8A). Marker NYU27 is ∼3 cM away from 119C and is estimated to reside near the same syntenic block on chromosome 24 in Danio. Marker 229B is linked to a number of the traits we analyzed, including right-sided SO1+2 fusion, symmetric SO2 bony area, symmetric SO3 bony area, and symmetric SO4 bony area phenotypes (Figure 8A). This marker is predicted to be ∼0.83 cM away from the gene alpha-A-crystallin (crystaa) within a 38.4-Mb syntenic block of Danio chromosome 1.
The microsatellite marker 209A (Figure 8A), linked to symmetric SO2 bony area, resides in a syntenic block on chromosome 20 in D. rerio, at approximate position 55.9 Mb. Microsatellite marker 55B, which is linked to the symmetric SO3 bony area phenotype, comprises a short (∼0.1 Mb) syntenic region on Danio chromosome 2 along with the marker 135C (Figure 8A). Marker 206A, on integrated linkage group (LGI) 23, which produces the peak LOD score for right-sided SO3 fragmentation, is predicted to reside within a 27.2-Mb syntenic block of Danio chromosome 17 (Figure 8, A and B). This marker is located ∼1.78 cM away from bone morphogenetic protein 4 (bmp4). Marker NYU53, also linked to right-sided SO3 fragmentation and right-sided SO5 area, was not evaluated in our prior publication. However, this marker is ∼2 cM away from marker 216C in our current linkage map, which is anchored to a syntenic block of Danio chromosome 22 (Figure 8A).
Microsatellite marker 214F, linked to symmetric SO4 bony area, resides near a block of synteny with Danio chromosome 22 (Figure 8A). The marker 112A, linked to left-sided SO4+5, fusion maps between two genetic markers (224C, ccng1), both of which are anchored to Danio chromosome 14 (Figure 8A). Marker 2B, linked to symmetric SO5 bony area, maps near a syntenic region shared between Astyanax LGI 12 and chromosome 22 (Figure 8A).
Genetic asymmetry in the cavefish craniofacial complex
Here we demonstrated genetic asymmetry in SO3 fragment number and fusion of the SO1+2 and SO4+5 bones. Prior work has established a role for fluctuating asymmetry in the generation of numerous aberrant morphologies, including wing vein patterns in Drosophila (Klingenberg and Mcintyre 1998), male sexual ornaments in swallows (Møller 1990), and mandibular bone landmarks in mice (Leamy 1993). This phenomenon is frequently attributed to developmental instability (Klingenberg and McIntyre 1998) caused by any of a number of factors such as environmental stress (Palmer 1994), interspecific biotic stresses such as pathogens (Watson and Thornhill 1994), intraspecific biotic stress such as competition (Manning and Chamberlain 1993), or genetic stress occurring as a result of hybridization (Clarke et al. 1986; Hochwender and Fritz 1999; Bourguet 2000; Dongen 2006), which leads to random-sided morphological aberrations. Some of the craniofacial changes we report here may be fluctuating asymmetrically. In our mapping cross, all of the craniofacial traits for which we observed a genetic effect appear with roughly equal frequency on both sides of the head (Figure S1). Interestingly however, the genetic effects for three of these craniofacial traits were detected only on one side. We observed a genetic effect on strictly the right side for two traits (SO3 fragment number and SO1+2 fusion) and a genetic effect on strictly the left side for one trait (SO4+5 fusion).
A few prior studies have reported the genetic basis of craniofacial asymmetry. Leamy et al. (1997) evaluated directional asymmetry in a cross of Large (LG/J) and Small (SM/J) inbred mouse strains and discovered 16 significant QTL associated with nine mandibular characters (Leamy et al. 1997). Portuguese water dogs commonly develop hip dysplasia, which arises from both heritable and environmental factors (Chase et al. 2004). Chase et al. (2004) evaluated the genetic basis for this trait through quantitative autoradiographic measurements of hip joint laxity of the hind legs. They discovered two QTL, one affecting the right hip (explaining 16% of the heritable variation) and one affecting the left hip (explaining 14% of the heritable variation). There was significantly greater laxity for the left vs. the right hip, and this polarity was conserved in 80% of the dogs they evaluated. This directional asymmetry could be linked to behavioral or physiological preferences (Chase et al. 2004). For instance, a genetically influenced behavior (“right-” or “left-footedness”) may lead to undue joint strain and result in a more severe phenotype on one (but not the other) side of the animal.
Across many generations, this form of heritable variation may lead to more significant morphological asymmetries [e.g., claw size in fiddler crabs (Morgan 1923; Leamy et al. 1997)]. A fascinating example of morphological asymmetry has been described in cichlids. The Lake Tanganyikan scale-eating fish, Perissodus microlepis, harbors two morphotypes that differ in the (left–right) direction of their mouth opening (Hori 1993). These fish attack and feed on the scales of other living fish based on their handedness, which was originally reported as a simple Mendelian trait (Hori 1993). Recently, Stewart and Albertson (2010) identified a jaw laterality locus segregating in an F2 cross of two herbivorous species from Lake Malawi based on length variation of the retroarticular processes of the jaw. The patterns of inheritance for this were consistent with the genetic effects predicted by earlier studies (Hori 1993; Hori et al. 2007). The ratio of right- and left-handed morphs within the population is determined by a mechanism of frequency-dependent natural selection, mediated by the alertness of the prey (Hori 1993). Thus, natural selection appears to have favored the evolution of these distinct jaw asymmetries to facilitate predation success (Stewart and Albertson 2010). Within the wider context of the Perissodini tribe, these asymmetries likely arise in response to shifts in habitat and predation strategy (Stewart and Albertson 2010).
It is difficult to imagine how such “handedness” could be adaptive in blind Mexican cavefish. However, cranial neuromasts (which mediate responses to water movements) densely populate the SO3 bone (Yoshizawa et al. 2012), which demonstrated the most severe fragmentation phenotype in our hybrid pedigree (Figure 1). Perhaps increased cranial neuromast density leads to bone fragmentation by interfering with normal condensation and ossification of the underlying cranial mesenchyme (Stensiö 1947; Lekander 1949; Reno 1969). If so, a consistent directional asymmetry affecting cranial neuromast number might also be predicted. At present, a connection between increased neuromast numbers and severity of fragmentation, or an attendant directional asymmetry in cranial neuromast number in natural populations, remains unknown.
Novel craniofacial QTL in the context of past studies
Several of the markers identified from our analysis overlapped with QTL positions identified by Protas et al. (2007) (Table 1). For instance, four microsatellites—234B, NYU14, NYU27, and 229B—colocalized with map locations for previously identified eye and lens QTL. In addition, two genomic markers (119C, 206A) mapped to positions associated with melanophore number, one marker (55B) mapped to a linkage group position associated with taste bud number, one marker (112A) is also associated with the maxilla, and two markers (229B, 2B) mapped to a region associated with tooth phenotypes (Protas et al. 2007). We did not observe colocalization of our craniofacial traits with the SO3 width trait discovered by Protas et al. (2008). However, one of the six SO3 width QTL reported by Protas et al. (2008) maps to marker 203E, which neighbors marker 55B. SO3 area phenotypes also map with a peak LOD score to this marker (Protas et al. 2008). Thus, this marker may broadly indicate a genetic basis for bone size differences between morphotypes, and our scoring methods may have redundantly reported the same underlying phenotypic difference.
Of the 11 craniofacial phenotypes that yielded QTL in our study, only four linked markers failed to map to regions associated with previously identified troglomorphic phenotypes. This may indicate that certain regions of the cavefish genome, populated with multiple different genes, have undergone accelerated evolution through the process of cave adaptation (Yokoyama and Yokoyama 1990). Alternatively, the craniofacial QTL identified, along with previously discovered QTL, may be mediated by the same pleiotropic genes (Protas et al. 2008) or be closely linked to the causative gene(s) controlling other cave-associated traits. Future studies evaluating a larger pedigree along with a higher-resolution physical map of the genome will enable us to discern between these two possibilities.
Yamamoto et al. (2003) demonstrated that eye removal results in both dependent and independent effects on the craniofacial complex. One trait concluded to be independent of eye removal was the number of SO3 bony fragments (Yamamoto et al. 2003). Our results support this conclusion; the two SO3 fragment number QTL we identified (linked markers 206A and NYU53) do not overlap with the positions of any eye or lens QTL (Protas et al. 2008). However, 206A maps to a previously discovered QTL for melanophore number (Protas et al. 2007). Moreover, the genetic effect we discovered was detected only on the right side, implying a genetic asymmetry in SO3 bone fragmentation. Yamamoto et al. (2003) hypothesized that SO3 area may be reduced as an indirect consequence of olfactory sense expansion in cavefish. In this study, we measured the area of the external opening to the olfactory epithelium to determine whether nares expansion is associated with SO3 fragmentation, but found no QTL associated with this trait.
Interestingly, effect plots for SO3 bone QTL revealed homozygous cave alleles are associated with a larger SO3 bone area (Figure 5, H and J). In contrast, Yamamoto et al. (2003) observed that the size of the circumorbital bones correlated with presence of the eye in their experimental studies; i.e., eyeless experimental subjects had smaller circumorbital bones (Yamamoto et al. 2003). Our differing results could be attributed to at least two explanations. First, in experimental manipulations, bone growth may be critically reliant upon the developing eye. However, under nonexperimental conditions, the genetic basis for SO3 bone size increase may be fundamentally distinct from that for eye loss in cavefish. Second, different-aged specimens were evaluated in Yamamoto et al. (2003) compared to the present study. We therefore wanted to determine whether larger-sized bones were observed, in part, due to allometric size increases in this bone that may be revealed by the older specimens used in our study.
A regression analysis of the individuals within our pedigree indicated positive relationships between bone area measurements (for both the right and left sides of the cranium) and standard length (Table 3). The SO2, SO3, and SO5 bones are all positively associated with standard length in our hybrid pedigree, with R2 values ranging from 0.37 to 0.75, suggesting that larger-sized specimens harbor larger circumorbital bones. We further investigated this relationship by scanning for QTL based on residual values obtained for each of our circumorbital bone measurements. This was carried out following correction for standard length of each specimen (Materials and Methods). Significant results were obtained only for the SO2, the SO3, and the left SO5 bony area measures (Table 3).
These residuals analyses allowed us to dissect apart the genetic effects associated with suborbital bone areas that are not indirectly influenced by the overall size of each specimen. For example, larger fish harbor larger bones, and therefore QTL identified using raw values alone may indirectly report associations based on overall size differences. Indeed, for the SO2 and SO3 bones (Table 3), our residuals analyses identified new QTL (i.e., on different linkage groups) compared to those detected using raw, uncorrected values. Interestingly, however, these genetic effects were largely symmetric, consistent with our analysis of raw area measurements. Thus, normalization based on body size provided us with additional sensitivity to detect relevant genetic loci involved in complex craniofacial features evolving in cave-dwelling fish. These results further underscore the complex genetic regulation of craniofacial diversity in this species.
Potential genes affecting bone fragmentation and fusion
Leamy et al. (1997) observed QTL associated with directional asymmetry in mandibular dimensions evaluated in a cross of two inbred mouse strains. The authors speculated that directional asymmetry might be governed by genes encoding hormones, growth factors, and/or hormone receptors (Leamy et al. 1997, 1998). They noted the positions of three candidate genes, including growth arrest-specific 1 (gas1) on mouse chromosome 13 and thyrotropin-releasing hormone receptor (trhr) and platelet derived growth factor, endothelial cell (ptgfec) on mouse chromosome 15 (Leamy et al. 1997). A subsequent study identified QTL associated with fluctuating asymmetry, many of which localized to similar chromosomal regions to those identified for directional asymmetry (Leamy et al. 1998). Our preliminary data do not indicate that any of these genes reside close to the QTL identified in our study, with the exception of the genes ptgfec and htr4, which reside near marker 112A, and bmp4, which maps close to marker 206A. Overall, this observation may indicate that the causative genetic lesions associated with directional asymmetry are divergent in our model system or may be governed by one or more paralogous genes in Astyanax, whose homologous genomic positions differ significantly from the genomic position of these genes in mice.
One prospective candidate gene, bmp4, plays a key role in the evolution and craniofacial development of a broad set of species (Schubert et al. 2000; Terai et al. 2002; Trainor et al. 2003). For instance, bmp4 is associated with adaptive jaw shape differences in cichlid fish (Albertson and Kocher 2006) as well as the depth and breadth of upper beaks in Darwin’s finches (Abzhanov et al. 2004). Alongside its role in craniofacial evolution and development, bmp4 has also been implicated in laterality of the heart and gut (Monteiro et al. 2008). This wider role in establishing left–right asymmetry in other organ systems during development renders bmp4 a particularly attractive candidate molecule for mediating asymmetric phenotypes in our system.
While the genes controlling suborbital bone fragmentation in cavefish remain unknown, bmp4 may play a critical role in this abnormal process. The SO3 dermal bone in cavefish develops from condensations of mesenchymal cells derived from the embryonic cranial neural crest (Northcutt and Gans 1983; Gross and Hanken 2008). Therefore, the origin of directional asymmetry in SO3 bone fragmentation may ultimately reflect aberrant migration, proliferation, or survival of the osteogenic cranial neural crest (in a directionally biased manner) during embryonic development.
Another prospective candidate gene mediating craniofacial asymmetry in our system is fgf8, which influences laterality in the zebrafish and mouse craniofacial skeletons (Albertson and Yelick 2005, 2007; Griffin et al. 2012). Specifically, zebrafish juveniles with attenuated fgf8 activity demonstrate consistent left–right asymmetry in craniofacial structures (Albertson and Yelick 2005). Using the transgenic aceti282a/fgf8 heterozygous zebrafish line, Albertson and Yelick (2007) further demonstrated a role for fgf8 in mediating a variety of abnormalities, including facial asymmetry, altered craniofacial geometry and sutural pattern, and ectopic bone formation. These abnormalities, which resemble certain human disorders (such as craniosynostosis and hemifacial microsomia), may originate from increased bone metabolism observed in the aceti282a/fgf8 transgenic line (Albertson and Yelick 2007). A study of fgf8 attenuation in mice, combined with gain-of-function mouse–chicken chimeric studies, suggests a deeply conserved role for fgf8 in craniofacial development and symmetry. Griffin et al. (2012) recently demonstrated a key role for fgf8 in midfacial integration and development of the optic capsule. Thus, fgf8 is an excellent candidate gene for future analyses in cavefish given that fgf8 signaling is a likely target for evolutionary selective pressures (Griffin et al. 2012). The gene fgf8b is predicted to reside near marker 229B, associated with SO1+2 fusion in the present study (Figure 8).
Alternative explanations for genetic asymmetry
Serotonergic signaling has been implicated in consistent left–right asymmetry in vertebrate embryos (Fukumoto et al. 2005; Levin et al. 2006). Serotonin is a neurotransmitter molecule with critical functions in a variety of processes, including physiology (Fuller 1992), cognition (Canli and Lesch 2007), and circadian rhythms (Moore and Klein 1974). In Xenopus embryos, proper serotonergic signaling through the 5-HTR3, 5-HTR4 receptors and the metabolism of monoamine oxidase (MAO) is essential for establishment of asymmetric gene expression and organ development (Fukumoto et al. 2005). Serotonin is functional very early in development through maternal deposition, demonstrates dynamic spatial localization, and exerts its action as early as the four-cell stage (Fukumoto et al. 2005). Fukumoto et al. (2005) did not evaluate craniofacial bone formation, which occurs later in development. However, prior studies in mice have indicated serotonergic signaling through the 5-HT2B receptor is critical for craniofacial morphogenesis (Choi et al. 1997). Thus, altered expression patterns of serotonin, its receptors, or related enzymatic pathways (e.g., through MAO) may have similarly evolved in the cavefish lineage. Our syntenic analysis does not indicate that these particular receptors map near the critical region of our QTL (Figure 8); however, we cannot rule out the role of a paralogous gene or other member of the serotonergic pathway in causing one or more directionally symmetric phenotypes in cavefish.
Another possibility is the asymmetric activation of a microRNA locus (He and Hannon 2004). In Caenorhabditis elegans, a functional asymmetry is generated during development through the two-step activation of lsy-6 in a postmitotic pair of sensory neurons (Chang et al. 2004). lsy-6 is first “primed” in the precursor of the left (but not right-sided) neuronal cell and then boosted to higher expression levels several divisions later (Chang et al. 2004). The downstream consequence of functionally relevant levels of expression leads to asymmetric functional identity, based on the expression of distinct guanylate cyclase receptors (Chang et al. 2004). In cave-dwelling vertebrates, such altered patterns of expression for molecules leading to directionally asymmetric craniofacial phenotypes may similarly be at play. These changes might have arisen randomly through drift or as an indirect consequence of selection for another unrelated function.
The evolutionary mechanism of regressive craniofacial trait evolution
The craniofacial traits we describe here may not evolve as a (direct or indirect) consequence of selection pressure. In the complete darkness of the cave, there are relaxed selection pressures associated with visually mediated traits (Borowsky and Wilkens 2002). In the above-ground environment, craniofacial symmetry may be under intense sexual selection, as an indicator of mate quality (Rhodes et al. 2001). Conversely, in the subterranean environment perhaps total darkness and loss of vision lead to a relaxation of this selective constraint, resulting in the asymmetries we observed in this hybrid cross.
Here we provide further evidence that organismal development can be influenced differentially on the left and right sides, leading to asymmetric phenotypes. Such asymmetries, being controlled by particular loci, may provide a substrate upon which natural selection can act to elaborate phenotypic differences along the left–right axis of the organism. We also demonstrate here that certain craniofacial QTL do not colocalize with previously identified cave-associated traits, while other QTL are associated with a variety of traits such as vision and pigmentation (Table 1). It remains unclear whether the latter result indicates pleiotropy or the close physical linkage of the causative genes for the phenotypes evaluated in this study. Future studies directed toward identifying the precise genetic lesion(s) explaining these traits will enable functional analyses. Functional validation will, in turn, clarify whether multiple phenotypes are affected by gain or loss of encoded gene products. Ultimately, the evolutionary mechanism (selection vs. drift) driving craniofacial changes in cavefish may differ, depending on the trait.
We thank the members of the Gross laboratory for helpful discussions. We especially thank Meredith Protas, Richard Borowsky, and Cliff Tabin for providing the cleared-and-stained specimens evaluated in this study. K. Dougherty and L. Bruns assisted with phenotypic measurements during an earlier phase of this project. We also thank Krista Nichols and two anonymous reviewers for providing helpful suggestions on an earlier draft of this manuscript. This project was supported by National Institutes of Health (National Institute of Dental and Craniofacial Research) grant DE022403 (to J.B.G.).
Communicating editor: K. M. Nichols
- Received October 24, 2013.
- Accepted January 28, 2014.
- Copyright © 2014 by the Genetics Society of America