Phenotypic innovation in one tooth induced concerted developmental evolution in another

Serial appendages are similar organs found at different places in the body, such as fore/hindlimbs or different teeth. They are bound to develop with the same pleiotropic genes, apart from identity genes. These identity genes have logically been implicated in cases where a single appendage evolved a drastically new shape while the other retained an ancestral shape, by enabling developmental changes specifically in one organ. Here, we showed that independent evolution involved developmental changes happening in both organs, in two well characterized model systems. Mouse upper molars evolved a new dental plan with two more cusps on the lingual side, while the lower molar kept a much more ancestral morphology, as did the molars of hamster, our control species. We obtained quantitative timelines of cusp formation and corresponding transcriptomic timeseries in the 4 molars. We found that a molecular and morphogenetic identity of lower and upper molars predated the mouse and hamster divergence and likely facilitated the independent evolution of molar’s lingual side in the mouse lineage. We found 3 morphogenetic changes which could combine to cause the supplementary cusps in the upper molar and a candidate gene, Bmper. Unexpectedly given its milder morphological divergence, we observed extensive changes in mouse lower molar development. Its transcriptomic profiles diverged as much as, and co-evolved extensively with, those of the upper molar. Consistent with the transcriptomic quantifications, two out of the three morphogenetic changes also impacted lower molar development. Moving to limbs, we show the drastic evolution of the bat wing also involved gene expression co-evolution and a combination of specific and pleiotropic changes. Independent morphological innovation in one organ therefore involves concerted developmental evolution of the other organ. This is facilitated by evolutionary flexibility of its development, a phenomenon known as Developmental System Drift. AUTHOR SUMMARY Serial organs, such as the different wings of an insect or the different limbs or teeth of a vertebrate, can develop into drastically different shapes due to the position-specific expression of so-called “identity” genes. Often during evolution, one organ evolves a new shape while another retains a conserved shape. It was thought that identity genes were responsible for these cases of independent evolution, by enabling developmental changes specifically in one organ. Here, we showed that developmental changes evolved in both organs to enable the independent evolution of the upper molar in mice and the wing in bats. In the organ with the new shape, several developmental changes combine. In the organ with the conserved shape, part of these developmental changes are seen as well. This modifies the development but is not sufficient to drastically change the phenotype, a phenomenon known as “Developmental System Drift”, DSD. Thus, the independent evolution of one organ relies on concerted molecular changes, which will contribute to adaptation in one organ and be no more than DSD in another organ. This concerted evolution could apply more generally to very different body parts and explain previous observations on gene expression evolution.


INTRODUCTION
6￼ this remains generally untested. Together this indicates that the independent evolution of 140 serial organs does cope with pleiotropic genes, but does not fully explain how. 141

142
Here we chose a new model and a different approach to address this question, focusing on 143 developmental dynamics. We studied the independent evolution of a drastically new shape 144 in the mouse upper molar, which is nicely described in the fossil record. Since it occurred 145 relatively recently, we could compare closely related species. Another strong advantage of 146 this model is that the development of mouse molars is very well understood 147 mechanistically, from years of developmental genetics and morphogenesis modeling.  Davis, 2011;Hillson, 2005). Hamsters are today's good 161 representative of the basal "cricetine" rodents. In the golden hamster lineage, both molars 162 kept the ancestral cusp number and organization. We can make the reasonable assumption 163 that the hamster presents ancestral developmental features, and in this study we use this 164 species as a phylogenetic control.    Mouse and hamster molars develop at different paces. We predicted developmental age 268 from embryonic weight in each species and aligned temporal series between species with 269 homologous start and end points of first molar morphogenesis (Figures 1 and S1). We then To detect genes carrying this conserved variation between upper and lower molars, we 298 modelled their temporal profiles with polynomials ( Figure 2A). In each species, we fitted 299 two models: one with two distinct curves, one per tooth, and another with a single curve, 300 common for both teeth. By comparing the fit of these two models, we detected 301 differentially expressed genes between the two teeth (   show an epithelial bulge (arrowhead) never seen in the lower molar.

332
The set of 550 genes with significant and consistent upper/lower bias in both species is 333 highly relevant: it contains the expected jaw-identity genes known for mouse (Nkx2.3, 334 Pou3f3, Dlx1, but also Dlx5 and 6 which are no longer lower-jaw specific at this stage, but 335 show a lower molar bias as already described in mouse (Pantalacci et al., 2017)), a fifth of 336 the genes whose mutant shows a phenotype in the lower molar (21 out of 87 "keystone 337 genes" from (Hallikas et al., 2021)) and key transcriptional regulators of molar 338 morphogenesis (Barx1, Msx1, Pitx1 Figure S2). Overall, these genes are strongly enriched 339 for transcriptional regulators. They are involved in epithelial and mesenchymal 340 development, cell adhesion, proliferation and differentiation and signalling, especially 341 WNT, BMP and NOTCH (enrichment for Gene Ontology terms, Figure S3). Among the 550 342 genes with a consistent upper/lower bias, only 165 display evolutionary conserved temporal 343 dynamics in both teeth (stars Figure S2, Table S2). 344

