Mary-Anne Maverick

Mathieu Noël

Jean-Jacques Vasseur

Carine Baraguey

Françoise Debart
email: francoise.debart@umontpellier.fr

Michael Smietana
email: michael.smietana@umontpellier.fr

Hypervalent Iodine-Mediated Preparation of 5-Hydroxuridine from 2',3'-O-isopropylidene N4-Acetylcytidine

Keywords: Hypervalent iodine, nucleoside, cyclouridine, oxidative cyclization, 5-hydroxiuridine

Whereas BAIB-TEMPO oxidation of 2',3'-O-TBDMS-N4-acetylcytidine results in the expected 5'-carboxylic acid nucleoside, its 2',3'-Oisopropylidene analogue reacts in a radically different way. We have demonstrated here that hypervalent iodine(III) in water triggers an unprecedented oxidative cyclization leading to a mixture of C5-substituted O6,5'-cyclo 5,6-dihydro uridines. This mixture of cyclouridines can be opened under basic conditions and, after deprotection, yields 5-hydroxyuridine, an important post-transcriptional modification of uridine at the wobble position (U34) of bacterial tRNA. NMR experimental values and calculations were performed to provide further insight on the specific reactivity of 2',3'-Oisopropylidene N4-acetylcytidine.

Introduction

5-Hydroxyuridine (5-OHU) is a major RNA lesion affecting ribosomal translation that can be generated from pyrimidines by the action of reactive oxygen species and free radicals. 1 Interestingly, 5-OHU is also involved in the formation of 5oxyalkylated modifications of the wobble position of bacterial tRNA that are believed to play key roles in the translational machinery. 2 From a synthetic point of view, the most common method to synthesize 5-OHU involves treating uridine with a halogenation reagent in basic aqueous media. 3 Based on their work on uridine 5'-diphosphate, Eaton and Hutchinson suggested that 5-hydroxy pyrimidine nucleoside was generated from a bromohydrin intermediate in which the bromine atom at C5 position was subsequently substituted by a hydroxyl group, followed by the elimination of a water molecule. 4 It is worth noting that, due to the presence of the 5'-diphosphate group, this proposed mechanism eliminates the possibility of the formation of an O6,5'-cyclouridine intermediate usually generated by anchimeric assistance of the 5′-hydroxyl function. 5,6 These authors also demonstrated that the corresponding 5-bromopyrimidine nucleoside was not an intermediate in the reaction and that when performed on a cytidine analogue, deamination occurred to give 5-OHU as the major product, although in modest yield (34%). 7 The scope of enzymatic reactions known to modify cytosine by deamination, methylation or oxidation is such that researchers see cytosine as a genomic "wild card" able to impact coding sequences and genetic expression. 8 To the best of our knowledge, the only mention of oxidative transformation of cytosine to 5acetoxyuridine was reported by Maki and co-workers in 1984. 9 First, the treatment of 2',3'-O-isopropylidene N4benzoylcytidine 1a with lead tetraacetate generated a mixture of 5'-O,6-cyclo-2',3'-O-isopropylidene N4-benzoylcytidine 2 (15%) and 5-acetoxy-5'-O,6-cyclo-5,6-dihydro-2',3'-Oisopropylidene uridine 3 (50%), by oxidative cyclization (Scheme 1a). Heating 3 in dimethylformamide for an unspecified time gave 5-acetoxy-2',3'-O-isopropylidene uridine 4, although the obtained yield was not reported. In the course of our study dedicated to the synthesis of C-5'-acyl-modified nucleosides, 10 we had to prepare several 5'-carboxylic acid nucleosides 11 following the efficient 2,2,6,6-tetramethyl-1piperidinyloxyl (TEMPO)-[bis(acetoxy)-iodo]benzene (BAIB) method described by Epp and Widlanski. 12 C-5'-oxidation proceeded smoothly for all studied pyrimidine nucleosides, except for the cytidine protected by N4-acetyl-2',3'-Oisopropylidene groups, which led to a complex mixture of compounds. This observation was all the more surprising since, under the same conditions, N4-acetylcytidine 2',3'-O-protected by silyl groups gave the expected 5'-carboxylic acid cytidine. We have now demonstrated that in the case of 2',3'-Oisopropylidene N4-acetylcytidine, a mixture of two nonchromogenic cyclouridines 3 and 5 is formed during the hypervalent iodine-mediated reaction and that both compounds can be readily converted to 5-OHU (Scheme 1b). Here are the results of our investigations.

Results

Early studies showed that under standard oxidation conditions (BAIB 2.2 eq, TEMPO 0.2 eq, acetonitrile-water (1:1)), 2',3'-Oisopopylideneuridine was oxidized to 5'-COOH as described in 80% yield. Under the same conditions, 2',3'-O-isopopylidene-N4-acetylcytidine led mainly to a mixture of UV-inactive compounds that could only be detected by vanillin staining, with only traces of the expected carboxylic acid derivative. This Template for SYNLETT Thieme observation suggests that the aromaticity of the nucleobase was lost during the reaction.

Scheme 1. Synthetic approaches toward metal-free oxidative cyclization of N4-protected 2',3'-O-isopropylidene cytidine Two major compounds, 3 and 5, were separated and isolated by column chromatography in 47 and 4% yield respectively, and were analyzed by NMR to unequivocally determine their structures. Analysis of the 1 H and 13 C NMR spectra of 3 clearly indicates a change in the structure of the nucleobase. First of all, the two doublets belonging to the vicinal H5 and H6 protons of the base were found to be linked to carbon atoms with chemical shifts too low (70.1 and 87.3 ppm) to be sp2-hybridized atoms. Moreover, heteronuclear correlations were observed between C6 and the methylene protons in the 5' position of the ribose. These data support a break in the π-bond between C5 and C6 and a sugar-base cyclization involving the 5'-extremity and the C6 carbon, as depicted in Scheme 1. Such a change in the hybridization of C5 and C6 leads to several stereoisomers and through-space correlations were used to discriminate between the different possibilities. The R configuration at C6 would bring H6 closer to H1', H2' and H3', but the corresponding correlations were not observed on the ROESY spectrum. Instead, H6 and H3' showed correlations peaks each engaging a different proton of the methylene (H6-H5' and H3'-H5''). This enabled these protons to be assigned stereo-specifically (H5' and H5'' being Pro(S) and Pro(R) respectively) and, above all, proved that H6 and H3' are opposite the mid-plane of the ring formed by the O5'-C6 bridge, thus confirming the absolute S configuration of C6 (Figure 1). Compound 5 showed very similar 1 H and 13 C spectra to compound 3, except that the acetyl group on the base was missing. Consequently, the most important variations in chemical shifts when comparing the spectra of the two compounds concerned the H5 and H6 protons, as well as the C4 and C6 carbons. The NMR data indicate that 5 is a deacetylated analogue of 3, bearing a hydroxyl group on C5. Bridging the 5'extremity to the C6 of the base pushes the nucleosides into a rigid conformation where the base and the dihedral angle around C4'-C5' are locked into anti and gauche-gauche conformations respectively, which in turn also constrains the sugar pucker. As dihedral angles and coupling constants are closely related, such conformational constrain impacts the coupling constants, and explains the simplified multiplicities in the 1H NMR spectra of 3 and 5, where three of the sugar's coupling constants are zero. Comparable values have been reported for 8,5'-anhydro-8-oxo-2',3'-O-isopropylidene adenosine (Table 1), 13 thus confirming that 3 and 5 share the ribose West conformation (P 240-280°) of this adenosine analogue.

On the basis of these structures, we started to evaluate the various parameters likely to influence the reaction. The addition of two equivalents of NaHCO3 to the solution as described by Epp and Widlanski 12 on unprotected 2',3'-O-isopropylidene cytidine, resulted in a sharp decrease in the yield of both compounds 3 and 5 (Table 2, entry 2). Interestingly, removal of TEMPO from solution increased the yields of 3 and 5 to 70 and 8% respectively (Table 2, entry 3). However, extending the reaction time to 12 h did not have a significant impact (Table 2, entry 4). The decrease or increase the number of BAIB equivalents led in both cases to a decrease in the formation of cyclodihydrouridines 3 and 5 in different proportions (Table 2, entries 5 and 6). Similarly, the removal of one of the two cosolvents had a detrimental effect on the production of 3 and 5 (Table 2, entries 7 and 8). However, nucleoside 5 was more abundant when the reaction was carried out in water only. Interestingly, when BAIB was replaced by phenyliodine(III) bis(trifluoroacetate) (PIFA), the yield of compound 5 raised to 55% with only 9% of derivative 3. No trace of an analogue bearing the trifluoroacetyl group in position 5 was observed. The slight presence of 3 suggests that the nitrogen protecting group may have been partially transferred to position 5 of the base, but the formation of a 5-trifluoroacetyl nucleoside intermediate, which would be rapidly hydrolyzed under these acidic conditions, cannot be ruled out (Table 2, entry 9). Finally, when 2',3'-O-isopopylidene N4-benzoylcytidine was subjected to BAIB, no transfer of the benzoyl group to position 5 was observed and compounds 3 and 5 were obtained in low yields, indicating that the main source of acetyl group to yield 3 comes from BAIB (Table 2, entry 10). adenosine [13] 0.0 5.6 0.0 2.3 0.7 -12.9

Template for SYNLETT Thieme a The coupling constant cannot be determined since H2' and H3' are accidentally isochronic. To gain more insight about the mechanism of formation of compounds 3 and 5, we protected the 5'-hydroxyl group of 1b with a benzoyl group and subjected the resulting compound to 2 equivalents of BAIB. As expected, no reaction occurred after 24 h and the starting material was fully recovered, confirming the crucial role of the 5'-hydroxyl group in the cyclization. Similarly, in the presence of BAIB alone, 2',3'-O-isopropylidene uridine also failed to react, suggesting the necessary involvement of the NH-Ac function.

Based on the above-mentioned experiments, we propose a possible reaction mechanism for the formation of cyclouridines 3 and 5 (Scheme 2). In line with the mechanistic proposal made by Suryavanshi and co-workers 14 for the formation of αacyloxyketones from benzylamine, the reaction between 1b and iodobenzene diacetate would generate an ammonium intermediate A that would trigger the cyclization via an attack of the 5'-hydroxyl group at position C6 to give the common intermediate B. In path a, intramolecular acetylation at position C5 generates imine species C by removal of iodobenzene. C is then attacked at the alpha-position of the imine by an acetate anion to form a hemiaminal ether intermediate D, which is hydrolyzed in water to yield compound 3. In path b, B is directly attacked by water at position C6, triggering the removal of iodobenzene and acetic acid to form imine E, which undergoes the same reaction sequence as D to give compound 5.

Scheme 2. Mechanism proposal for the formation of 3 and 5

Compounds 3 and 5 proved stable for weeks when stored in powder form at room temperature, and even when heated in DMF at 80°C for several hours, no degradation or transformation was observed. This is not surprising considering the numerous pyrimidine cyclonucleosides that have already been described. [5] Usually prepared by halogenation, dehalogenation or under radical conditions, 5,6 oxidative cyclization of cytidine via hypervalent iodine has, to our knowledge, never been reported. Hypervalent iodine chemistry has gained considerable importance for the development of a wide range of organic reactions, 15 but its potential in aqueous media has not yet to be fully exploited. 16 The use of hypervalent iodine reagents in nucleoside chemistry is still mainly limited to glycosylation reactions. 17,18 To easily access to 5-OHU, a mixture of 3 and 5 was stirred for 24 h in a MeOH/H2O/Et3N mixture (5:1:1), resulting in the unique formation of product 6 in 60% yield after purification on silica gel. The 2',3'-O-isopropylidene protecting group was then removed under acidic conditions to yield 5-OHU in 90% yield (Scheme 3). Thus, under optimized conditions (Table 2, entry 4), 5-OHU was obtained from 1b in 43% yield in three steps.

However, one question remains unanswered. Indeed, when 2',3'-O-isopropylidene protection is replaced by 2',3'-O-TBDMS groups, BAIB-TEMPO conditions lead to the formation of the expected 5'-carboxylic acid cytidine without any trace of cyclonucleosides. This is also the case with the 2'-deoxy-3'-O-TBDMS-N4-acetylcytidine analogue. To gain further insight, we performed a thorough NMR analysis of the conformation of 1b and its 2'-3'-O-TBDMS analogue 1c, dissolved into a CD3CN/D2O mixture (9/1 v/v) in order to approximate synthesis conditions. The conformation of a nucleoside is globally defined by three factors: the base and 5'-extremity orientation as well as the sugar pucker.

The determination of the sugar geometry is usually defined according to the pseudorotation concept, 19 in which the five-Template for SYNLETT Thieme membered ring shape is totally described by just two parameters: the phase angle P (0 to 360°) and the pucker amplitude ɸm. Most of the nucleosides and nucleotides occupy two regions of the pseudorotation wheel in particular: the North, centered around P SX).

