Contrasting patterns of RUNX2 repeat variations are associated with palate shape in phyllostomid bats and New World primates

Establishing the genetic basis that underlies craniofacial variability in natural populations is one of the main topics of evolutionary and developmental studies. One of the genes associated with mammal craniofacial variability is RUNX2, and in the present study we investigated the association between craniofacial length and width and RUNX2 across New World bats (Phyllostomidae) and primates (Catarrhini and Platyrrhini). Our results showed contrasting patterns of association between the glutamate/alanine ratios (Q/A ratio) and palate shape in these highly diverse groups. In phyllostomid bats, we found an association between shorter/broader faces and increase of the Q/A ratio. In New World monkeys (NWM) there was a positive correlation of increasing Q/A ratios to more elongated faces. Our findings reinforced the role of the Q/A ratio as a flexible genetic mechanism that would rapidly change the time of skull ossification throughout development. However, we propose a scenario in which the influence of this genetic adjustment system is indirect. The Q/A ratio would not lead to a specific phenotype, but throughout the history of a lineage, would act along with evolutionary constraints, as well as other genes, as a facilitator for adaptive morphological changes.

region appears to play an important role in the repression and activation of essential proteins for craniofacial development 18 . Clinical genetic studies have demonstrated that haplo-insufficiency of the RUNX2 gene is associated with cleidocranial dysplasia (CCD) in humans, an autosomal dominant disease characterized by systemic skeletal anomalies, causing shortening of the face, abnormalities in the permanent dentition, as well as stature effects in affected individuals [19][20][21][22][23][24] . Fondon and Gardner 25 evaluated the correlation between different amino acid repeat motifs present in developmental genes and their effects across dog breeds with diverse skeletal and craniofacial phenotypes. They found a positive correlation between the ratio of polyglutamines (Q) to polyalanines (A) in the Q/A repeat domain of the RUNX2 gene and craniofacial length. Similar patterns of significant positive correlation between the RUNX2 repeat and facial length was also observed in other mammals belonging to the order Carnivora 10 . In the same study, based on functional tests, Sears et al. 10 found significant correlations between the ratios of polyglutamines to polyalanines in the RUNX2 tandem repeat region and RUNX2 transcriptional activity. Other studies also provided evidence for the polyQ and polyA regions as being a mechanism of fine tuning in the activation and repression of the RUNX2 protein 10,16 . Conversely, Pointer et al. 11 found no correlation between RUNX2 sequence and face length across a variety of non-carnivoran placental mammals. Also, recently, Newton et al. 12 demonstrated that RUNX2 does not drive craniofacial diversity in marsupials. Taken together, these previous studies suggest that Q/A repeat domain of RUNX2 might play a role in craniofacial evolution in some mammalian lineages but not others. Therefore, the role of the Q/A repeat domain of RUNX2 in craniofacial evolution lineage must be investigated independently for each lineage.
In the present study we combine phylogenetic comparative methods and morphological integration analysis to investigate the association between craniofacial length and width with RUNX2 Q/A ratios across two mammalian lineages displaying large variation in skeletal terms, the New World leaf-nosed bats (Family Phyllostomidae) and Primates (Catarrhini and Platyrhini). The investigation of these taxa is of interest to understand the developmental mechanism involved in the expression of RUNX2.
The New World leaf-nosed bats, Family Phyllostomidae, is a Neotropical lineage that evolved during the last 30 million years 26 and represents one of the most ecologically diverse families of mammals 27 . They present a marked diversity in feeding specialisations, including species that consume insects, blood, small vertebrates, fruits, nectar and pollen. Feeding specialisation is accompanied by morphological adaptations, such as conspicuous craniofacial variability (Fig. 1A) [27][28][29] . Few studies have explored the genetic basis underlying diversity in craniofacial phenotype in bats. One of the rare studies to consider this subject was conducted by Phillips et al. 30 , who examined the evolutionary dynamics of the Pax9 transcription factor and craniofacial diversification in bats. Subsequently, Ball et al. 16 found an association between the RUNX2 variation and cranial bone density in bats and other mammals. On the other hand, the craniofacial morphological variation across phyllostomid lineages has been widely documented through biomechanical, functional, and developmental studies 27,29,31-42 . From a comparative perspective, primates show remarkable morphological variation (Fig. 1B). The Catarrhini evolved during the last 35 million years 43 and represents a monophyletic lineage encompassing all Old World monkeys (OWM) and apes, including humans. Catarrhines exhibit a diverse range of life histories, habitat use, and dietary strategies, and their geographic distribution includes tropical ecosystems in Africa and Asia 44 . The evolutionary diversification of the Platyrrhini (NWM) occurred over a period of more than 30 million years, driven by diet and size-related selection 45,46 . Species from both Catarrhini and Platyrrhini display diverse cranial morphologies and remarkable variation in body size. Among the studies with hominids, Green et al. 47 hypothesized that an evolutionary change in RUNX2 affected the upper body and cranial morphology of modern humans. Additionally, as demonstrated by Lindskog et al. 48 , humans and macaques showed differences in protein and RNA expression levels in the RUNX2 functional network genes (e.g. RUNX1, SMAD3, SOX6, and SATB). Other genomic studies showed a correlation between RUNX2 and nose bridge breadth in current human populations 49 . However, in several other primate lineages, such as the NWM, little is known about the genetic basis that controls the fine-tuning of craniofacial phenotypic variation.