15￼ 346
The earliest phase of cusp patterning gathers many features that distinguish lower and 347 upper molars in both species 348 349 We next looked for criteria that would distinguish cusp patterning dynamics in lower and 350 upper molar in both species, and may have been conserved from the common ancestor of 351 mouse and hamster. 352 353 First, we noticed that in both species, SEK formation is initiated (transition to 1-SEK stage 354 and 2-SEK stage) and completed (transition to 6-SEK stage) later in the upper molar as 355 compared to the lower molar of the same embryo. This delay was also obvious in the 356 transcriptome: the upper molar transcriptome looks "younger" than the lower molar 357 transcriptome of the same embryo ( Figure S4).

400
Example of Sfrp2 gene expression levels (grey dots), significantly better modelled by two curves in 401 mouse and one in hamster (purple curves), compared with one curve per species (grey curves). Right:

402
Number of genes detected by this tooth-specific model in mouse (purple) and in hamster (green, LRT 403 with adjusted p < 0.05). Barplots show their frequencies for the "total" gene set, genes from 404 developmental "pathways", genes from "bite-it" database, and genes with a mild or strong 405 phenotype in tooth mutants ("Dispensable" or "Keystone"). Size of each gene set into brackets.

18￼
Pou3f3 is the only TF specific of the upper molar in both species. Its expression showed a 407 twofold increase in mouse upper molar. Dlx1/2 genes are expressed in both molars, but are 408 essential only for upper molar formation in mouse (Qiu et al., 1997). Their expression levels 409 were more than twice increased in mouse molars and the ratio, slightly in favor of the upper 410 molar in hamster, is increased in mouse. Barx1 is a key molar-specific TF whose levels have 411 been correlated with cusp number in mammalian molars (Miletich et al., 2011). The bias in 412 favor of the upper molar was markedly increased in mouse. This effect is selective since 413 Msx1, another TF which cooperates with Barx1 (Miletich et al., 2011), showed similar bias in 414 the two species ( Figure S2). 415 416 Surprisingly, such changes were not restricted to the upper molar. The expression of 417 Nkx2.3, the specific TF of the lower molar, showed an almost twofold increase in mouse. 418 For Pitx1, a shared TF whose mutation impairs more specifically lower molar development 419 (Mitsiadis & Drouin, 2008), the ancestral bias in favor of the lower molar was increased. Dlx1 420 expression levels were also twice increased in the lower molar. This is not true for the Dlx5-421 6 genes, which were more specifically associated with lower jaw identity (Depew et al., The dissimilarity of upper/lower molar transcriptomes is increased in mouse 427 Next, we asked whether the transcriptomes would capture a general increase of 428 differentiation of molar development in mouse as compared to hamster. 429 As exemplified above, a number of genes have a marginally or small significant bias in 430 hamster, which is increased in mouse. This is in agreement with the multivariate analysis 431 presented in Figure 2C. Along an axis that captures a lower/upper molar variation common 432 to both species, the variation is higher for mouse than for hamster. Thus, an ancestral 433 dissimilarity is exaggerated in mouse. 434

19￼
We detected genes whose expression profiles differ between upper and lower molars in 436 mouse but not in hamster, by a dedicated model based on the 4 molars altogether ( Figure  437 3B). As a control, we built a reciprocal model to detect genes with tooth specific profiles in 438 the hamster. We found 2.5 times more genes with tooth-specific profiles in mouse as 439 compared to hamster. Even after removing the effect of baseline expression levels, which 440 may be impacted by differences in cell composition (Pantalacci et al., 2017), we still 441 observed 1.6 times more tooth-specific profiles in mouse. The function of the genes with 442 tooth specific profiles differs markedly between mouse and hamster ( Figure S6). In mouse, 443 genes are linked to cell adhesion and migration, as well as cell cycle and mitosis. This is 444 consistent with our previous findings suggesting that morphogenetic movements are 445 We previously showed that the upper molar germ of the mouse contains more 461 mesenchyme relative to the epithelium than the lower germ since early cap stage. Tooth 462 engineering studies suggest that this higher proportion may help to form the 463 supplementary cusps. Indeed, in artificial teeth made by reassociating a varying amount of 20￼ mesenchymal cells to a single epithelium, the number of cusps formed increases with the 465 number of mesenchymal cells (Hu et al., 2006). To control whether this higher 466 mesenchyme:epithelium ratio is specific to the mouse, we extracted mesenchyme and 467 epithelium-specific marker genes from tissue-specific transcriptomes ( Figure 1C), and used 468 in silico deconvolution to estimate the mesenchyme proportions from whole tooth germ 469 transcriptomes ( Figure 4A). The proportion of mesenchyme was indeed significantly higher 470 in the upper molar in mouse, but not in hamster (Wilcoxon tests, p < 2e-16 and p = 0.152).     To look for evidence of persistent molecular B/L polarity during cusp morphogenesis, we 516 collected and analyzed mouse transcriptomes of buccal and lingual halves at the early 1-517 SEK stage ( Figure 1C). We found that Osr2 and Sfrp2 are still expressed with a strong 518 lingual bias at 1-SEK stage, and in the timeseries, their expression is maintained at high 519 levels during the period of bucco-lingual development of the tooth germ ( Figure S7). To 520 determine if the tooth germ is polarized at the 1-SEK stage, we estimated the levels of molars, but very weakly on the lingual halves, which therefore still appears as a naive tissue 528 ( Figure 4C). In the upper molar, the lingual half looked even more naive than in the lower 529 molar, with even lower levels of pathway activities, consistent with twice higher levels of 530 Sfrp2 expression, and delayed downregulation ( Figures 4C, S2, S3). In the buccal half, 531 activation of the BMP4 and WNT pathways in the mesenchyme is stronger in the upper 532 molar, which thus appears as more polarized than the lower molar. These findings support 533 the idea that the BMP4/OSR2 antagonism is still acting during early mouse molar 534 morphogenesis to set up the B/L polarity of the tooth. This polarity maintains naive tissue 535 on the lingual side of the germ at 15.0, which grows faster and shows delayed cusp 536 formation relative to the buccal side. 537

538
We reasoned that the larger the proportion of this naive lingual tissue, the stronger the 539 germ growth potential and its capacity to form cusps on its lingual side. We therefore 540 quantified this proportion in mouse and hamster tooth germs, by deconvoluting the 541 timeseries dataset with buccal and lingual tissue transcriptomes (Gong & Szustakowski, 542 2013). As expected due to progressive cusp formation, we found that the proportion of 543 naive lingual tissue decreases during morphogenesis in both species ( Figure 4D). But in 544 mouse molars, and even more markedly in the mouse upper molar, the initial proportion of 545 24￼ naive tissue is larger, and diminishes more slowly. We noted that the naive tissue proportion 546 correlates with levels and temporal pattern of Sfrp2 expression in mouse and hamster 547 molars ( Figure S7). 548

549
In summary, a B/L polarity of the tooth germ is maintained during the first steps of cusp 550 formation. This is especially true in mouse molars whose proportion of lingual naive tissue is 551 increased, and correlates with higher levels of Sfrp2 expression. This polarization is further 552 exaggerated in the mouse upper molar.

26￼ 594
To gain insight into the function of Bmper, we obtained a mouse null mutant and studied 595 the shape of its molars. Since the homozygous Bmper mutant are lethal at birth, this had to 596 be done by reconstructing the enamel-dentin junction at 19.5 days, carefully matching them 597 with controls of similar developmental age (see material and methods). The upper molar is 598 modified: one of the buccal cusp is absent or poorly grown ( Figure 5C, arrowheads) and 599 one lingual supplementary cusp is more prominent (arrows). The lower molar is also 600 modified: the central lingual cusp, which is determined by formation of the 2-SEK following 601 lingual growth, is more prominent. Thus, the loss of Bmper modifies the bucco-lingual 602 equilibrium, favoring the lingual side of the molar at the expense of the buccal side. Since 603 Bmper modulates the BMP4 pathway, this mutant phenotype reinforces the idea that the 604 BMP4 pathway regulates B/L polarity during cusp formation, and that Bmper might have a 605 causative role in displacing the BMP4/OSR2 balance in mouse. 606

607
In summary, we show that the BMP4/OSR2 antagonism persists in mouse to regulate the 608 B/L polarity of the tooth during morphogenesis. As compared to hamster upper molar, the 609 mouse upper molar has an asymmetrical expression of Bmper (buccal side) and an 610 increased expression of Sfrp2 (lingual side). This is associated with an increased and 611 persistent proportion of naïve lingual tissue. These differences seem very consistent with 612 the newly evolved lingual cusps of this tooth. We were very intrigued however that 613 qualitatively, these morphogenetic changes are also seen in the lower molar, although to a 614 lesser extent and without change in number and respective size of buccal/lingual cusps. We     (Hallikas et al., 2021). A-C, Gene categories as in Figure 2A with 654 numbers into brackets, genes from developmental "pathways" further splitted. We first examined Bmp4, since finding this essential gene for tooth development among 669 co-evolving genes was a surprise ( Figure 6B). Bmp4 is expressed in both the epithelium and 670 the mesenchyme, but the epithelial domain is so small relative to the mesenchymal domain 671 (see later Figure 7 and S9), that the latter will dictate the bulk transcriptomic profile. We Wif1 rises earlier in mouse transcriptomes, and its mesenchymal expression is seen both 682 earlier and in a larger territory in mouse tooth germs (Figures 6C, S9). Dkk1 expression 683 transiently decreases in mouse transcriptomes, which coincides with an earlier relocalisation 684 of its expression at future cusp tips beneath the SEK ( Figure S9). This expression pattern 685 suggested to us that cusp formation might in fact be anticipated in mouse as compared to 686 hamster. 687 We thus turned back to our quantification of cusp formation dynamics, and realized that 688 both mouse molars quickly transition to 1-SEK after a rather short PEK stage (Figures 1, 7A). Since the lower molar phenotype has been much more conserved during evolution, the 752 lower molar developmental phenotype captured by the temporal profiles should be more 753

conserved. Otherwise, this is an indication of DSD. 754 755
We scored the divergence between mouse and hamster upper and lower molars by 756 modelling temporal profiles with polynomials (LRT with adjusted p < 0.05). We found that 757 for 22.0% of genes, the profiles diverged in the lower molar, which is even more than in the 758 upper molar (17.3%, Figure 6D, Table S2). This is true as well for genes relevant for tooth 759 development and phenotype ("bite-it", "keystone", "pathways"; Figure 6E). 760 Put together, these observations suggest that the development of the lower molar has 761 drifted while co-evolving with the upper molar.

35￼
As in mouse molars, co-evolution is pervasive in bats limbs and concerns genes whose 800 expression evolution was key for the independent phenotypic evolution of the forelimb. study (Pantalacci et al., 2017)). Since these findings were made in hamster and mouse, 820 these specificities of lower and upper molar were likely present in their common ancestor, 821 but we suspected they may even date from early mammals. molars. This also constitutes a case of "recapitulation" since early ontogeny of cusp 838 formation recapitulates phylogeny (Gould, 1977). to be essential to tooth development, such as Msx1, Barx1, Pitx1 also showed a conserved 861