In 2',3'-OH ribonucleos(t)ides undergoing a classical northsouth equilibrium (with a pucker amplitude of 39°), the cisoïdal J2'-3' coupling constant and the sum of the two transoïdal couplings constants (J1'-2' + J3'-4') are found around 5.0 and 9.5 Hz respectively. An abnormally puckered ribose is related to a lower value of J2'-3' associated with a higher value of J1'-2' + J3'-4'. In contrast, the reverse observation is diagnostic of a flattened sugar. 19 For 2',3'-O-TBDMS-C Ac , J2'-3' is found around 4.2 Hz, i.e. about 1 Hz below the usual value of ~5.0 Hz observed for nonconstrained ribose, while the value for J1'-2' + J3'-4' (9.3 Hz) remained quite close to the standard value. In contrast, for the 2',3'-O-isopropylidene nucleoside 1b, the value for J2'-3' was found about 1.5 Hz above the ~5.0 Hz value, while J1'-2' + J3'-4' dropped down to 5.3 Hz, indicating a flattening of the ribose and/or a shift of the phase angle (P) values corresponding to less usual conformations.Both compounds also showed different behavior with respect to increasing temperature. For 2',3'-O-TBDMS-C Ac increasing the temperature results in an increase of J1'-2' (J  +0.6 Hz) and in a decrease of J3'-4' (J  -0.7 Hz). The variations observed between the two extreme temperatures suggest the existence of an equilibrium between two conformations. For compound 1b, the coupling constants are much less temperature-dependent than for the silylated analogue 1c, J2'-3' and J3'-4' increasing in a lesser extent (J 0.1 and 0.3 Hz respectively) while J1'-2' remains unchanged (Table 3).

A graphical means of exploring the sugar conformation, in particular when a two-state equilibrium occurs, is to plot J3'-4' as a function of J1'-2'. 22 Vicinal coupling constants were calculated from Karplus-type equations using the Altona parametrization defined for β-ribose. 23 Plots of J3'-4' versus J1'-2' calculated for P = 0 to 360° and for two values of ɸm, together with experimental values, are shown in Figure 2.

For 2',3'-O-TBDMS-C Ac 1c (Figure 2, left), the experimental points are aligned on a straight line which suggests an equilibrium between two conformations.

Figure 2. Plots of J1'-2' vs J3'-4'. Curves represent the dependence of J1'-2' and J3'-4' calculated for P values ranging from 0 to 360° and for two values of ɸm using the Karplus-type equation parametrized by Altona. [24] Dots are the experimental data. Top: for compound 1c with ɸm = 35° (■) and 45° (). Bottom: for compounds 1b, with ɸm = 20° (■) and 35° (). Relevant P values are reported on the curves.

The reverse dependence of the two experimental coupling constants with temperature (Table 3) is characteristic of a classical north-south equilibrium and the low value of J1'-2' as well as the position of the points on the figure indicate that the north conformer slightly prevails. Its population can be approximated at 55-63% from 10 x 3 J3'-4'. [19] As the low value of 3 J2'-3' (4.2 Hz) suggests an increase in pucker amplitude, theoretical curves were plotted for ɸm = 35° and 45°. From the intercept of the line with the two curves, PN and PS were estimated to -20° to +20° (C2'-exo to C3'-endo) and 160° to 190° (C3'exo to C2'-endo) respectively.

Crystallographic data for 2',3'-O-isopropylidene analogues revealed that 2',3'-cyclization leads to a reduced pucker amplitude for the ribose ring, with ɸm ranging between 20 and 35 degrees for solid-state non-planar structures. [20][21] Curves of J3'-4' against J1'-2' for 1b were therefore built for these two ɸm values (Figure 2, right). The position of the data points in the lower left part of the graph indicates a predominance of a South/West conformation. The experimental points are clustered in a shorter area, suggesting a much lower flexibility compared to the 2'-3'-OTBDMS analogue 1c, but since it was nevertheless possible to fit a straight line through the experimental data, we concluded that compound 1b was still involved in a two-state equilibrium. From the intersection of this straight line with the two curves, it emerges that the compound mainly adopts a West-type conformation (240° < PW < 270°) in equilibrium with a less favorable North-type

Template for SYNLETT Thieme conformation (20° < PN < 50°). Orientation of the 5' extremity was then assessed, from analysis of the coupling constants implying protons H5' and H5'' after stereospecific assignment of the latter (see Supporting Information). The  torsion angle about the C4'-C5' bond (O5'-C5'-C4'-C3') usually adopts three main staggered conformations g + (60°), t (180°) and g -(-60°). The populations of each rotamer was estimated from the experimental coupling constants 3 JH4'-H5' and 3 JH4'-H5'' as follow: 80 % of g + , 18% of t and 2% of g -for 2',3'-O-TBDMS-C Ac 1c; 53 % of g + , 35% of t and 12% of g -for 1b. The g + conformation, which is found in the rigid cyclized dihydro compounds, is the one that brings the 5'-O extremity closer to the base. Nevertheless, this conformation is much more represented in the silylated compound than in the isopropylidene analogue 1b, which is nevertheless the only one capable of forming the 5'-O-C6 bridge.

In both compounds, the 5'-O extremity is flexible enough to adopt any of the three staggered conformations which is probably more important than the ability of the 5'-OH to preorganize into the final g+ conformation.

Finally, orientation of the base relative to the sugar and defined by the  torsion angle (O4'-C1'-N1-C2) was also considered using the three-bond heteronuclear couplings between H1' and the carbons C2 and C6 (see Supporting Information), and showed as expected the prevalence of the anti conformation for both 1b and 1c. It is very likely that the conformational constraints induced by the isopropylidene group at the 2'-O-and 3'-O positions are responsible for the difference in chemical reactivity between this compound and the 2',3'-O-TBDMS analogue. The creation of a bond between the base and the 5' extremity implies that the O5' and C6 atoms can be close enough in space for the bond to be formed. The conformational study showed that the most important difference between 1c and 1b is the sugar pucker, the dioxolane ring fused at the 2'-O and 3'-O positions pushing the ribose towards a more rigid and preferred West conformation. Whatever the 2'-3' protecting groupisopropylidene or TBDMS -in both cases, the base and 5'-OH preferentially adopt the ideal conformation where 5'-OH points in the C1'-C2' direction and C6 points in the O4'-C4' direction. This might explain the observed stereoselectivity at C6 in anhydro (ou dihydro?) compound 3. Interestingly, the g+ conformation is more favored in the silylated compound 1c than in the isopropylidene analogue 1b, yet only the latter led to the formation of an anhydro nucleoside. Moreover, it is well known that the sugar conformation directly affects the base and the 5' extremity's orientations. This suggests that although in both compounds the base and the 5' extremity are able to place themselves in the right position, in 1c where a standard northsouth equilibrium takes place, O5' and C6 are probably too far apart while in 1b the north-west equilibrium would bring these atoms close enough to allow bonding.

Conclusion

In summary, we have developed an efficient and unprecedented metal-free and halogen-free method for the synthesis of 5hydroxuridine. Its access involves sequential oxidation mediated by hypervalent iodine(III) followed by mild deprotection reaction conditions. An in-depth NMR study showed that the presence of the 2',3'-O-isopolylidene moiety had a major influence on the sugar puckering, which was crucial for promoting oxidative cyclization. In addition to providing a rapid access to 5-OHU, a major RNA lesion, these results also extend the use of hypervalent iodine(II) in aqueous media.

Figure 1 .

 1 Figure 1. Selected region of the 2D NOESY spectrum of 3 (left) showing the correlation peaks used to characterize the stereochemistry of C6. These correlations were reported on the molecular structure (right).

  '-anhydro-8-oxo-2',3'-O-isopropylidene

Scheme 3 .

 3 Scheme 3. Synthesis of 5-OHU from a mixture of nucleosides 3 and 5

Table 1 .

 1 Coupling constants (J, Hz) for the ribose measured in CDCl3

Table 2 .

 2 Optimization of reaction conditions for the synthesis of cyclodihydrouridines 3 and 5

	Entry	R	Oxidant (eq)	Additives (eq)	Time (h)	3 (%)	5 (%)
	1	Ac	PhI(OAc)2 (2)	TEMPO (0.2)	3	47	4
	2	Ac	PhI(OAc)2 (2)	TEMPO (0.2) + NaHCO3 (2)	3	10	5
	3	Ac	PhI(OAc)2 (2)	-	3	70	8
	4	Ac	PhI(OAc)2 (2)	-	12	71	11
	5	Ac	PhI(OAc)2 (1.5)	-	3	67	11
	6	Ac	PhI(OAc)2 (5)	-	3	34	6
	7 a	Ac	PhI(OAc)2 (2)	-	3	37	5
	8 b	Ac	PhI(OAc)2 (2)	-	3	47	16
	9	Ac	PhI(OCF3)2 (2)	-	3	9	55
	10	Bz	PhI(OAc)2 (2)	-	3	18	2