Results
Morphological variation and RUNX2 repeat variability. Craniofacial variation is a conspicuous trait of phyllostomid bats. We observed marked differences in palate shape across species representing all feeding strategies (electronic supplementary material, Table S1 and Fig. 1A Similarly, great variation in RUNX2 Q/A repeats and flanking control regions was evident both in phyllostomid bats and primates (Table 1). Overall, phyllostomid bats presented greater variation in number of glutamate and alanine residues than did primate species (Fig. 1C). Frugivorous phyllostomids presented the widest range, with a Q/A ratio varying from 1.21 to 2.33 (Table 1 and Fig. 1C).
Correlation between Q/A ratios and palate length and width. Spearman's correlation test for the phyllostomid bats data showed a significant negative correlation between Q/A ratio and palate length r = −0.51 (p = 0.01) and significant positive correlation between palate width and Q/A ratio, r = 0.55 (p < 0.01) ( Fig. 2A,B). Spearman's correlation test for the primates data (NWM and OWM analysed together), and for OWM (analysed separately) showed no significant association between the Q/A ratio and craniofacial morphological traits (NWM + OWM: palate length × Q/A ratio: r = 0.44, p = 0.12; palate width × Q/A ratio: r = 0.07, p = 0.80; OWM: palate length × Q/A ratio: r = 0.58, p = 0.17; palate width × Q/A ratio: r = 0.29, p = 0.52). For NMW (analysed separately), we found a positive marginally significant association between palate length and Q/A ratio (r = 0.79, p = 0.05) (Fig. 2C).
Bayesian PCM was also used to estimate the correlation coefficient between palate and Q/A ratios in both groups, corrected for phylogenetic relatedness, evolutionary rates and divergence times. In accordance with the previous findings, these analyses showed a positive correlation between Q/A ratio and palate width (r = 0.2910; pp = 0.94), and a negative correlation between Q/A ratio and palate length (r = −0.2260; pp = 0.13) for phyllostomid bats ( Table 2). The craniofacial variation of primates in general did not correlate with Q/A ratios. However, when the dataset is split in two subgroups (NWM and OWM), we found a significant correlation between palate length and Q/A ratios for NWM (r = 0.583, pp = 0.93) For phyllostomids, the PGLS results showed a phylogenetic signal (Pagel's lambda = 1) associated with palate length and width, and with zygomatic width ( Table 3). The regression between phenotypic traits and Q/A ratios showed that palate width is positively associated with RUNX2 variation (adjusted R² = 0.19, p = 0.01). For primates (NWM and OWM analysed together) the PGLS results showed a phylogenetic signal associated with palate length and width. In that case, we found no significant association between phenotypic traits and Q/A ratios (Table 3). For NMW (analysed separately), the PGLS results presented phylogenetic signal associated with palate length and width, and a positive marginally significant association between palate length and Q/A ratio (adjusted R² = 0.55, p = 0.06) ( Table 3). Finally, in the OWM, we found phylogenetic signal associated with palate length and no association between phenotypic traits and Q/A ratios in the PGLS analyses (Table 3).

Correlation between Q/A ratios and morphological integration index.
The morphological integration index (ICV), obtained for phyllostomid bats, exhibited strong negative correlation with Q/A ratios (r = −0.6, p = 0.02; Fig. 3A). Additionally, the evolutionary flexibility index was strongly correlated with Q/A ratios (r = 0.7, p < 0.01; Fig. 3B), and the constraints index exhibited negative correlation with Q/A ratios (r = −0.6, p = 0.03; Fig. 3C). Spearman's correlation test for the primates data analysed together and separately showed no significant association between the Q/A ratio and ICV, flexibility and constraints indexes. The lack of correlation observed for primates data might be due to the relatively small sample sizes for both NWM and OWM clades.

Discussion
Frugivorous phyllostomid bats: a hotspot of variation. The Phyllostomidae family includes subfamilies that are recognised mainly by their feeding habits and associated morphological adaptations 27,28,51 . Phytophagous phyllostomids (frugivorous and nectar feeding species) evolved in nested clades among ancestral animalivorous species, and previous studies have proposed that the insectivorous ancestor of all phyllostomids also fed on plants 28,31,51,52 . In the present study, the largest variation in the Q/A ratios was found in frugivorous phyllostomids (Q/A ratio varying from 1.21 to 2.33; see Table 1). The feeding specialisation in the subfamily Stenodermatinae varies from primary frugivory (feed on fruits and complement diet with other items, e.g. Artibeus, Carollia, Sturnira,

Magnitude of skull integration in bats and primates: flexibility to evolve.
In contrast to the remarkable stability observed in the patterns of covariation in the mammalian skull 54-57 the magnitudes of inter-trait association, estimated by the average squared phenotypic correlation between traits (r 2 ), are highly variable 56 . For the same cranial measurements, primates (chimpanzees, humans, and gorillas) and Chiroptera   56,58 . The high variability in the magnitude of inter-trait correlations has potential important consequences for the evolution of these clades 59 . Higher magnitudes of skull integration are associated with constrained evolutionary responses to selection, which would reduce adaptive flexibility to evolve in directions other than those driven by the tight inter-trait associations, at least on short-term 46,59,60 . Alternatively, lower magnitudes indicate that most traits are weakly connected and that there is more potential to evolve across several axes of variation. Species with lower magnitudes of trait association could evolve in more directions, therefore presenting less potential evolutionary constraint and greater evolutionary flexibility 59 . Thus, the different integration magnitudes in phyllostomid bats in relation to the RUNX2 gene ( Fig. 3A-C) together with the probable differences in the reaction norms of this gene within this clade, could explain the differences in the genotype-phenotype relationships observed, as well as the contrasting patterns described for other lineages of mammals.
The mice genome informatics database (http://www.informatics.jax.org/vocab/mp_ontology/MP:0020039) lists several known skeleton phenotypes compiled so far and associated with RUNX2 mutants. In this context, we could assume that delayed or increased bone ossification (both endochondral and intramembranous) might explain the opposed association of RUNX2 variation and/or integration and facial forms we observe in different groups of mammals. By affecting bone ossification the locus potentially can impact allometric size variation within population, which in turn would usually impact mostly facial traits due to its late development (see 61,62 for craniofacial development in mammals). Consequently it would lead to an increase or decrease in correlation among facial traits and thus to a larger or smaller morphological integration. Thus, the strong positive correlation found between evolutionary flexibility index and the RUNX2 gene (Fig. 3B) suggests a potential association between the high potential ability for craniofacial traits evolution in phyllostomid bats and RUNX2 variation. The increased ability to promote bone tissues formation might have played a fundamental role during transitions to feeding specializations, thus facilitating evolutionary change. Finally, since the ratio of glutamine to alanine residues in the RUNX2 protein might influence the regulation of bone development, our results indicate that Q/A ratio can act as a key feature for natural selection to operate in these species.

Long, short, or broad faces: what is the impact of RUNX2?
Our results clearly showed contrasting patterns of association between tandem repeats of RUNX2 Q/A and palate shape in the highly diverse groups of phyllostomid bats and primates (Tables 1, 2 and Fig. 2). Until recently, studies on marsupials and placental mammals, including domestic dog breeds, evidenced two contrasting patterns: i) a correlation of Q/A ratios and face lengths 10 ; or ii) a non-association of these traits 11,12 . Sears et al. 10 suggested a pattern of increase in face length of carnivores due to a higher Q/A ratio. In contrast, the results obtained in our study point to a reduction in length and a widening of the maxilla, which would lead to shorter and broader faces in phyllostomid bats. However, in NWM, the result is consistent with that described for the order Carnivora. Overall, our results support the main idea proposed by Sears et al. 10 , reinforcing the role of the Q/A ratio as a flexible genetic mechanism that would rapidly change the time of ossification of the skull throughout development.
In this context, the findings of the present study propose a complex scenario, where the permissiveness of this genetic adjustment system created by the variability in the Q/A ratios of the RUNX2 gene is indirect, and does not lead to a specific phenotype. However, throughout the history of a lineage (at the species or family level), RUNX2 variability acts with external factors and evolutionary constraints, together with other genes, as a facilitator for the adaptive morphological changes currently found. In the case of phyllostomids, these changes seem to have been selected by dietary constraints, corroborating the phenotype-environment correlation 27,32,33,41,63 . Overall, it should be noted that many genetic, epigenetic, and environmental attributes modulates the expression of genetic variance and interact with developmental networks to produce the available phenotypic variation. Therefore, the  Table 3. Results of phylogenetic generalized least squares (PGLS) models exploring the putative association between RUNX2 Q/A ratios and relative palate-length and width in phyllostomid bats and primates, depicting the magnitude of phylogenetic signal (Pagel's lambda). results of the present study seem to agree with Alberch's 64 view on the conceptualisation of genotype-phenotype association, i.e. genes do not directly specify the development or the shape of an organism; they are only one of the several causal factors that jointly determine the phenotype. Thus, developmental events are able to ambiguously affecting, or be affected, by gene expression.

Material and Methods
RUNX2 sequence data. Samples from 31 phyllostomid bat species from a wide range of feeding habits and craniofacial morphologies were surveyed in this study (electronic supplementary material, Table S1). Tissues and skulls individually paired are deposited in the Department of Zoology collection at the Fundação Universidade Regional de Blumenau (FURB Qiagen, Germantown, MD, USA) was used to extract and isolate genomic DNA from muscle samples for the 37 species. The polymerase chain reaction (PCR) was carried out to amplify the target RUNX2 exon, using primers and conditions described by Pointer et al. 11 , standardising the PCR conditions for the different species surveyed. Shrimp alkaline phosphatase and exonuclease I enzymes (GE Healthcare, Chicago, IL, USA) were used to purify the PCR products. Amplicons were sequenced in an ABI3730xl (Applied Biosystems Inc., Foster City, CA, USA) sequencer, using the same primers employed in the PCR for both DNA strands. Sequences were analysed using the ClustalW algorithm, implemented in the Codon Code Aligner software (CodonCode Corporation, Dedham, MA, USA) 65 . A prior screening test was performed, which considered only the sequences presenting conserved regions flanking the Q/A motif. These regions were then used as reference markers for the Q/A motif. The Q/A ratio for each specimen was calculated by dividing the number of repeated glutamine residues by the number of repeated alanine residues in tandem 10,11 (Table 1).
Samples and linear measurements. We measured 31 phyllostomid skulls (the same specimens sequenced for RUNX2) using a digital calliper, and obtained the primate data set employing a combination of craniometric data recorded by Marroig and Cheverud 66 and Oliveira et al. 54 using a Microscribe 3DMX. From these data, we obtained five linear measurements from the ventral crania of all specimens: palate length, palate width, skull length, skull width, and zygomatic breadth (Fig. 4, see also the electronic supplementary material, Table S2, for detailed measurement descriptions). These linear measurements were chosen to reflect variations in skull form among species and to represent the dimensions of palate width and length both in phyllostomid bats and primates 54,63,66,67 .
Additionally, evolutionary alterations in palate morphology are thought to underlie much of the feeding specialisations in bats 50 . We calculated the geometric mean for all measurements and standardised the palate length and width for each specimen by dividing each measurement by the geometric mean (GM) ( Table 1). Only skulls from adult specimens were measured. For both primates and bats, specimens were considered adults when the basisphenoid and basioccipital joint was completely fused.

Statistical analysis.
In an initial set of analysis, we performed Spearman's correlation test to investigate the association between RUNX2 tandem repeat ratios and palate length and width. Additionally to this traditional statistical method, a Bayesian phylogenetic comparative method (PCM) approach was implemented in CoEvol 1.4 software 68 to investigate the correlation between RUNX2 Q/A ratio and palate length and width, taking into account phylogenetic relationships, evolutionary rates, and divergence time under a Brownian motion model of evolution. In the CoEvol results, a posterior probability (pp) close to 1 represents a strong statistical support for a positive correlation, and a pp close to 0 indicates a supported negative correlation. Two data sub-sets for each taxonomic group (bats and primates) were analysed: I, morphometric measures for each species, corrected by geometric mean of all measures and II, phylogenetic relationships estimated from neutral markers (i.e., mitochondrial genes obtained from public databases). We ran CoEvol for 30,000 MCMC (total number of cycles) and considered the analysis finished when all effective sample size values reached at least 300 68 .
Finally, we performed a phylogenetic generalised linear model (PGLS) using the R function (pgls) of the "caper" R package 69 to investigate the magnitude of phylogenetic signal associated with phenotypic (palate length and width -dependent variables) and genotypic (Q/A region -independent variable) changes. This function fits a linear model controlling for phylogenetic non-independence between data. We also used this parameter for comparative analyses with the results of previous studies on carnivores 10 and other placental mammals 11 . We assume that all the three analytical methods employed in our study are integrative and represents different levels of complexity, ranging from a more general exploratory analysis (Spearman's correlation) to a more complex phylogenetic comparative methods (PGLS and Bayesian).

Q/A ratios and morphological integration index.
The magnitudes of integration among traits were obtained by using the coefficient of variation of the eigenvalues (ICV 70 ;). This index was calculated for primates 56,59 and phyllostomids (Rossoni et al. manuscript in preparation) as the standard deviation of the eigenvalues (σ (λ) ) divided by the average of the eigenvalues (λ). The flexibility and constraints index were also provided from the author's cited above. The evolutionary flexibility index captures the ability of a species to respond in direction of selection, while the evolutionary constraints index captures the relative influence of the line of least resistance (first principal component) on the direction of the evolutionary responses 59 . These measurements are presented in the electronic supplementary material, Table S3 (see also 58,59,70 for a detailed description of these index estimates). We performed Spearman's correlation test in order to investigate the association between RUNX2 tandem repeat ratios and ICV, flexibility and constraints index.