bias. It remains to be tested if this bias is directly controlled by identity genes. 862
Our results show that during evolution, the details of developmental interactions in serial 863 organs diverge extensively, but some developmental specificities of one organ with respect 864 to the other are conserved (e.g. the relative order and timing of appearance of the 3 first 865 cusps and the period of maximal transcriptomic divergence, the delayed development of 866 the upper molar). These specificities could be encoded in a conserved relative dose of the 867 key transcription factors specifying an organ (here a molar or a limb), and this conserved 868 relative dose could be controlled by identity genes. Altogether this forms an ancestral 869 molecular and developmental identity for the two teeth. 870 871

Reinforcement of the ancestral molecular and morphogenetic identity in mouse molars 872
We found that the molecular identity was reinforced in mouse, in the upper molar but also 873 more surprisingly in lower molar. Expression levels doubled for the upper-molar specific TF 874 Pou3f3 and the lower-molar specific TF Nkx2-3, and the ancestral bias of many genes was 875 exaggerated, whether in favor of upper or lower molar (e.g. Barx1 and Dlx1; Pitx1, 876 respectively). Consistent with these changes in individual TFs, we found at genome-wide 877 scale that the temporal profiles in the mouse were exaggeratedly different from each other. 878 As for specific TF, divergence is seen both in upper and lower molar. "recapitulation" looks superficially like a case of "terminal addition", a mechanism for 885 evolutionary change whereby development is incremented with one more step (Gould, 886 1977). However, species divergence peaks early in development and we point to three 887 features that concern mouse early upper molar development, but could pave the way for 888 the supplementary cusps (Figure 8). The second feature is the stronger polarisation of the upper molar field along the 918 bucco-lingual axis, associated with a precocious transition from PEK to 1-SEK and a very 919 long 1-SEK stage. As a consequence, a larger undifferentiated field is present on the lingual 920 side, where activation-inhibition loops can pattern SEKs. This feature seems highly relevant 921 because it could explain why the increase in cusp number is focused on the lingual side of 922 the tooth, while the size of buccal cusps is reduced. It seems to exploit an ancestral 923 specificity of the upper molar as compared to the lower molar, which produces a longer 1-924 SEK stage. This specificity could combine together with a novelty in mouse responsible for 925 shortening the PEK stage in both teeth, and produce an upper molar with a very long 1-SEK 926 stage and a large undifferentiated field, while change remains more modest in the lower 927