a acetonitrile only; b water only

Primary Data

NO.

Funding Information

The GAVO consortium is gratefully acknowledged for financial support.

Supporting Information

YES (this text will be updated with links prior to publication)

Conflict of Interest

The authors declare no conflict of interest.


Epigenetic control is involved in molecular dialogue in plantmicrobe symbiosis

Keywords: ectomycorrhizal fungus, epigenetic control, host tree, hypomethylation, Laccaria bicolor, mycorrhizal symbiosis, poplar, transposable element ectomycorrhizal fungus, epigenetic control, host tree, hypomethylation, Laccaria bicolor, mycorrhizal symbiosis, poplar, transposable element

HAL is

2 Following the huge amount of insight in gene regulation at the transcriptional level by omic expression studies, revealing transcription factors and gene networks, a new step of whole genome regulation levels has added again more complexity to our understanding of gene expression control. Indeed, epigenetic processes affect also gene expression without changes in DNA sequences, as the chromatin structure is the molecular target of these mechanisms. In eukaryotes, the most common epigenetic modification is a DNA cytosine methylation that consists in the addition of a methyl group (CH 3 ) to a cytosine nucleotide. In plants, genomic DNA (gDNA) methylation is an extensively studied epigenetic mark that can occur in three sequence contexts, namely at the sites CG, CHG, and CHH (where H is A, T, C). Moreover, it is known that transposable elements (TEs) are presumed to be densely methylated and are thus inactive. Within last years, molecular mechanisms of DNA methylation have been extensively studied in model plants and recently reviewed [START_REF] To | Crosstalk among pathways to generate DNA methylome[END_REF][START_REF] Kakoulidou | DNA methylation remodeling in F1 hybrids[END_REF]. Methodological knowhow for such analyses has been remarkably progressed using whole genome bisulfite sequencing (WGBS) performed on gDNA for full methylome analyses [START_REF] Sow | RNAi suppression of DNA methylation affects the drought stress response and genome integrity in transgenic poplar[END_REF]. Moreover, machine learning for systematic detection of differentially methylated regions (DMRs) was developed and reported recently [START_REF] Hüther | MethylScore, a pipeline for accurate and context-aware identification of differentially methylated regions from population-scale plant whole-genome bisulfite sequencing data[END_REF]. However, few data are so far known regarding the link of epigenetics and plant-microbe interactions [START_REF] Satgé | Reprogramming of DNA methylation is critical for nodule development in Medicago truncatula[END_REF][START_REF] Ramos-Cruz | Epigenetics in plant organismic interactions[END_REF]. Regarding the formation of root nodules as new organs within the legume-Rhizobium symbiosis, [START_REF] Satgé | Reprogramming of DNA methylation is critical for nodule development in Medicago truncatula[END_REF] identified DNA methylation to play a determining role in Medicago truncatula nodules. A demethylase gene, called DEMETER (DME), is involved in the regulation of a high number of genes implicated in plant and bacterial cell differentiation required for nodule organogenesis in symbiotic interactions. Later, Ramos-Cruz et al. ( 2021) have stated that in addition to host genotypes and the environment, also the host epigenotype will determine the outcome of biotic interactions, potentially leading to plant defense, symbiosis or parasitism. In the case of arbuscular mycorrhizal (AM) symbiosis, first studies have shown that epigenetic modifications occur in seeds of Geranium sylvaticum or in roots and leaves of Geranium robertianum in symbiosis with the AM fungus Funneliformis mosseae [START_REF] Varga | Paternal arbuscular mycorrhizal fungal status affects DNA methylation in seeds[END_REF][START_REF] Varga | Arbuscular mycorrhizal fungi change host plant DNA methylation systemically[END_REF]. In this issue of New Here, the Authors gathered new insight concerning the effect of a modified epigenetic background of a host plant on its symbiotic relationship with a mycorrhizal fungal partner. [START_REF] Vigneaud | DNA hypomethylation of the host tree impairs interaction with mutualistic ectomycorrhizal fungus[END_REF] took the symbiotic model plant poplar (Populus tremula x Populus alba) and its ECM fungus Laccaria bicolor and used previously produced six independent epigenetic-modified plant lines obtained by overexpression or RNAi. These transgenic poplar lines were two overexpressinglines for the demethylases Demeter-Like DML (OX-dml1 and OX_dml2) and two RNAi lines for the chromatin remodeler DDM1 (Decreased in DNA Methylation 1). DDM1, conserved in animals, yeast, and plants, is known to facilitate methylation by providing access for DNA methyltransferase. A strong decrease of DDM1 expression levels in the model plant Arabidopsis thaliana resulted in a significant reduction of DNA methylation [START_REF] Lippman | Role of transposable elements in heterochromatin and epigenetic control[END_REF]. All these four poplar lines are hypomethylated lines. The two other used poplar mutants were RNAi-mediated hypermethylated DML lines [START_REF] Sow | RNAi suppression of DNA methylation affects the drought stress response and genome integrity in transgenic poplar[END_REF]. The level of DNA methylation of the host tree being manipulated in these transgenic plant lines, the authors observed differential gene and TE methylation together with differential gene expression, linked to reduced mycorrhization rates only for the hypomethylated lines. Comparison of wild-type poplar with poplar lines with modified DNA methylation levels showed not only differentially methylated plant genes (86) and TEs (288) and, remarkably, differentially methylated fungal genes (120) and TEs (1441), but also that such manipulation will affect in turn mycorrhizal interactions. In fact, the fungus methylome will become altered when the colonizing fungus interacts with poplar hosts with hypomethylated gDNA, suggesting that disturbing the control of poplar methylation will induce also the remodeling of the fungal genomic methylome. Moreover, hypomethylated poplar results in differential gene (205) expression. From these results, the authors suggested poplar gene candidates potentially involved in physiological processes that could be linked to mycorrhizal interaction, as root initiation, immune responses, hormonal pathways (ethylene and jasmonate), or terpenoid metabolism. Such candidates will need further analyses to rule out their specific roles and regulation. The study by [START_REF] Vigneaud | DNA hypomethylation of the host tree impairs interaction with mutualistic ectomycorrhizal fungus[END_REF] concludes that a number of host genes under epigenetic control is involved in the mutualistic interaction highlighting a new level of gene regulation to the rather complex interaction between host plants and microbial symbionts. In the context of ECM symbiosis playing a major role for nutrition, health and stress tolerance of host trees, these new findings increase our knowledge about the tightly regulated crosstalk between symbiotic partners.

Interestingly, the impact of DNA methylation for plant genomic immunity has been demonstrated [START_REF] Kim | DNA methylation as a system of plant genomic immunity[END_REF][START_REF] Parker | Epigenetics: a catalyst of plant immunity against pathogens[END_REF][START_REF] Huang | Coordinated epigenetic regulation in plants: a potent managerial tool to conquer biotic stress[END_REF]. Plants have evolved sophisticated and tightly regulated methylation and demethylation pathways to protect their genes against damage by transposon invasions [START_REF] Kim | DNA methylation as a system of plant genomic immunity[END_REF]. In response to biotic stress, epigenetic regulation, including DNA methylation, histone and chromatin structure modifications, has been reported to be involved in plant immunity against interacting pathogens [START_REF] Huang | Coordinated epigenetic regulation in plants: a potent managerial tool to conquer biotic stress[END_REF]. As plants as sessile organisms without specialized cells for immune responses need to tolerate and resist a number of abiotic and biotic stress conditions, among them pathogen attacks, these epigenetic regulation tools are of high importance for survival and adaptation. Moreover, epigenetic modifications can be heritable and thus contribute to the defense priming and transgenerational memory. In fact, biotic stress will alter DNA methylation and modification of the epigenetic regulation would affect resistance to plant diseases [START_REF] Parker | Epigenetics: a catalyst of plant immunity against pathogens[END_REF]. Research findings in this field have clearly demonstrated a link between plant immunity and TE hypomethylation (Wilkinson et al., 2019).

As beneficial plant-microbe interactions occur at a border of plant defense and establishment of symbiotic relations, it could be expected that a tight epigenetic control between host and symbiont would also be involved and needed in this process. Moreover, comparison between beneficial symbiotic interactions and pathogenic attacks might give further clues on the specific impact of epigenetic control between "friend and foe". In perspective, a comparison of epigenetically regulated plant genes as reported in the context of drought stress [START_REF] Sow | RNAi suppression of DNA methylation affects the drought stress response and genome integrity in transgenic poplar[END_REF] without and with the fungal ectomycorrhizal partner in the wild-type situation would be interesting to advance the knowledge on the reciprocal crosstalk between both symbiotic partners. In addition, an interesting open question is 5 whether epigenetic regulatory mechanisms might be involved in induced systemic resistance (ISR) by beneficial microbes [START_REF] Pieterse | Induced systemic resistance by beneficial microbes[END_REF].

Altogether, [START_REF] Vigneaud | DNA hypomethylation of the host tree impairs interaction with mutualistic ectomycorrhizal fungus[END_REF] convincingly and nicely demonstrated with their contribution the impact of the host plant epigenetic machinery on the interaction with a beneficial ectomycorrhizal fungus. A series of further challenging questions concerning the crosstalk between host plant and beneficial microbes, the specificity of beneficial interactions compared to pathogen attacks, and the involvement of epigenetic regulation in reported ISR might be raised in future research. Finally, more generally, these questions will become important to master plant-microbe interactions, priming of defense and/or better use of beneficial associations.

1

  Figures: 0

   analyzed gene regulation at the epigenetic level interfering with ectomycorrhizal (ECM) interaction.


Malak Salame

Crystel Bonnet

Ely Cheikh

Mohamed Moctar

Selma Mohamed Brahim

Abdallahi Dedy

Abdel Ledour

Fatimetou Vetah

Cheikh Tijani Veten

Christine Hamed

Ahmed Petit

Houmeida

Identification a novel pathogenic LRTOMT mutation in Mauritanian families with nonsyndromic deafness

Keywords: SNHI, Deafness, LRTOMT, Mutation, Mauritania

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Deafness is the most common sensorineural disability in developing countries as it affects approximately 5% of the population [START_REF] Wilson | Global hearing health care: New findings and perspectives[END_REF][START_REF] Sheffield | The Epidemiology of Deafness[END_REF]. In children, 1 to 6 per 1,000 newborns are suffering from a form of hearing impairment [START_REF] Kimberly | Hearing Loss in Children[END_REF]. Despite the lack of reliable statistics on the affected populations in developing countries, the proportion of hearing loss was estimated to be greater than in states with higher income [START_REF] Prasansuk | Incidence/prevalence of sensorineural hearing impairment in Thailand and Southeast Asia[END_REF][START_REF] Elahi | Paediatric hearing loss in rural Pakistan[END_REF]. In a random WHO survey, the prevalence of deafness was indeed found to be closely affected by widespread issues such as deprived life environment, lack of proper care services but also genetic background including the effect of consanguinity [START_REF] Smith | WHO activities for prevention of deafness and hearing impairment in children[END_REF]7]. SNHI were related to about 134 identified causal genes with most cases passed down through variations in GJB2 gene with autosomal inheritance (Hereditary Hearing loss Home page: Last update August 2021). Several variations of LRTOMT gene, mostly in exon 7, have also been associated with deafness in different populations including cohorts from North African populations making of this gene the second most frequent, after GJB2, involved in congenital severe to profound hearing impairment [8,[START_REF] Mosratima | Deep analysis of the LRTOMT c.242G>A variant in non-syndromic hearing loss North African patients and the Berber population: Implications for genetic diagnosis and genealogical studies[END_REF][START_REF] Ichinose | Novel Mutations in LRTOMT Associated With Moderate Progressive Hearing Loss in Autosomal Recessive Inheritance[END_REF][START_REF] Charif | The c.242G>A mutation in LRTOMT gene is responsible for a high prevalence of deafness in the Moroccan population[END_REF][START_REF] Ahmed | Mutations of LRTOMT a fusion gene with alternative reading frames cause non syndromic deafness in humans[END_REF][START_REF] Riahi | Whole exome sequencing identifies new causative mutations in Tunisian families with non-syndromic deafness[END_REF]. In a previous work using Sanger sequencing of GJB2 gene, we have detected two pathogenic variants c.35delG (p.Gly12Valfs*2) and c.94C>T (p.Arg32Cys) in only 5 out of the 53 families with sporadic or hereditary deafness investigated [8]. This finding left thus about 95% of the families with other putative genetic cause to be assigned. The aim of this study was to identify the molecular etiology of nonsyndromic hearing loss not linked to GJB2 mutations in a large cohort of Mauritanian patients using a targeted mutation screening of exon 7 in LRTOMT gene.

Materials and methods

Patients

The cohort was composed of 128 unrelated Mauritanian children (1 child per family) respectively from 83 families with multiple affected siblings and 45 families with no history of deafness (sporadic cases), recruited from two deaf schools in Nouakchott, Mauritania.

Forty-one children of the cohort were identified as negative for GJB2 variants in a previous study we carried out on non-syndromic hearing impairment associated GJB2 variants in Mauritania [8] and 87 children were newly recruited. All patients were diagnosed with severe-to-profound non-syndromic congenital hearing impairment following medical examination. This study was approved by the ethics committee of the University of Nouakchott, Mauritania as stated by the ethics clearance letter N o 002/2020/CE/UNA. Written informed consent was obtained from all parents or legal guardians of the children. A questionnaire was used to collect demographic data and medical history of the subjects.

Molecular Analysis

Blood samples were collected in EDTA tubes and used to extract genomic DNA by Qiagen DNA Blood Minikit procedure (QIAGEN Genomic DNA Handbook.2015). PCR reactions (20 µl) contained 1 µl (20-30 ng) of genomic DNA, 1 µl (10µM) of each specific primer (5'-AGGATAATAATTGCTACTGGCAAAA-3'and5' ATCCCAAATATTCCTTCACTGTCTT-3'), 10µl of AmpliTaqGold 360 Master Mix (Life Technologies) and 7 µl of distilled water. The PCR program presented a DNA denaturation step at 95°C for 10 min followed by 35 amplification cycles (denaturation at 94°C for 30s, annealing at 53°C for 35s and extension at 72°C for 40s) tailed by a 5 min final extension at 72°C. PCR products were sequenced on capillary ABI3730 Genetic Analyzer (Applied Biosystems, California, USA). Data obtained were then matched with reference sequences of the LRTOMT gene (NM_001145309.4) using Seqscape3.0 software program package (Gene Codes, MI, and USA). Amino acid change was considered as potentially deleterious if predicted by (PolyPhen2, http://genetics.bwh.harvard.edu/pph2), Sorting Intolerant from Tolerant (SIFT, http://sift.jcvi.org/) and Mutation Taster (http://www.mutationtaster.org/). In families with multiple cases of deafness, DNA of one proband was sequenced. Sanger sequencing was then extended to available parents and all affected siblings carrying the mutation in homozygous state.

LRTOMT protein 3-D structure modeling

The predicted tertiary structure of LRTOMT protein (AF-Q8WZ04-F1) was obtained using AlphaFold protein structure database (https://alphafold.ebi.ac.uk/). Similarity between wildtype and mutated 3D structures of LRTOMT protein was explored by using UCSF Chimera version 1.15.

Results

Screening of exon7 in LRTOMT

Analysis of sequencing data from exon 7 of LRTOMT gene of the 128 selected patients showed only one biallelic missense variant (c.179T>C;p.(Leu60Pro) in four (3.1%) unrelated patients. The genomic coordinates of the newly identified variant were Chr11(GRCh37):g.71817077T>C; NM_001145309.4:c.179T>C. Sequence alignment revealed that this variant was not reported before. Notably, the arginine residue was highly conserved through evolution among multiple-species (from Human to opossum considering 6 species) (Fig. 1). The variant was predicted to have deleterious effect on gene function by SIFT (score: 0.01) and a disease-causing effect by Mutation Taster (prob: 1). Exploration of effect of p.Leu60Pro on the predicted three-dimensional (3D) structure of the encoded protein showed a disruption from an organized α helix in the wild protein, to a random conformation in the mutant form (Fig. 1). The healthy parents of the four patients in families (F7-DFN, F28-DFN, F31-DFN, F81-DFN) carried the mutation in heterozygous state and all affected children were homozygous for the variant. Three of the four patients had other affected relatives and one presented no family history of HI. The four affected families originated from the Trarza region in southern Mauritania (Fig. 2 and3).

Cochlear implant outcome

Based on clinical assessment and audiological data, the proband of family F28-DFN, a 7year-old boy, diagnosed with profound bilateral congenital deafness, received a cochlear implant fitted in the right ear. Now, aged 9 years, signs of auditory recovery are largely perceptible (Fig. 3) on the basis of his score on APCEI scale used to assess his communication skill progress.

Discussion

In this study, we identified a novel biallelic missense pathogenic variant in exon 7 of LRTOMT gene in 3 Mauritanian children with family history of congenital deafness and one child with no other relative being affected. The variant was found in both girls and boys and only carriers at homozygous state showed the phenotype which supported an autosomal recessive segregation of the variant. Because of their higher number of inter-residue contacts, helices are known to tolerate sequence amino-acid changes without affecting their secondary structure [START_REF] Abrusán | Marsh Alpha Helices Are More Robust to Mutations than Beta Strands[END_REF]. As a result of this relative robustness, compared to strands and coils conformations, any variation in the helix organization was implied to be enough significant to cause alteration of the protein structure and therefore likely to affect its function. For instance, the wild residue Arg81 formed a salt bridge and hydrogen bonds between the helix and following loop, central for the protein stability. Its substitution with glutamine c.242G>A in exon 7: (p.Arg81Gln) was found to disrupt the alpha helix resulting in local destabilization of the protein tertiary structure [START_REF] Ahmed | Mutations of LRTOMT a fusion gene with alternative reading frames cause non syndromic deafness in humans[END_REF]. Although involving a radical structurally distinct, our variant (p.Leu60Pro), located in the same exon 7, appeared to be inducing a comparable disorganizing effect on the protein alpha helix which likely resulted in the hearing impairment observed in the four patients of our cohort and also predicted for the (p.Arg81Gln) mutation [START_REF] Charif | The c.242G>A mutation in LRTOMT gene is responsible for a high prevalence of deafness in the Moroccan population[END_REF]. Since the first reported case of direct involvement of proline in human diseases [START_REF] Schafer | Familial hyperprolinemia, cerebral dysfunction and renal anomalies occuring in a family with hereditary nephropathy and deafness[END_REF], mutant proteins for this amino acid have been implicated in the etiology of number of genetic disorders including non-syndromic recessive hearing impairment [START_REF] Batissoco | A novel missense mutation p.L76P in the GJB2 gene causing non-syndromic recessive deafness in a Brazilian family Brazilian[END_REF]. For instance, different studies showed that p.Ser50Pro and p.Thr56Pro variants, both consisting of replacement of wild amino acid by proline residue in connexin 50 (Cx50, GJA8), affected the protein contribution to the lens normal survival through angular restriction of α helix or β sheet secondary structure [START_REF] Adam | The cataract causing Cx50-S50P mutant inhibits Cx43 and intercellular communication in the lens epithelium[END_REF][START_REF] Hadrami | A novel missense mutation of GJA8 causes congenital cataract in a large Mauritanian family[END_REF]. Proline is indeed a distinctive amino acid with a nitrogen atom covalently locked within a ring inducing a unique constrained φ angle of approximately -65° [START_REF] Alexander | Proline: The Distribution, Frequency, Positioning, and Common Functional Roles of Proline and Polyproline Sequences in the Human Proteome[END_REF]. This irregular geometry disrupts protein secondary structure by steric hindrance i.e. in inhibiting the backbone to conform into a normal alpha-helix turn conformation. In addition, the lack of hydrogen on proline's nitrogen (proline cannot donate protons) prevents it from participating in hydrogen bonding essential for the helix structure. The apparent deleterious consequence we observed (SIFT score (0.01) and disease-causing effect (prob: 1) by Mutation Taster) for this novel mutation was therefore likely due to a change of the protein secondary structure induced by the substitution to proline residue consisting in an alpha helix breaking. The wild type leucine residue was also located in an evolutionary highly conserved region [START_REF] Ahmed | Mutations of LRTOMT a fusion gene with alternative reading frames cause non syndromic deafness in humans[END_REF] suggesting an important functional position of the residue (Fig 1B). Modeling, using WHATIF server (http://swift.cmbi.ru.nl) on rat COMT crystal structure as template for catechol-Omethyltransferase domain of the human LRTOMT2 (NM_001145309.4) showed that common LRTOMT variants associated with non-syndromic recessive hearing impairment led to hair cell degeneration process [START_REF] Mosratima | Deep analysis of the LRTOMT c.242G>A variant in non-syndromic hearing loss North African patients and the Berber population: Implications for genetic diagnosis and genealogical studies[END_REF][START_REF] Charif | The c.242G>A mutation in LRTOMT gene is responsible for a high prevalence of deafness in the Moroccan population[END_REF][START_REF] Gibriel | Analysis of p.Gly12Valfs*2, p.Trp24* and p.Trp77Arg mutations in GJB2and p.Arg81Gln variant in LRTOMT among non-syndromic hearing loss Egyptian patients: implications for genetic diagnosis[END_REF]. Close to p.Leu60Pro variant, various changes in the same exon 7 were also associated with cases of moderate to profound congenital deafness making of this segment a mutation hot spot in the LRTOMT gene [START_REF] Ichinose | Novel Mutations in LRTOMT Associated With Moderate Progressive Hearing Loss in Autosomal Recessive Inheritance[END_REF]. Besides the above homozygous c.242G>A (p.Arg81Gln) variant found in Tunisian, Libyan and Egyptian families [START_REF] Mosratima | Deep analysis of the LRTOMT c.242G>A variant in non-syndromic hearing loss North African patients and the Berber population: Implications for genetic diagnosis and genealogical studies[END_REF][START_REF] Gibriel | Analysis of p.Gly12Valfs*2, p.Trp24* and p.Trp77Arg mutations in GJB2and p.Arg81Gln variant in LRTOMT among non-syndromic hearing loss Egyptian patients: implications for genetic diagnosis[END_REF] we may also cite missense variants c.122G>A (p.Arg41Gln), c.193T>C (p.Trp65Arg) in Morocco and Tunisia [START_REF] Charif | The c.242G>A mutation in LRTOMT gene is responsible for a high prevalence of deafness in the Moroccan population[END_REF][START_REF] Ahmed | Mutations of LRTOMT a fusion gene with alternative reading frames cause non syndromic deafness in humans[END_REF] and a nonsense mutation c.208C>T (p.Arg70*), in Tunisia [START_REF] Charif | The c.242G>A mutation in LRTOMT gene is responsible for a high prevalence of deafness in the Moroccan population[END_REF][START_REF] Ahmed | Mutations of LRTOMT a fusion gene with alternative reading frames cause non syndromic deafness in humans[END_REF][START_REF] Riahi | Whole exome sequencing identifies new causative mutations in Tunisian families with non-syndromic deafness[END_REF] (Table 1). In this ethnic context, we noticed that all four affected families carrying the LRTOMT mutation (c.179T>C;p.(Leu60Pro) identified in this report belonged to the Moors (Maures), a predominant racial group of the Mauritanian population. This group, of Berber-Arab origin, ethnically and culturally self identifies with the neighboring North Africa populations [START_REF] Vernet | La Mauritanie, des Origines au Début de l'Histoire[END_REF].

In a previous study, we detected GJB2 (MIM 121 011) variants in only 5 out of the 53 (9.4%) investigated families with congenital hearing impairment [8]. The present work showed that LRTOMT was also involved in hearing impairment in our population (3.1%; 4/128). However, a limitation of this study was that we have only screened exon 7 of LRTOMT. This screening only in exon 7 was mostly guided by very scarce findings of deafness associated variants in other exons reported by various studies in the region [START_REF] Mosratima | Deep analysis of the LRTOMT c.242G>A variant in non-syndromic hearing loss North African patients and the Berber population: Implications for genetic diagnosis and genealogical studies[END_REF][START_REF] Ichinose | Novel Mutations in LRTOMT Associated With Moderate Progressive Hearing Loss in Autosomal Recessive Inheritance[END_REF][START_REF] Charif | The c.242G>A mutation in LRTOMT gene is responsible for a high prevalence of deafness in the Moroccan population[END_REF][START_REF] Ahmed | Mutations of LRTOMT a fusion gene with alternative reading frames cause non syndromic deafness in humans[END_REF][START_REF] Riahi | Whole exome sequencing identifies new causative mutations in Tunisian families with non-syndromic deafness[END_REF].

Indeed, a total of 14 pathogenic variants in exon 7 out of 18 causal variants have presently been reported. Search in these sequences along multi-gene hear panel is underway to locate other deafness related variants in our families.

Conclusions

This study identified a novel biallelic predicted pathogenic variant (c. Screening for other parts of this gene or using multigene panel in a larger population study is also recommended.

  179T>C;p.(Leu60Pro) in exon 7 of LRTOMT gene likely causing a localized alpha helix breaking in the LRTOMT protein. Being the single exon 7 variant identified in our patients underlined the significance of a genetic screening for this variant in the context of an accurate molecular diagnosis of hearing impairment or in assessing the transmission of the disease in our population and, broadly, among Northern African populations.Further hearing tests are still required to clearly establish the effect of early cochlear implant fitting on audition rehabilitation in patients carrying mutations in the LRTOMT gene.

Figure 1 :

 1 Figure 1: Sequencing data, alignment and prediction of protein secondary structure of the LRTOMT protein. (A) Electropherogram plots showing the reference sequence and the mutation in heterozygous and homozygous states. (B) Sequence alignment presenting conservation of the region of LRTOMT protein containing (Leu60) among various species. (C) LRTOMT protein tertiary structure modeling performed with Chimera, based on the AlphaFold prediction of the transmembrane O-methyltransferase (AF-Q8WZ0F1). (D) Localized structural change in the wild type alpha helix (Phe57, Leu59, Leu60, Val61, Thr63, Val61, Val64 and Leu65) (left) showing the breaking of the alpha helix induced by substitution of Leu60 by Prolin (right).

Figure 2 :

 2 Figure 2: Pedigrees of affected children with family history. Segregation of the mutation pattern was compatible with autosomal recessive inheritance. All cases originated from southern Mauritania (map top right).

Figure 3 :

 3 Figure 3: Pure-tone audiogram of a patient from family F28-DFN who received cochlear implant. Proband, indicated by an arrow, had profound hearing loss in both ears.

Table 1 :

 1 Overview of known pathogenic variants in the LRTOMT Gene and associated phenotypes. Report from Deafness variation Database (https://deafnessvariationdatabase.org/gene/LRTOMT)

	Exon	Variants	Amino acid	Phenotype	Population	Reference
	number		change
	7	c.96_109del p.Thr33HisfsTer57	Non syndromic-HL	USA	[22]

Acknowledgments

We thank all the families for their participation to the study. This work was supported by Fondation pour l'Audition (FPA-IDA05) and the ANRSI-Mauritania.

Funding

The authors declare that no funds or grants were received during the preparation of this manuscript.

Availability of data and materials

The datasets on variants generated during the current study are available from the corresponding author on reasonable request.

Declarations

Ethics approval and consent to participate

Research in this study has been performed in accordance with the Declaration of Helsinki and was approved by an appropriate ethics committee: the Ethics Committee of the University of Nouakchott, Mauritania (ref. n o 002/2020/CE/UNA). We confirm that all methods were performed in accordance with the relevant guidelines and regulations.

Consent for publication

Informed consent was obtained from all study participants. The informed consent of all patients was also obtained for the data publication.

Competing interests

The authors declare that they have no competing interests.


Z Lisicic

A-C Tadic

M Fabre

A Fagan

Herrel

M Taverne
email: maxime.taverne@mnhn.fr

P J Watson²

H Dutel

R Boistel

D Lisicic

Z Tadic

A-C Fabre

M J Fagan

Form function relationships underlie rapid dietary changes in a lizard

Keywords: Adaptation, diet, bite force, lizard, finite element models, multibody dynamics

Introduction

Macroevolutionary changes such as changes in habitat use or diet are often associated with convergent, adaptive changes in morphology. Previous studies have suggested that for macroevolutionary changes to occur, directional selection driven by consistent changes in the environment is needed. Consequently, Carroll and colleagues (2007) suggested that "macroevolution may thus be nothing more than an aggregate of many small events". Despite the initial assertion that macroevolutionary patterns cannot be predicted from processes at the population level [START_REF] Eldredge | Phylogenetic Patterns and the Evolutionary Process[END_REF], subsequent authors have demonstrated that patterns of variation among taxa can be predicted using population genetics theory [START_REF] Arnold | The adaptive landscape as a conceptual bridge between micro-and macroevolution[END_REF]. Moreover, in some cases, variation in morphology can be rapid and may drive the subsequent evolutionary trajectory of a population [START_REF] Herrel | Rapid large-scale evolutionary divergence in morphology and performance associated with exploitation of a different dietary resource[END_REF]. Indeed, variation in fitness-relevant traits has the potential to rapidly drive a population across a valley of low fitness to a new adaptive peak [START_REF] Arnold | Morphology, performance, and fitness[END_REF][START_REF] Schluter | The ecology of adaptive radiation[END_REF][START_REF] Arnold | The adaptive landscape as a conceptual bridge between micro-and macroevolution[END_REF][START_REF] Svensson | The adaptive landscape in evolutionary biology[END_REF], thus potentially driving rapid and seemingly punctuated changes in morphology [START_REF] Gingerich | Rates of evolution: effects of time and temporal scaling[END_REF][START_REF] Eldredge | Time Frames: The evolution of punctuated equilibria[END_REF]. Yet, for small-scale population-level variation to facilitate or drive rapid shifts in ecology and potentially allow organisms to reach these new adaptive peaks, this variation in morphology needs to have a significant impact on function. As the link between morphology and function is often non-linear [START_REF] Koehl | When does morphology matter[END_REF], even small differences in form may potentially give rise to significant differences in function. For example, since the force output of a musculoskeletal system scales to the second power of linear dimensions, small changes in head dimensions can generate significant differences in bite force [START_REF] Herrel | Ontogenetic scaling of bite force in lizards and turtles[END_REF], which might in turn allow species to access novel resources [START_REF] Aguirre | The implications of food hardness for diet in bats[END_REF][START_REF] Herrel | Omnivory in lacertid lizards: adaptive evolution or constraint?[END_REF]Herrel et al., , 2006[START_REF] Herrel | Rapid large-scale evolutionary divergence in morphology and performance associated with exploitation of a different dietary resource[END_REF].

Herbivory is an attractive ecological strategy that, at least in mammals, has resulted in fast species diversification [START_REF] Price | Tempo of trophic evolution and its impact on mammalian diversification[END_REF] with over 38% of all mammals being considered herbivores [START_REF] Atwood | Herbivores at the highest risk of extinction among mammals, birds, and reptiles[END_REF]. In contrast to mammals, squamate herbivores are rare with merely 2% of all species being considered herbivores [START_REF] Cooper | Distribution, extent, and evolution of plant consumption by lizards[END_REF]. Among lizards this value may be as high as 5% [START_REF] Meiri | Traits of lizards of the world: Variation around a successful evolutionary design[END_REF], but herbivory remains an uncommon dietary niche. Among the reasons that have been suggested to have prevented lizards from becoming herbivorous are their small body size [START_REF] Szarski | Some remarks on herbivorous lizards[END_REF][START_REF] Ostrom | Further comments on herbivorous lizards[END_REF][START_REF] Pough | Lizard energetics and diet[END_REF], unspecialized dentition [START_REF] Hotton | A survey of adaptive relationships of dentition to diet in the North American Iguanidae[END_REF], and lack of complex food processing [START_REF] King | Reptiles and herbivory[END_REF]. However, it has since been demonstrated that even small lizards can be herbivorous [START_REF] Espinoza | Recurrent evolution of herbivory in small, cold-climate lizards: breaking the ecophysiological rules of reptilian herbivory[END_REF] and that complex cusped teeth have evolved associated with an herbivorous diet [START_REF] Lafuma | Multiple evolutionary origins and losses of tooth complexity in squamates[END_REF]. Moreover, herbivorous lizards typically have higher bite forces than insectivorous species [START_REF] Herrel | Omnivory in lacertid lizards: adaptive evolution or constraint?[END_REF][START_REF] Herrel | Lizard Ecology: The evolutionary consequences of foraging mode[END_REF] which may allow them to crop leaves from a larger plant. The herbivorous niche in lizards thus appears to be an adaptive peak that is rather difficult to attain and needs to be accompanied by a suite of distinct anatomical and functional specializations.

A notable exception to the general idea that herbivory is difficult to achieve in lizards is the previously documented rapid evolution of a largely herbivorous diet in a population of Italian wall lizards (P. sicula) roughly three decades after its introduction onto a small islet in the Adriatic [START_REF] Nevo | Competitive exclusion between insular Lacerta species (Sauria, Lacertidae). Notes on experimental introductions[END_REF][START_REF] Herrel | Rapid large-scale evolutionary divergence in morphology and performance associated with exploitation of a different dietary resource[END_REF]. Indeed, lizards that were introduced onto the islet of Pod Mrčaru from the neighboring islet of Pod Kopište switched to a diet composed of up to 60% of plants in summer and show physiological adaptations, as well as differences in their microbiome, allowing them to more efficiently extract energy from a plant-based diet [START_REF] Herrel | Rapid large-scale evolutionary divergence in morphology and performance associated with exploitation of a different dietary resource[END_REF][START_REF] Wehrle | Rapid dietary shift in Podarcis sicula resulted in localized changes in gut function[END_REF][START_REF] Lemieux-Labonté | Proximate drivers of population-level lizard gut microbial diversity: impacts of diet, insularity, and local environment[END_REF]. Moreover, in the 36 years since the introduction of these lizards on Pod Mrčaru, they developed larger heads, muscles, and bite forces and changed the shape of their cranium and mandible [START_REF] Herrel | Rapid large-scale evolutionary divergence in morphology and performance associated with exploitation of a different dietary resource[END_REF][START_REF] Taverne | From micro to macroevolution: drivers of shape variation in an island radiation of Podarcis lizards in the Adriatic[END_REF]. As the two islets are similar in size and show similar vegetation types and high lizard densities [START_REF] Herrel | Rapid large-scale evolutionary divergence in morphology and performance associated with exploitation of a different dietary resource[END_REF] the ecological drivers of the observed changes in diet remain unclear. Moreover, despite the documented changes in head shape and muscle architecture it remains unclear how these lizards were able to gain a great enough functional advantage to allow them to occupy a new adaptive peak in the fitness landscape (i.e., a plant-based diet). Here, we use dissections, geometric morphometrics, and mechanical engineering tools, including multibody dynamics analysis (MDA) and finite element modeling (FEM), to better understand the functional advantages provided by the rather subtle changes in skull shape and muscle architecture previously documented in these two populations [START_REF] Taverne | From micro to macroevolution: drivers of shape variation in an island radiation of Podarcis lizards in the Adriatic[END_REF]. We demonstrate how an intricate co-evolution of skull and mandible shape with muscle architecture can allow for a more efficient mechanical transfer of forces from the muscles to the jaws, and a more resistant skull configuration when biting. These results provide insights into how subtle phenotypic variation may give rise to fitness-relevant changes in function allowing the rapid transition towards new adaptive peaks and the occupation of novel trophic niches in lizards.

Material and methods

Quantification and comparison of head shape and muscular anatomy

Thirteen male specimens from Pod Kopište and 14 male specimens from Pod Mrčaru were captured by hand or by noose at the end of the summer of 2013. Their body size (SVL: snout-vent length) and head size was measured with a Mitotoyo digital caliper (± 0.01 mm) and their bite force was measured using a custom-designed bite force set-up (Herrel et al., 1999b). They were sacrificed by an intramuscular injection of pentobarbital under a permit of the Croatian Ministry of the Environment. Specimens were preserved in a 10% aqueous formaldehyde solution for 48h, rinsed and stored in a 70% aqueous ethanol solution.

The heads of these 27 specimens were scanned at the University of Poitiers on an Easytom micro-CT (at a voxel size of 24.90 μm with the following parameters: X-ray voltage, 90 kV; X-ray intensity, 70 μA; exposure time, 2000 ms; number of projections, 2500). Scans were imported into Avizo 9.0 (Thermo Fisher Scientific) to segment the mandible and skull which were exported as PLY files. A set of landmarks and semi-landmarks on curves were used to quantify the skull and mandible shape (see [START_REF] Taverne | From micro to macroevolution: drivers of shape variation in an island radiation of Podarcis lizards in the Adriatic[END_REF]Figure S1, Table S1). Semi-landmarks on curves were slid while minimizing the bending energy and all landmarks were aligned by a Procrustes superimposition with the function "gpagen" ("geomorph" package). The function "prcomp" from the stats package was used to run a Principal Component Analysis (PCA) on the Procrustes coordinates. A MANOVA on the principal components (PC) cumulatively explaining at least 85% of the variance was run to test for shape differences between populations. The theoretical shapes corresponding to the extremes of the PCs that distinguish populations were visualized using the functions "tps3d" and "shade3d" ("Morpho" package).

Five additional male specimens per population were included to quantify variation in jaw muscle architecture (resulting in a total 18 specimens from Pod Kopište and 19 specimens from Pod Mrčaru). Each muscle bundle was extracted by dissecting the left side of the head. Muscle bundles were blotted dry and weighed with a digital balance (Mettler AE100; ± 0.1 mg). The connective tissue surrounding the muscles was digested by submerging the muscles in a 30% aqueous nitric acid solution for 24h. Next, the nitric acid was removed, and a 50% glycerol solution was added to arrest the muscle digestion. Muscle fibers were drawn under a binocular scope (Leica) with camera lucida (see [START_REF] Taverne | From micro to macroevolution: drivers of shape variation in an island radiation of Podarcis lizards in the Adriatic[END_REF] and measured using Image J (Rasband 1997). Muscle volume was calculated as the ratio between muscle mass and muscle density (1.06 g.cm -3 ; see [START_REF] Mendez | Density of fat and bone mineral of mammalian muscle[END_REF]. The physiological cross-sectional area (PCSA) of each muscle bundle was subsequently calculated by dividing muscle volume by the mean fiber length and was subsequently corrected for pennation angle. The muscles were grouped into four functional groups: the external adductors, the pseudotemporalis group, the adductor posterior, and the pterygoid group (Table S2). The jaw depressors and the constrictor dorsalis muscles were not considered since they are not involved in jaw closing.

Differences in the muscular architecture (muscle volume, PCSA, fiber length) between populations were tested by means of a multivariate analysis of covariance (MANCOVA) with the island as factor and SVL as co-variable, using the function "mancova" from the package "jmv". Subsequent univariate ANOVAs were used to test which muscle groups were responsible for the observed differences.

Multibody dynamics model design

The skull and the mandible of two males of similar size, one from the island of Pod Kopište (PK), one from Pod Mrčaru (PM), were segmented using Avizo 9.0 (Thermo Fisher Scientific). The 3D models of the cranium and jaw of each specimen were used to build two MDA models. The location of the origin and insertion of each jaw muscle was assessed based on dissections and each muscle was virtually divided into several distinct muscle bundles according to the size of the muscle. The number of bundles used was based on the muscle volume and area of insertion (Table 1). The coordinates of the sites of origin and insertion of each muscle bundle were determined using Avizo 9.0 (Thermo Fisher Scientific). The mandible was converted into a parasolid format allowing to compute the inertial properties based on a bone density of 1.5 g.cm -3 [START_REF] Sellers | Using sensitivity analysis to validate the predictions of a biomechanical model of bite forces[END_REF]. The skull surface, the mandible parasolid, and the muscle bundle coordinates were imported into MSC ADAMS multibody dynamic simulation software (Santa Ana, CA, USA). The muscle bundles were modelled as contractile springs. Where necessary, muscle bundles were wrapped around the bone to increase the accuracy of modelling [START_REF] Gröning | The importance of accurate muscle modelling for biomechanical analyses: a case study with a lizard skull[END_REF]Fig. 1). Muscle PCSAs were corrected by both the pennation angle and the typical 10% loss in muscle volume due to tissue preservation [START_REF] Kikuchi | Differences in muscle dimensional parameters between non-formalinfixed (freeze-thawed) and formalin-fixed specimen in Gorilla (Gorilla gorilla)[END_REF]. The maximum force of each muscle (Table 1) was calculated by multiplying the PCSA by an intrinsic muscle stress of 40 N/cm² [START_REF] Gröning | The importance of accurate muscle modelling for biomechanical analyses: a case study with a lizard skull[END_REF]. The MDA models simulated opening and closing of the jaw assuming maximal activation of the bundles when the system was at equilibrium. A food particle was created and aligned perpendicularly to the toothrow of the upper jaw to generate a reaction force due to jaw closing (simulated bite force). To ensure consistency in the output bite force generated by the MDA models, the location of the contact between the food particle and the teeth was prescribed to match the location of the contact between the teeth and the plates of the force transducer used to measure in vivo bite force in the field [START_REF] Taverne | Proximate and ultimate drivers of variation in bite force in insular lizards Podarcis melisellensis and Podarcis siculus[END_REF]. To run further simulations, the location of the item was then standardized at the middle of the maxillary toothrow, as observations revealed that lizards typically crush prey at that location [START_REF] Taverne | Relationships between dietary breadth and flexibility in jaw movement: a case study of two recently diverged insular populations of Podarcis lizards[END_REF]. The stiffness of the food particle was intentionally set beyond the hardness of prey typically consumed by the lizards to ensure that the gape angle did not change during a biting simulation and to obtain the maximal bite force.

Finite element model design

Two finite-element meshes consisting of about 1.5 million tetrahedral elements (Pk: 1,527,268 elements; PM: 1,321,278 elements) were generated in Avizo and imported into ANSYS (ANSYS, Inc., Canonsburg, PA, USA) for finite element analysis. Bone was assumed to have isotropic homogenous material properties with a Young's modulus of 17 GPa and a Poisson's ratio of 0.3 (Cowin, 2001;[START_REF] Rayfield | Using finite-element analysis to investigate suture morphology: a case study using large carnivorous dinosaurs[END_REF]. The meshes were constrained at the ventral base of the quadrates (the right side in all three directions, the left in A-P and D-V only) and at two symmetrical bite points in the D-V direction (Fig. 1). The coordinates of these locations were directly exported from ADAMS, with the three-dimensionally resolved forces exported from the MDA solutions then directly applied onto the skull model. The quadrato-jugal and the epipterygoid-parietal ligaments were modelled as tensiononly links with a stiffness of 250 N/mm 2 and a cross-sectional area of 1mm² (Fig. S2). Sensitivity tests were carried out using the Pod Kopište model biting at a 20° gape to quantify the impact of the presence / absence and variation in stiffness of the ligaments (50, 250, and 500 N/mm²) on the stress distribution (Figs. S3,S4). The stress values associated with each element of the mesh were exported into an element table for postprocessing.

Simulations

Four different MDA and FEA models were built: two "natural" models in which the specimens were modelled with their own musculature (PK_PK and PM_PM); and two "theoretical" models in which muscle PCSAs were swapped (i.e., PK_PM: the morphology of Pod Kopište with the musculature of Pod Mrčaru, and PM_PK the morphology of Pod Mrčaru with the musculature of Pod Kopište). These models allowed us to test the impact of changes in cranial morphology and muscle anatomy on the calculated bite force by comparing the results of the simulations. The calculated bite force in each of the four models was calculated for 10 different gape angles, from 0° (closed jaw) to 45° (maximum gape typically attained at the onset of fast closing, [START_REF] Taverne | Relationships between dietary breadth and flexibility in jaw movement: a case study of two recently diverged insular populations of Podarcis lizards[END_REF], thereby effectively varying the size of the prey item. For each MDA model, the conversion rate of the total muscle force into bite force was calculated at every gape tested by dividing bite force by the total muscle force and was used to assess the efficiency of a muscle to translate intrinsic muscle force into bite force. The total amount of change in bite force between the two natural models was quantified as it provides direct information on how much the mechanical properties of the masticatory system of the two specimens differ. The amount of change in bite force due to the change in muscle anatomy or to the change in cranial morphology alone was also quantified. This allowed us to estimate the relative contribution (in %) of changes in cranial morphological or in muscular anatomy to generate the observed differences in bite force between the specimens. The role of each muscle in bite force generation was assessed by running MDAs with only one muscle activated at a time.

The results of the MDA models were then imported into ANSYS, after having been multiplied by five to obtain realistic absolute loads corresponding to in vivo bite forces (i.e. the calculated bite forces were five times lower than the in vivo forces). For each model, FEAs were performed at three gape values that mimicked three ecologically relevant scenarios: biting a thin object like a leaf (0°), biting a prey of medium size (20°) and biting a large prey item or a conspecific during for example male-male interactions (45°). Values of stress in each skull element were saved. Specifically, we focused on the von Mises' stress (vMS) and recorded the mean stress values for the whole skull and the variation in the stress along the skull. To do so we divided the geometry into 10 segments of equal length (Fig. 1) and averaged the vMS values within each section. We then calculated the 'advantage' of the PM morphology relative to the PK morphology as the percentage of increase/decrease of the mean vMS accumulation along the skull. Finally, we estimated the homogeneity of the distribution of stress by calculating the Shannon's entropy [START_REF] Shannon | A mathematical theory of communication[END_REF] based on the density distribution of vMS values. A higher entropy reflects a low redundancy and a high disparity in the data. Thus, the 'advantage' provided by the PM morphology was estimated as the proportional difference in entropy between the two models.

Results

Morphological differences

Lizards from the two islets differ in snout-vent length (Pod Kopište: 62.3 ± 3.7 mm for males; 58.2 ± 3.0 mm for females; Pod Mrčaru: 68.3 ± 3.2 mm for males; 62.72 ± 2.9 for females), head length (Pod Kopište: 14.9 ± 0.7 mm for males; 13.0 ± 0.6 mm for females; Pod Mrčaru: 16.22 ± 0.8 mm for males; 13.8 ± 0.6 for females), and bite force (Pod Kopište: 17.3 ± 4.6 N for males and 9.6 ± 1.9 N for females; Pod Mrčaru: 27.01 ± 3.2 N for males and 14.8 ± 2.7 N for females).

The MANOVA including the first nine axes of the PCA carried out on skull shape (cumulative variance: 85.2%) revealed an effect of island (Wilks' lambda = 0.22; F 1,25 = 6.861; P < 0.001), especially on PC3 (variance = 9.3%; F 1,25 = 26.48; P < 0.001) and PC6 (variance = 3.3%; F 1,25 = 4.83; P = 0.037). The theoretical shape deformation along PC3 showed that the population of PM had a slightly shorter snout and a more curved quadrate (Figure S5). The MANOVA including the first 9 axes of the PCA carried out on mandible shape (cumulative variance: 85.0%) also revealed an effect of island (Wilks' lambda = 0.40; F 1,25 = 2.86; P = 0.029), especially along PC1 (variance = 24.4%; F 1,25 = 7.50; P = 0.011). The theoretical shape deformation along PC1 showed that the mandible of the PM population was overall more ventrally curved (Figure S5). The coronoid was medially and caudally thicker, the lateral crest, which serves as the insertion site of external adductors was wider, the retroarticular process shorter, and the joint surface area larger.

The MANCOVA carried out on the muscle variables detected significant differences between islands (Table 2). All variables (muscle mass, fiber length, and PCSA) were greater in specimens from Pod Mrčaru. When accounting for body size, specimens from Pod Mrčaru had heavier external adductors and pterygoid muscles, shorter fibers in the AMP, and a greater PCSA of the external adductors and the AMP (Table 3).

Mechanical basis of variation in bite force

The MDAs showed that models with the musculature of the PM specimen always resulted in higher bite forces, and that models with the shape of PM individuals resulted in higher bite forces irrespective of variation in muscle architecture (Fig. 2, Table 4). This held irrespective of the gape tested, and was even exacerbated with an increase in gape. Simulated bite force reached two peaks in the four models, at closed gapes (0° or 5° in PK_PM) and wide gapes (45°) and was minimal at a gape of around 30° / 35°. Knowing that the total bite force of the PK and the PM individuals modeled in the present study was 19.1 N and 25.0 N respectively, the simulations showed that the conversion rate of the PM morphology is always higher than in the PK morphology (Fig. 2). The total amount of change in simulated bite force between the two natural models varied between 30.5% (0°) and 37.5% (40°). On average, the relative contribution from the change in skull shape to total variation in bite force increased with gape and varied between 2.41% (5°) and 12.49% (45°), the remainder being explained by variation in muscle architecture (Fig. 3).

In all models, the external adductors contributed the most to bite force (on average approximately 50%), then the pseudotemporalis (35%), and the pterygoids (15%) (Fig. S6, Fig. S7, Table S3). The relative contribution of the muscles varied with gape; the contribution of the external adductors decreased with gape, while that of the pterygoid increased. The pseudotemporalis group had the highest conversion rate of muscle force into bite force (approximately 30% on average), followed by the external adductors (22%) and the pterygoids (12%) (Figure S7). The conversion rate of the adductors decreased with gape in all models, while the conversion rate of the pterygoids increased, and that of the pseudotemporalis decreased from gape 0° to 30° and increased from 30° to 45°. Compared to the morphology of PK, the PM morphology was associated with higher conversion rates of the adductors for all gapes, of the pseudotemporalis for gapes wider than 10°, and for the pterygoids for gapes wider than 30°.

Structural performance of the skull

The comparison of the two natural models showed that, on average, the von Mises stress magnitude (vMS) in the cranium of PM_PM model was higher than in that in the PK_PK model (from 6% to 17%) (see Fig. 4). Although PM_PM showed lower vMS magnitudes (except for the two most posterior sections of the skull) at a 0° gape, more stress was detected at 20° and 45°. Interestingly, the amount of vMS observed significantly increased with gape in PM_PM, whereas this was not the case in PK_PK. Specifically, the increase in vMS in the PM_PM natural model affected the whole skull (Fig. 4). An increase of 2% to 6% was also observed in the overall entropy in the PM_PM model (Table 5).

The PM morphology enabled the mean vMS to decrease by 15% compared to a PK skull model with the PM musculature. The advantage provided by the PM morphology, estimated in percentage of loss in vMS, increased at low gape and was more pronounced in the anterior part of the skull than the posterior part (see Fig. 5). Indeed, the advantage of a PM morphology was 13% to 41% in the snout, 9% to 20% in the region of the jugal, and became null or slightly negative in the two most posterior regions (where stress magnitude increased). Specifically, the unloading allowed by the PM skull morphology concerned the maxilla, the vomer, the frontal, the anterior part of the pterygoid, the epipterygoid, the parietal, and the quadrate. On the contrary, the areas that showed higher stress magnitudes included the wings of the prootic and the posterior tip of the pterygoid. Although less obvious when reaching a 45° gape, this pattern of stress distribution was similar at all gapes (Fig. 5). On average, the vMS magnitude increased (from 19% at 0° to 27% at 45°) when a PM musculature rather than a PK musculature was applied to a given morphology. In all cases, the overall level of entropy increased when a PM musculature rather than a PK musculature was applied on a given morphology (from 6% at 0° to 7% at 45°). The PM morphology provided only a small advantage in reducing overall entropy (differences < 5%; Table 5).

Discussion

A significant difference in cranial anatomy and muscle architecture exists between the individuals from Pod Kopište and Pod Mrčaru. The latter showing stronger jaw adductors associated with their larger size and omnivorous diet. Lizards from Pod Mrčaru are bigger than those from Pod Kopište, but the variation in muscle cross-sectional area is not only due to differences in size. This is in line with the results of a previous study comparing multiple populations of Podarcis lizards [START_REF] Taverne | From micro to macroevolution: drivers of shape variation in an island radiation of Podarcis lizards in the Adriatic[END_REF]. Our results suggest that other factors may impact the differences in muscle architecture. For example, lizards from Pod Mrčaru have stronger pterygoid muscles which have a low moment arm at low gape [START_REF] Herrel | Sexual dimorphism of head size in Gallotia galloti: testing the niche divergence[END_REF]. As plant consumption typically involves biting at low gape, this suggests that diet might not underpin the differences in the cross-sectional area of pterygoid muscles observed between lizards from the two populations. We rather suggest that biting in the context of intraspecific interactions may better explain this difference. The density of lizards on Pod Mrčaru is roughly five times higher than on Pod Kopište [START_REF] Herrel | Rapid large-scale evolutionary divergence in morphology and performance associated with exploitation of a different dietary resource[END_REF][START_REF] Vervust | The mystery of the missing toes: extreme levels of natural mutilation in island lizard populations[END_REF], increasing the probability of encounters and aggressive interactions. Biting congeners involves biting at large gape which may drive the observed differences in the pterygoid muscle [START_REF] Donihue | Feed or fight: testing the impact of food availability and intraspecific aggression on the functional ecology of an island lizard[END_REF]. Populations also differed in cranial and mandibular shape. Our results show that lizards on Pod Mrčaru have a mandible that is more ventrally curved, with a wider lateral insertion area for adductor muscles, and a more robust coronoid. The results for skull morphology are less clear and mainly show differences in the curvature of the quadrate and a shorter snout in animals from Pod Mrčaru. The fact that the morphological differences in the cranium are less clear is likely a consequence of the multiple functions that the cranium fulfills (e.g., protection of the central nervous system and of the sensory organs), and the associated constraints and trade-offs. Non-adaptive processes including founder effects or genetic drift after introduction may also have driven the phenotypic differences between the two populations observed here [START_REF] Kolbe | Founder effects persist despite adaptive differentiation: a field experiment with lizards[END_REF].

Our modeling results demonstrate that differences in bite force between the two morphotypes are principally driven by the difference in the total muscle PCSA (Table 2). The in-silico estimates of bite force generated by the MDA simulations suggested that the lizards from Pod Mrčaru bite harder than those of Pod Kopište. Yet, subtle differences in cranium and mandible shape also contributed to the observed variation in bite force between the two populations. The comparison of the natural and hybrid models showed that the morphology of Pod Mrčaru enables greater bite forces to be generated in all cases. Moreover, the correlation between bite force and gape parallels the observed relation between conversion rate and gape (Fig. 2). Hence, variation in bite force at different gape angles depends on the proportion of muscle force converted into bite force, which in turn is impacted by head shape. This is further supported by the relative contribution of musculature and morphology to the differences in bite force. Consequently, the relatively subtle shape changes in the mandible and cranium between these two recently diverged populations are responsible for a considerable part of the variation in bite force. This is achieved by modulating the lever-arms associated with the muscle bundles that depend on the variation in skull and mandible geometry as well as differences in muscle size and cross-sectional area.

The comparison of the natural models PK_PK and PM_PM highlighted differences in structural performance between the phenotypes of the two populations and revealed that the phenotype of Pod Mrčaru accumulates more stress. Given that the PM model includes a greater total muscle PCSA, this is expected. Moreover, our models confirm that applying greater forces on a given morphology increases the level of VM stress over the skull. But the magnitude of differences in stress due to musculature is greater (19% to 27%) than that observed between the natural models (from 6% to 17%) suggesting that morphology itself might help dissipate stress. Indeed, the PM morphology conferred a clear advantage in dissipating stress compared to the PK morphology (Fig. 5). Similarly, the comparison of the natural models revealed that the phenotype of the Pod Mrčaru individuals had a greater level of overall entropy compared to that of Pod Kopište (from 2% at 0° to 6% at 45°). Again, this is mainly due to an increase in the total muscle force with the application of a PM musculature to a given morphology leading to an overall increase in entropy (from 6% at 0° to 7% at 45°). Interestingly, our results suggest that an increase in the total muscle force is accompanied with an increased disparity in the distribution of the stress over the skull, yet the phenotype of Pod Mrčaru maintained a more homogenous stress distribution.

Unexpectedly, the PM morphology was more advantageous than the PK morphology in producing bite force at wide gapes (as revealed by the MDA simulations), whereas it is more advantageous in dissipating stress at lower gapes (as revealed by the FEA simulations). Although this might first seem contradictory, we propose that these two results are not incompatible. We rather think that selection primarily acts upon the performance of the musculoskeletal system (i.e., bite force) within the limits imposed by the intrinsic capacity of the system to dissipate the associated stress and strain. Therefore, the MDA results suggest that the evolution of skull morphology towards a Pod Mrčaru phenotype enabled an optimization of bite force at wide gapes. This is likely driven by intraspecific competition, enhanced by a five-fold increase in population density on Pod Mrčaru [START_REF] Herrel | Rapid large-scale evolutionary divergence in morphology and performance associated with exploitation of a different dietary resource[END_REF][START_REF] Vervust | The mystery of the missing toes: extreme levels of natural mutilation in island lizard populations[END_REF]. The consumption of hard and fibrous material like plant matter requires repeated and frequent biting [START_REF] Herrel | Modularity complexity of the feeding repertoire in scincid lizards[END_REF] to reduce the item into smaller and more easily digestible bite-size pieces [START_REF] Bjorndal | Digestive fermentation in herbivores: effect of food particle size[END_REF][START_REF] Bjorndal | Body size and digestive efficiency in an herbivorous freshwater turtle: advantages of small bite size[END_REF]. The structural advantage provided by the PM morphology at low gapes might be important in this context. Indeed, the repeated loading occurring at low gapes experienced by the skull while feeding on plants is likely to be particularly constraining and might represent an important selective agent driving variation in skull shape. A scenario in which head shape in insular Podarcis lizards is initially driven by intraspecific competition, and maintained by functional demands associated with resource use, would be in line with the results of Donihue and collaborators (2016) concerning intraspecific variation in morphology and performance in another Podarcis species.

Previous comparative studies on the colonization of Podarcis melisellensis and Podarcis siculus in the Adriatic archipelago and their morphological evolution in relation to variation in ecological contexts [START_REF] Taverne | Proximate and ultimate drivers of variation in bite force in insular lizards Podarcis melisellensis and Podarcis siculus[END_REF][START_REF] Taverne | From micro to macroevolution: drivers of shape variation in an island radiation of Podarcis lizards in the Adriatic[END_REF] suggested that relationships between form and function were similar at different levels of integration, whether between isolated populations or between species, and even when accounting for phylogenetic relationships. Specifically, convergent evolution of muscle architecture and head morphology were found in similar ecological contexts. Here, we demonstrate that these form-function relationships can evolve on ecological time scales and be associated with an optimization of the mechanical output of the masticatory system. Consequently, subtle morphological variation may accumulate over time and give rise to macroevolutionary patterns.

Future directions

Further efforts might help make the finite element models even more realistic and circumvent some of the current limitations. First, the simulations revealed an unrealistic accumulation of vMS at the very back of the skull (Figs. 4,5) at the most lateral region of the neurocranium (paraoccipital process of the exoccipital). In the present study, the joint between the lateral process of the exoccipital and the quadrate was modeled as connected and fixed, hence possibly explaining why abnormal stress concentrations are located in this area instead of being more homogeneously distributed. It would be useful to model a ligamentous connection between these bones allowing a better stress dissipation.

Second, bone tissue was modeled as an isotropic material here, although it is known to be anisotropic [START_REF] Moazen | Mechanical properties of calvarial bones in a mouse model for craniosynostosis[END_REF]. A precise quantification of the bony material properties throughout the skull (Young's modulus) should be undertaken using a nano-indentation approach. A third way to gain accuracy would be to model the sutures and the osteoderms, since they have been suggested to impact how loads are distributed and dissipated [START_REF] Moazen | Assessment of the role of sutures in a lizard skull: a computer modelling study[END_REF][START_REF] Costantini | The relationship between cephalic scales and bones in lizards: a preliminary microtomographic survey on three lacertid species[END_REF][START_REF] Jones | The biomechanical role of the chondrocranium and sutures in a lizard cranium[END_REF].

Conclusions

Our results show significant differences in skull form and function in individuals of two populations of P. siculus. Following the introduction into a novel environment, the population on Pod Mrčaru rapidly evolved a largely herbivorous diet. Our analyses of shape and function provide a mechanistic underpinning of the link between relatively subtle differences in morphology and the observed changes in ecology. Our results demonstrate that relationships between form and function may drive variation in fitness-relevant performance traits resulting in changes in trophic ecology over a relatively short time scale.

Figure legends

Figure 1: MDA model design illustrating the muscle bundles used of each muscle and their wrapping. A: left lateral view of the skull, B: dorsal view, C: caudal view, D: ventral view. The jaw opener is represented in yellow, the external adductors in shades of blue, the pseudotemporalis in red and orange and the pterygoids in green. E: segments of the skull used to investigate the variation in stress accumulation over the structure (here represented by a ventral view) in subsequent FEA. Note that the bite points and the joints are indicated by white stars and white circles, respectively.

Figure 2 :

 2 Figure2: A: relationship between bite force (in Newtons), estimated by the MDA for the four models, and the gape (in degrees). B: relationship between the conversion rate of the total muscle force into bite force (in percentage) and the gape (in degrees) for all models. The two natural models (PK_PK: including the morphology and the musculature of Pod Kopište, PM_PM: including the morphology and the musculature of Pod Mrčaru) are represented by full lines, whereas the two hybrid models (PK_PM: including the morphology of Pod Kopište and the musculature of Pod Mrčaru, PM_PK: including the morphology of Pod Mrčaru and the musculature of Pod Kopište) are represented by dashed lines.

Figure 3 :

 3 Figure3: Contribution (in percent) of the differences in musculature (in light gray) and skull shape (in dark gray) in explaining the total difference in calculated bite force between individuals from the two islands (Pod Kopište and Pod Mrčaru). The contributions were estimated for every gape tested (in degrees) by comparing the change in calculated bite force induced by the change in musculature or morphology alone (theoretical models) with the change in calculated bite force between the two natural models.

Figure 4 :

 4 Figure4: von Mises stress distribution in the two natural models (Brown lines: PK_PK, green lines: PM_PM) for the three gapes tested (full lines: 0°, dashed lines: 20°, dotted lines: 45°). Stresses are plotted relative to the virtual skull section number (i.e. position along the antero-posterior axis for which results were averaged; see methods). The left lateral, ventral, caudal and dorsal views of the

Figure 5 :

 5 Figure 5: Consequences of skull shape change on the von Mises stress distribution. The graph represents the mean advantage (in %) provided by the PM morphology along the skull (full line: 0°, dashed line: 20°, dotted line: 45°).Here the models PK_PM and PM_PM (PK / Pod Kopište and PM / Pod Mrčaru in the figure) were compared for each of the three gapes tested (see FigureS8for the alternative combination PM_PK vs. PK_PK). Hotter colors are associated with higher von Mises stress magnitudes (in MPa).

Acknowledgements

We would like to thank two anonymous reviewers for constructive comments on the manuscript. This project was funded by a National Geographic Explorer Grant to AH and by a doctoral scholarship provided by Sorbonne Université (Idex SUPER SU-17-R-DOC-03). We also are grateful to the doctoral school ED227, to the Society for Experimental Biology, to the Société des Amis du Muséum and to the Carl Gans Collections & Charitable Fund for providing precious financial support that enabled MT to do field work and to attend congresses in the context of this study.

Data Accessibility

doi:10.5061/dryad.866t1g1vk sharing link: https://datadryad.org/stash/share/uy5-FARWE09mXDedc1SGivCLqz8aqQU6fK16_Ak9L_Y

Author contributions

MT, MJF and AH designed the study. MT, HD, RB, DL, ZT, ACF and AH generated the raw data. MT, PW and MJF built the models. MT analyzed the data and ran the models. All authors were involved in the interpretation of the results, the writing, and the revision of the manuscript.