molar. 928
The third feature is the narrower range of expression of signaling molecules in 929 mouse signalling centres, which is especially obvious for Bmp4. Reducing the range of Mapping mutations corresponding to these developmental phenotypes is out of the 937 scope of this study, but transcriptomics provided us with molecular mechanisms and some 938 candidate genes. No expression change was observed in two obvious candidates from the 939 literature (Fgf3 and Activin A, Figure S12, (Charles et al., 2009;Kwon et al., 2017). 940

941
In insects, the evolution of the dose of identity genes has been correlated with the 942 evolution of the size of serial organs (Paul et al., 2021). Pou3f3 is expressed specifically in 40￼ the upper molar mesenchyme of both species and its dose is twice increased in mouse. 944 Hence Pou3f3 is a good candidate to explain the larger proportion of mesenchyme 945 specifically in the mouse upper molar. Pou3f3 loss-of-function mutants miss some skeletal 946 elements of the upper jaw, but their upper molars showed "no major defects" (data 947 unshown in (Jeong et al., 2008)). This may deserve re-examination, or study in sensitized 948 backgrounds. The causative mutation may also be upstream in the regulatory network, in 949 particular in Dlx1/2, since Pou3f3 is regulated by Dlx1/2 in the early jaw (Jeong et al., 2008)  The three features that we observed hint at very complementary aspects of tooth 1037 development (Figure 8). The tooth literature shows it is difficult to increase cusp number in 1038 mouse molars: in vitro experiments have shown it can be necessary to play on multiple 1039 pathways, and mouse mutants show, at best, small accessory cusps, but no supplementary 1040 main cusps (Harjunmaa et al., 2012). None of these three changes should be sufficient on its 1041 own to induce the major changes in cusp size and proportions seen in mouse as compared 1042 to its ancestor. We therefore propose that the new phenotype involves combining 1043 mutations in at least two or three different genes, corresponding to these three features 1044  interdigital area in the wing, but not in the foot (from Figure 3H in (Weatherbee et al., 1116(Weatherbee et al., 2006). In contrast, the BMP inhibitor Grem1 expression persists in both limbs, with higher 1117 levels in the wing (from their Figure 3D, note the blue staining remaining around digits 1118 whereas more proximal parts of the limb are not stained at all). Thus, at this stage, specific 1119 and exaggerated shared gene expression changes seem to combine to pass the threshold 1120 for apoptosis suppression in the wing, but not in the foot (Figure 8). This evolutionary 1121 scenario of independent evolution is thus very similar to teeth, involving a combination of 46￼ specific, shared, and exaggerated shared expression changes and differential threshold 1123 shared effect could be used by adaptation. This is unexpected since it is commonly thought 1148 that adaptive mutations need to be modular at the DNA level to have organ-specific effects 1149 and thereby circumvent gene pleiotropy. The capacity of developmental systems to 1150 undergo DSD is another way of circumventing gene pleiotropy, and thus appears as a 1151 mechanism by which non-modular mutations can be selected in adaptation. We propose 47￼ this is the reason why independent evolution can be so frequently seen in nature despite 1153 gene pleiotropy. 1154 1155 Pleiotropy, concerted evolution and DSD 1156 Serial organs such as molars and limbs have a heavy pleiotropy load and for this reason, 1157 they are possibly especially prone to developmental co-evolution and DSD. We 1158 nevertheless believe that our results in serial organs illustrate a much more general 1159 correlation between pleiotropy and DSD at the organismal level, as suggested previously 1160 (Félix, 2007;Pavlicev & Wagner, 2012). 1161 The link between pleiotropy and DSD has been observed in experiments of in silico 1162 evolution (Johnson & Porter, 2007;Tulchinsky et al., 2014). It has also been observed in 1163 nematode genetics with a mutation increasing the fitness in laboratory conditions that has 1164 induced DSD in the vulva (Duveau & Félix, 2012). Finally, a link between pleiotropy and 1165 concerted transcriptomic evolution has already been suggested. In most multispecies 1166 transcriptomic analyses, samples of different organs tend to group by species (like molar 1167 samples in Figure 1C). This pattern, so-called "species signal", often dominates in samples