Claire Mérot

Vincent Debat

Yann Le Poul

Richard Merrill

Russell Naisbit

Adélie Tholance

Chris Jiggins

Mathieu Joron

Hybridization and transgressive exploration of colour pattern and wing morphology in Heliconius butterflies

Keywords: diversification, mimicry, wing shape, colour pattern, Lepidoptera, Transgression, geometric morphometrics Wing shape -Colour pattern -Transgression -Geometric morphometry -Lepidoptera -Mimicry -Diversification

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Recent studies on phenotypic diversification reveal a more ambivalent role for hybridization than was previously recognised [START_REF] Mallet | Hybridization, ecological races and the nature of species: empirical evidence for the ease of speciation[END_REF][START_REF] Hegarty | Invasion of the hybrids[END_REF][START_REF] Abbott | Hybridization and speciation[END_REF]Yakimowski & Rieseberg, 2014;[START_REF] Meier | Ancient hybridization fuels rapid cichlid fish adaptive radiations[END_REF]. Hybridization has classically been viewed as a force that homogenizes populations, acting against genetic differentiation, and preventing speciation [START_REF] Mayr | Animal Species and Evolution[END_REF][START_REF] Felsenstein | Skepticism towards Santa Rosalia, or why are there so few kinds of animals?[END_REF]. However, many species continue to hybridize throughout the process of speciation (Kronforst et al., 2006b;[START_REF] Wiens | Rapid diversification, incomplete isolation, and the "speciation clock" in North American salamanders (Genus plethodon): Testing the hybrid swarm hypothesis of rapid radiation[END_REF][START_REF] Mallet | Natural hybridization in heliconiine butterflies: the species boundary as a continuum[END_REF][START_REF] Kulathinal | The genomics of speciation in drosophila: diversity, divergence, and introgression estimated using low-coverage genome sequencing[END_REF][START_REF] Martin | Genome-wide evidence for speciation with gene flow in Heliconius butterflies[END_REF], and increasing evidence suggests that a limited level of hybridization can instead represent a source of novel genetic variation, with possibly greater effects than mutation [START_REF] Grant | Phenotypic and genetic effects of hybridization in Darwin's finches[END_REF][START_REF] Abbott | Hybridization and speciation[END_REF]Dittrich-Reed & Fitzpatrick, 2013).

Hybridization produces novel allelic combinations in hybrids, which may result in transgressive phenotypes that lie beyond natural variation expressed in parental populations [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF]Stelkens et al., 2009). Transgression has been reported for traits as diverse as mandible shape in mice [START_REF] Renaud | Mandible shape in hybrid mice[END_REF]2012), tolerance to temperature in copepods [START_REF] Pereira | Ecological novelty by hybridization: experimental evidence for increased thermal tolerance by transgressive segregation in Tigriopus californicus[END_REF], and transcription level in fish [START_REF] Czypionka | Transcriptome changes after genome-wide admixture in invasive sculpins (Cottus)[END_REF]. When novel phenotypes are deleterious, reduced hybrid fitness forms a post-zygotic barrier, which can enhance pre-mating isolation through reinforcement, leading to species differentiation and ultimately diversification (Butlin, 1987a;Butlin, 1987b;Servedio & Noor, 2003;Smadja & Butlin, 2006). When novel phenotypes bear an adaptive value, in contrast, hybrids may backcross to parents and thus contribute to increasing the available genetic and phenotypic variation. This may influence the evolutionary trajectories of phenotypes in parental species [START_REF] Selz | Relaxed trait covariance in interspecific cichlid hybrids predicts morphological diversity in adaptive radiations[END_REF], allowing, for instance, the colonization of new ecological niches [START_REF] Nolte | An invasive lineage of sculpins, Cottus sp (Pisces, Teleostei) in the Rhine with new habitat adaptations has originated from hybridization between old phylogeographic groups[END_REF][START_REF] Hermansen | Hybrid speciation in sparrows I: phenotypic intermediacy, genetic admixture and barriers to gene flow[END_REF][START_REF] Stemshorn | Rapid formation of distinct hybrid lineages after secondary contact of two fish species (Cottus sp.)[END_REF][START_REF] Stelkens | Hybridization facilitates evolutionary rescue[END_REF][START_REF] Meier | Ancient hybridization fuels rapid cichlid fish adaptive radiations[END_REF] or increasing within-and between-species phenotypic diversity [START_REF] Gilbert | Adaptive novelty through introgression in Heliconius wing patterns: evidence for shared genetic "tool box" from synthetic hybrid zones and a theory of diversification[END_REF][START_REF] Mallet | Rapid speciation, hybridization and adaptative radiation in the Heliconius melpomene group[END_REF].

Hybrid transgression can have diverse genetic origins, including epistasis (Eshed & Zamir, 1996), overdominance [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF], complementary action of additive alleles between multiple quantitative trait loci (Devicente & Tanksley, 1993;[START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF][START_REF] Rieseberg | The genetic architecture necessary for transgressive segregation is common in both natural and domesticated populations[END_REF]Albertson & Kocher, 2005;Stelkens & Seehausen, 2009;Stelkens et al., 2009) or even epigenetics [START_REF] Shivaprasad | Extraordinary transgressive phenotypes of hybrid tomato are influenced by epigenetics and small silencing RNAs[END_REF]. Theoretical predictions and empirical observations suggest that the degree of transgression depends on the phenotypic divergence between parental species. In particular, species that are phenotypically similar may be more likely than phenotypically more distinct species to accumulate complementary additive alleles which, summing up in hybrids, lead to transgression [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF].

We here investigate the relationship between hybridization, adaptation and diversification in Heliconius butterflies. Hybridization is frequent between recently diverged species within the genus and results in a large fraction of the genome being shared between closely-related species [START_REF] Mallet | Natural hybridization in heliconiine butterflies: the species boundary as a continuum[END_REF][START_REF] Martin | Genome-wide evidence for speciation with gene flow in Heliconius butterflies[END_REF]. Hybridization has been suggested as a critical factor in the evolution of Heliconius phenotypic diversity because inter-specific laboratory crosses generate a wide variety of wing patterns, a subset of which match natural patterns from other geographic areas [START_REF] Gilbert | Adaptive novelty through introgression in Heliconius wing patterns: evidence for shared genetic "tool box" from synthetic hybrid zones and a theory of diversification[END_REF]. The importance of gene flow in the evolution of Heliconius is supported in nature by molecular evidence showing adaptive introgression of mimicry loci [START_REF] Heliconius | Butterfly genome reveals promiscuous exchange of mimicry adaptations among species[END_REF]Pardo-Diaz et al., 2012), sometimes followed by a shuffling between major-effect loci enhancing overall pattern diversity [START_REF] Wallbank | Evolutionary Novelty in a Butterfly Wing Pattern through Enhancer Shuffling[END_REF][START_REF] Enciso-Romero | Evolution of novel mimicry rings facilitated by adaptive introgression in tropical butterflies[END_REF]. Here we use hybrids from laboratory crosses in Heliconius to investigate the effect of hybridization on wing shape and colour pattern, two traits involved in Müllerian mimicry but displaying a different range of natural variation and likely based on different genetic architecture.

Wing colour pattern is a strikingly diversified trait in Heliconius. Wing pattern is under strong natural selection through its role as a warning signal involved in Müllerian mimicry (Mallet, 1989;[START_REF] Merrill | Disruptive ecological selection on a mating cue[END_REF]. The genetic architecture of wing pattern is shared across species and typically involves a limited number of large-effect loci switching on or off the presence of black shutters or colour patches in different compartments of the wing [START_REF] Gilbert | Adaptive novelty through introgression in Heliconius wing patterns: evidence for shared genetic "tool box" from synthetic hybrid zones and a theory of diversification[END_REF][START_REF] Joron | A conserved supergene locus controls colour pattern diversity in Heliconius Butterflies[END_REF]Kronforst et al., 2006a;[START_REF] Baxter | Butterfly speciation and the distribution of gene effect sizes fixed during adaptation[END_REF][START_REF] Ferguson | Characterization of a hotspot for mimicry: assembly of a butterfly wing transcriptome to genomic sequence at the HmYb/Sb locus[END_REF][START_REF] Huber | Conservatism and novelty in the genetic architecture of adaptation in Heliconius butterflies[END_REF][START_REF] Wallbank | Evolutionary Novelty in a Butterfly Wing Pattern through Enhancer Shuffling[END_REF][START_REF] Enciso-Romero | Evolution of novel mimicry rings facilitated by adaptive introgression in tropical butterflies[END_REF]. Between species differing in colour pattern, hybridization usually reshuffles parental alleles, forming hybrids with a novel colour pattern that consists of a combination of elements of the parental patterns [START_REF] Gilbert | Adaptive novelty through introgression in Heliconius wing patterns: evidence for shared genetic "tool box" from synthetic hybrid zones and a theory of diversification[END_REF][START_REF] Naisbit | Mimicry: developmental genes that contribute to speciation[END_REF]. Such phenotypes generally do not sufficiently match either of the parental patterns that are recognised and avoided by predators, and so hybrids suffer from higher predation [START_REF] Merrill | Disruptive ecological selection on a mating cue[END_REF]. This process strongly contributes to isolation between closely-related species [START_REF] Jiggins | Reproductive isolation caused by colour pattern mimicry[END_REF][START_REF] Jiggins | Ecological speciation in mimetic butterflies[END_REF][START_REF] Mallet | Shift happens! Shifting balance and the evolution of diversity in warning colour and mimicry[END_REF]. Yet, if the novel recombined pattern resembles other local patterns, or if hybrids are produced at a nonnegligible frequency, hybridization may be a way to cross the fitness valley between adaptive peaks of mimicry, and by this process, hybridization may contribute to phenotypic diversification [START_REF] Gilbert | Adaptive novelty through introgression in Heliconius wing patterns: evidence for shared genetic "tool box" from synthetic hybrid zones and a theory of diversification[END_REF]Mavarez & Linares, 2008;[START_REF] Mallet | Rapid speciation, hybridization and adaptative radiation in the Heliconius melpomene group[END_REF][START_REF] Merrill | The diversification of Heliconius butterflies: what have we learned in 150 years?[END_REF].

Wing shape and size have received less attention than colour pattern because they appear more homogenous throughout the clade. However, experiments on flight and quantification of wing shape suggest that shape similarity between co-mimics is also favoured by selection [START_REF] Srygley | Locomotor mimicry in Heliconius butterflies: contrast analyses of flight morphology and kinematics[END_REF][START_REF] Jones | Wing shape variation associated with mimicry in butterflies[END_REF][START_REF] Mérot | Refining mimicry: phenotypic variation tracks the local optimum[END_REF]. Contrary to colour pattern, for which the genetic architecture is rather simple, shape is generally polygenic, involving interactions between several loci [START_REF] Zimmerman | Quantitative trait loci affecting components of wing shape in Drosophila melanogaster[END_REF]Navarro & Klingenberg, 2007;[START_REF] Klingenberg | Evolution and development of shape: integrating quantitative approaches[END_REF][START_REF] Pitchers | The power of a multivariate approach to genome-wide association studies: an example with Drosophila melanogaster wing shape[END_REF]. This makes it difficult to predict the effects of hybridization.

Geographical races of Heliconius allow us to contrast cases differing in the nature and degree of divergence in parental morphology. Here we consider two pairs of species belonging to sister clades: in Panama, H. cydno chioneus and H. melpomene rosina and in Peru, H. timareta thelxinoe and H. melpomene amaryllis (Fig. 1A). H. timareta and H. cydno belong to the same phylogenetic lineage, which diverged from H. melpomene about 2 million years ago [START_REF] Kozak | Multilocus Species Trees Show the Recent Adaptive Radiation of the Mimetic Heliconius Butterflies[END_REF]. Therefore, the two pairs involve similar duration of phylogenetic divergence. However, phenotypic distances differ: the Peruvian H. t. thelxinoe and H. m. amaryllis are co-mimics and exhibit a very similar colour pattern with a red forewing patch and a yellow hindwing bar. In contrast, the Panamanian H. c. chioneus and H. m. rosina differ in colour and pattern, mimicking two very distinct distantly-related Heliconius species (Fig. 1A). As a result, wing colour pattern, and likely wing shape, are under convergent selection in the Peruvian co-mimetic pair (H. m. amaryllis/H. t. thelxinoe) and under divergent selection in the Panamanian non-mimetic pair (H. m. rosina/H. c. chioneus). We can therefore investigate the effect of parental phenotypic divergence upon the extent of hybrid transgression.

Transgression is the occurrence of extreme or novel phenotypes relative to the parental phenotypes but, in the literature, it has been quantified in different ways, corresponding to subtly different definitions. Transgression has been considered (1) as a deviation from a linear intermediate between parental phenotypes, i.e. the mean hybrid value is not included in the interval between the two parental means [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF][START_REF] Renaud | Modularity as a source of new morphological variation in the mandible of hybrid mice[END_REF][START_REF] Renaud | Impact of Hybridization on Shape, Variation and Covariation of the Mouse Molar[END_REF], (2) or as the proportion of hybrids falling outside of the range of the parental phenotypes [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF][START_REF] Parsons | Hybridization Promotes Evolvability in African Cichlids: Connections Between Transgressive Segregation and Phenotypic Integration[END_REF]Holzman & Hulsey, 2017),

(3) or as the expansion of the hybrid range beyond the range of the parental phenotypes (Stelkens et al., 2009;Holzman & Hulsey, 2017;[START_REF] Husemann | Body shape differences in a pair of closely related Malawi cichlids and their hybrids: Effects of genetic variation, phenotypic plasticity, and transgressive segregation[END_REF].

We explore variation in wing size, shape and pattern, across parents, F1 and backcrosses in the two species pairs, to address the following questions: (1) What are the hybrid phenotypes and do they differ from parental morphologies? (2) To what extent do hybrids exhibit transgressive phenotypes and thereby explore novel wing morphologies? (3) Is the level of transgression dependent upon the phenotypic distance between parental species? In doing so, we explored three ways of measuring transgression and discuss the implications of choosing a given method (Fig. 1BCD).

Methods

Butterfly collection and crosses

Heliconius cydno chioneus (C) and Heliconius melpomene rosina (Mr) were collected from Gamboa, Republica de Panama. Heliconius timareta thelxinoe (T) and Heliconius melpomene amaryllis (Ma) were collected around Tarapoto, Peru (Fig. 1A). Field-caught mated females were used to establish stock populations of each taxon. For each sympatric pair, controlled crosses were performed to generate families of F1 hybrids (F1 C/Mr in Panama, F1 T/Ma in Peru) and backcross hybrids from crosses to both parental species (BC, backcross towards H. c. chioneus, BMr, backcross towards H. m. rosina; BT: backcross towards H. t. thelxinoe, BMa backcross towards H. m. amaryllis, Fig. 1A). Full details of the controlled crosses used here can be found elsewhere [START_REF] Naisbit | Mimicry: developmental genes that contribute to speciation[END_REF][START_REF] Merrill | Pervasive genetic associations between traits causing reproductive isolation in Heliconius butterflies[END_REF][START_REF] Mérot | Beyond magic traits: Multimodal mating cues in Heliconius butterflies[END_REF][START_REF] Merrill | Genetic dissection of assortative mating behavior[END_REF].

Measurement of wing size, shape and pattern

Images of ventral (v) and dorsal (d) forewings (FW) and hindwings (HW) were captured using a Nikon D90 digital camera with a Nikon micro 105/2.8G ED VR lens. We retained only specimens with undamaged wings, resulting in the collection of data, for shape and pattern respectively, from 67 and54 H In total, 20 FW landmarks and 18 HW landmarks were used to quantify wing shape.

Landmarks were placed at vein intersections and vein termini on the ventral side (Fig. S1) using TpsDig2 [START_REF] Rohlf | TpsDig, ver. 2. 16[END_REF]. Landmark coordinates were superimposed using a general Procrustes analysis [START_REF] Bookstein | Morphometrics tools for landmark data: geometry and biology[END_REF]Zelditch et al., 2004). A principal component analysis was applied to superimposed coordinates and the non-null PCs were used as shape descriptors in further analyses. Wing size was measured using centroid size [CS; see [START_REF] Bookstein | Morphometrics tools for landmark data: geometry and biology[END_REF]]. Dorsal wing colour pattern was analysed using Colour Pattern Modelling (CPM) (Le Poul et al., 2014). Briefly, wing outline was extracted individually from the background and, within this area, RGB colours were categorized into four colour classes (black, red, white or yellow) for each pixel. All individual wings were aligned by rotation, translation and scaling based on an iterative process maximizing overlap between pattern elements. This yields to wing surfaces characterized by a set of pixels with homologous position across specimens.

Principal component analysis was applied to the set of pixels and the PCs were used as colour pattern descriptors in further analyses.

Following [START_REF] Husemann | Body shape differences in a pair of closely related Malawi cichlids and their hybrids: Effects of genetic variation, phenotypic plasticity, and transgressive segregation[END_REF] and Stelkens et al. (2009), to remove sex-related effect, we performed a preparatory ANCOVA or MANCOVA with sex as factor and retained the residuals as variables for each of the six traits quantified (size, shape and colour pattern for the FW and HW). Because forewing and hindwing showed similar results, they were considered together, by summing the residuals for size, and by grouping all forewing and hindwing axes for either shape or pattern into two multivariate morphological matrixes (one for shape, one for pattern).

Differences in size between groups were investigated with a one-way analysis of variance (ANOVA), with p-values corrected for multi-test comparison following Benjamini and Hochberg (1995). Differences in shape and in pattern between groups were investigated with a one-way MANOVA with genotype as factor and respective PCs (shape or pattern) as dependent variables. PCA dimensionality was reduced following the Kaiser-Guttman criteria, i.e., by keeping all PCs whose eigenvalues were above the average eigenvalue [START_REF] Jackson | A user's guide to principal components[END_REF].

Phenotypic distances between parental species

To quantify phenotypic distances between parental species (T vs. Ma, C vs. Mr), we worked on a subset including only individuals from parental species. Between all possible pairs of parental specimens, for size, we computed absolute differences in wing size, and, for shape or colour pattern, we measured Euclidian distances in PCAs applied on the parental subsets, truncated following the Kaiser-Guttman criteria [START_REF] Jackson | A user's guide to principal components[END_REF]. To test whether parental distances differed between the mimetic and the non-mimetic pair, we applied a linear mixed model with individual distances as a variable, species pair (co-mimics/non-mimics) as a factor and identity of each sample as a random factor. Results were congruent when evaluated with Malahanobis distances between the four species (Table S1).

Analysis of hybrid phenotypic transgression

To quantify transgression in a meaningful morphospace, we translated the two morphological matrices (pattern and shape) into 12 subset morphospaces by applying a PCA on a subset with 3 genotypes, i.e., the target hybrid population H (F1 hybrids or either backcross hybrids) and its parental species (for instance the target hybrid population is BC and the parental populations are C and Mr), further truncated following the Kaiser-Guttman criteria [START_REF] Jackson | A user's guide to principal components[END_REF]. Transgression was then quantified with three indices corresponding to different definitions of transgressive segregation used in the recent literature for multivariate data.

 Transgression of the mean (Fig. 2B) [START_REF] Renaud | Modularity as a source of new morphological variation in the mandible of hybrid mice[END_REF]; [START_REF] Renaud | Impact of Hybridization on Shape, Variation and Covariation of the Mouse Molar[END_REF] proposed that in the absence of transgression, the average hybrid phenotype should be strictly intermediate between parental groups. In that case the sum of the distances between the mean hybrid and the two parents (𝑑(𝐻,𝑃 1 ) +𝑑(𝐻, ) should be equal to the distance between the parents . 𝑃 2 ) 𝑑(𝑃 1 ,𝑃 2 )

Transgression of the mean could thus be assessed as the departure of the mean hybrid phenotype from this expectation:

𝑇 𝑚 = 𝑑(𝐻,𝑃 1 ) + 𝑑(𝐻,𝑃 2 ) -𝑑(𝑃 1 ,𝑃 2 ) 𝑑(𝑃 1 ,𝑃 2 )
Where and are the Euclidian distances between the mean hybrid phenotype 𝑑(𝐻,𝑃 1 ) 𝑑(𝐻,𝑃 2 )

and the mean phenotype of the first and second parent, respectively, and is the

𝑑(𝑃 1 ,𝑃 2 )
Euclidian distance between the mean parental phenotypes.

 Range transgression (Fig. 2C) [START_REF] Parsons | Hybridization Promotes Evolvability in African Cichlids: Connections Between Transgressive Segregation and Phenotypic Integration[END_REF] considered transgression in relation to the frequency of transgressive hybrids, i.e., hybrids whose trait value exceeded a certain distance from the parental mean trait. Here, following traditional definition of transgression [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF], and the "proportional expansion index" defined in [START_REF] Holzman | Mechanical transgressive segregation and the rapid origin of trophic novelty[END_REF], we chose that distance to be the parental range. In the multivariate data, the range was defined as the maximal Euclidian distance in the morphospace between a parental phenotype (P i ) and its consensus (P 1 for parent 1 or P 2 for parent 2).

A hybrid was thus be considered transgressive when its Euclidian distance in the morphospace to both parental means was greater than the given parental range, i.e., if the distance between the hybrid H i and each of the parents (P 1 and P 2 ) meets the following two conditions:

𝑑(𝐻 𝑖 ,𝑃 1 ) > 𝑀𝑎𝑥 ( 𝑑(𝑃 𝑖 ,𝑃 1 )) & 𝑑(𝐻 𝑖 ,𝑃 2 ) > 𝑀𝑎𝑥 ( 𝑑(𝑃 𝑖 ,𝑃 2 ))
The index of range transgression (Tr) was calculated as the frequency of transgressive hybrids.

 Transgression strength (Fig. 2D) To compare the transgression index against a null model, each index was compared to a distribution of the transgression index built by simulating 1000 populations of hybrids built as a linear combination of the parental phenotype. More specifically, the trait of each simulated hybrid was calculated as the weighted mean trait of two randomly chosen parents:

Stelkens
𝐹 1 = 𝑃 1 + 𝑃 2 2 .
The simulated hybrid population was of the same sample size

; 𝐵𝑃 1 = 3𝑃 1 + 𝑃 2 4 ; 𝐵𝑃 2 = 𝑃 1 + 3𝑃 2 4
as the corresponding experimental population. Then, to evaluate the effect of unbalanced sample size, we also re-calculated each index 1000 times by drawing randomly 30 individuals in each experimental group (P1, P2, H) and reported the mean and standard deviation of this distribution. This was compared to the null distribution of the transgression index calculated again on 1000 populations of theoretically non-transgressive hybrids but with sample size fixed to 30.

To visualise the directions of maximal variation between parents and hybrids, as well as the direction of transgression, we used a between-group PCA (bg-PCA) [START_REF] Boulesteix | A note on between-group PCA[END_REF]Mitteroecker & Bookstein, 2011).

All analyses were performed in R 3.1 (R Core Team, 2014), using the libraries lme4 [START_REF] Bates | lme4: Linear mixed-effects models using Eigen and S4. R package version 1[END_REF], ade4 (Dray & Dufour, 2007) and Rmorph [START_REF] Baylac | Rmorph: a R geometric and multivariate morphometrics library[END_REF].

Results

Phenotypic distances in parental species

The parental species in both pairs differed significantly in size, shape and pattern (Table S1).

H. timareta thelxinoe and H. cydno were larger than H. melpomene (Fig. 2). They had a more elongated forewing with a shorter discal cell and a wider hindwing compared to H. melpomene (Fig. 3A-B, Fig. S3). Despite being co-mimetic, H. melpomene amaryllis and H. timareta thelxinoe differed in pattern, with subtle variation in the position and morphology of the colour patches (Fig. 3C, Fig. S4). Overall phenotypic distances were significantly smaller between mimetic H. m. amaryllis/H. t. thelxinoe than between non-mimetic H. m. rosina/H. cydno, by a factor of 1.3 for wing size (F 1,135 = 16, p < 0.001), 1.2 for wing shape (F 1,135 = 60, p < 0.001) and 3.1 for pattern (F 1,118 = 3896, p < 0.001), indicating not only similarity of pattern between co-mimics but also increased similarity of wing shape and size (Fig. S2, Table S1).

Hybrid phenotype and transgression

Wing size

For both the Peruvian and the Panamanian sets, forewing and hindwing centroid size of hybrids were intermediate between the smaller H. melpomene and the larger H. cydno/H. timareta (Fig. 2). Hybrid sizes were generally ordered according to the genetic contribution of parental species.

Transgression of the mean was null for all hybrid groups, except for BT, which exhibited a positive transgression of the mean (0.12). Few hybrids fell outside of parental wing size ranges. Range transgression and transgression strength were null, except for some very weak values in BMr (Tr = 0.01, Ts = 0.1) and BMa (Tr = 0.01, Ts = 0.01).

Wing shape

Hybrid wing shapes were significantly different from parental shapes, both for the mimetic and the non-mimetic crosses (Table S2). In both crosses, shape variation in forewing and hindwing was dominated by the differentiation of the parental species (bgPC1, 73-76% of shape variance; Fig. 3A,B). H. melpomene forewing and hindwing had a longer discal cell and a shorter, rounded distal region of the wing compared to H. cydno or to H. timareta. Along this axis of variation, driven by inter-specific parental divergence, wing shapes of hybrids were intermediate, staying within the overlapping parental ranges, with backcrosses closer to their majority parent. Along the second and third axes of the bgPCA describing 20-25% of variance, a slight transgressive effect was observed, with the group of hybrids deviating from the expected inter-parental position. This shifted position is associated with positive transgression of the mean wing shape, significantly higher than expected in simulated non-transgressive hybrids (Table 1).

Transgression of the mean was higher in backcrosses between the co-mimics (BMa: Tm = 0.27, BT: Tm = 0.44) than in backcrosses between non-mimics (BMr: Tm = 0.16, BC: Tm = 0.16), and this difference held true when controlling for parental distance ratio (Table S4). In contrast, the range transgression index was very low and similar between all groups, with less than 1% of hybrids outside the parental ranges (except F1 C/Mr : Tr = 0.03).

Transgression strength was also low and comparable between all groups when summing all dimensions (Ts ranges between 0.02 and 0.06), but significantly higher than in the simulated hybrids. The strength of transgression varied widely depending on the PC considered (Table S3), from 0, when hybrids did not exceed parental ranges, to 0.08-0.38, indicating that hybrids exceeded parental range by 8-38% in some directions of the shape space, for instance on the second or third PCs (representing about 10-15% of variance).

Wing colour pattern

For the mimetic Peruvian pair, parents and hybrids present the same general features, a black background with a red forewing patch and a yellow hindwing band. Despite this, there is significant variation in the morphology and the position of colour patches between parents and hybrids (Table S2, Fig. 3C, Fig. S4). Most of the between-group variation, represented by bg-PC1 (65% of variance), corresponds to differences between parents. The H. timareta thelxinoe forewing red patch is narrower, with a somehow zigzagging shape, and its yellow hindwing bar is narrower than that of H. melpomene amaryllis [START_REF] Mérot | Genetic differentiation without mimicry shift in a pair of hybridizing Heliconius species (Lepidoptera: Nymphalidae)[END_REF]. Along this axis of variation, hybrids are mostly intermediate and within the partially-overlapping parental range. Along the second and third bg-PC, hybrids displayed a transgressive pattern that is not observed in parents.

For the non-mimetic Panamanian pair, as previously described [START_REF] Naisbit | Mimicry: developmental genes that contribute to speciation[END_REF], F1 hybrids and some of the backcross individuals displayed a pattern that was never found in parents, recombining white and red in the forewing bar and showing a completely black hindwing (Fig. 3D, Fig. S4). This transgressive pattern falls between the two non-overlapping parental ranges on bg-PC1 and dominates the second bg-PC. Beyond the presence/absence of colour patches, there was also large variability within backcross groups in the position and size of colour patches (Fig. 3D, Fig. S4).

For hybrids from both sets of crosses, transgression of the mean colour pattern was significantly higher than expected in the simulated non-transgressive hybrids (Table 1). Transgression of the mean was higher in backcrosses between co-mimics (BMa: Tm = 0.41; BT: Tm = 0.42) than in backcrosses between non-mimics (BMr: Tm = 0.22; BC: Tm = 0.21), but this difference was not observed when controlling for the ratio of parental distances (Table S4). Range transgression is higher for wing pattern than for shape in both pairs and much higher for the non-mimetic pair than for the mimetic pair (Table 1). Only 5-10% of hybrids between the mimetic pair T/Ma fall outside of the parental ranges, while for the nonmimetic C/Mr pair, 71% to 83% of the backcrosses and 100% of F1 lie outside of the parental ranges. Transgression strength is also much higher in the non-mimetic pair than in the mimetic pair. Again, within the overall index, it is worth noting that, along certain directions of variation (several PCs explaining 2-10% of variance), mimetic hybrids exceed the parental range by up to 44% and non-mimetic hybrids by up to 330% (Table S3).

Discussion

Overall, the precise quantification of morphology in hybrids demonstrates the variable consequences of hybridization. Wing size is largely intermediate between parental phenotypes; hybrid wing shape deviates from the expected intermediate position between the parental shapes; colour pattern also deviates from an intermediate position and sometimes displays complete novelty compared to parents, particularly in the cross between non-mimetic species. The level of transgression therefore clearly varies among traits. This variation may reflect trait complexity, genetic architecture, and the phenotypic distances between parents. Furthermore, different methods of defining and quantifying transgression show different patterns and influence our conclusions with regard to the degree to which hybridization generates transgression and the influence of parental phenotypic similarity.

Transgression and novelty in hybrid phenotypes

In complex traits such as many morphological traits, transgression is expected and frequently observed in hybrids [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF][START_REF] Rieseberg | The genetic architecture necessary for transgressive segregation is common in both natural and domesticated populations[END_REF]Stelkens & Seehausen, 2009;Stelkens et al., 2009;[START_REF] Parsons | Hybridization Promotes Evolvability in African Cichlids: Connections Between Transgressive Segregation and Phenotypic Integration[END_REF][START_REF] Renaud | Modularity as a source of new morphological variation in the mandible of hybrid mice[END_REF][START_REF] Selz | Relaxed trait covariance in interspecific cichlid hybrids predicts morphological diversity in adaptive radiations[END_REF]Holzman & Hulsey, 2017;[START_REF] Husemann | Body shape differences in a pair of closely related Malawi cichlids and their hybrids: Effects of genetic variation, phenotypic plasticity, and transgressive segregation[END_REF][START_REF] Renaud | Impact of Hybridization on Shape, Variation and Covariation of the Mouse Molar[END_REF]. Here, we have shown that size, a univariate trait, shows nearly no transgressive effect, whereas multivariate traits such as shape and pattern show much more marked transgression. This may perhaps be due to an increasing range of possibilities for transgression when a greater number of dimensions of variation contribute to the phenotype.

The extent of hybrid transgression also depends on which definition of transgression is used. Indeed, indices that consider the mean and the range actually address different hypotheses.

Transgression of the mean evaluates the extent to which the mean phenotype of all hybrid offspring is different from an intermediate phenotype lying between the parental means. It actually tests a deviation from a genetic additive linear model, and whether a preferential direction of deviation is explored by the hybrids as a group but not individual novelty. In contrast, the range transgression, and transgression strength, test the extent to which the phenotypic range of hybrid exceeds parental ranges, in other words to what extent novel phenotypes are produced in hybrids.

The different indices provided congruent values of high transgression only for colour pattern

of non-mimetic hybrids. In contrast, wing shape (both pairs) and colour pattern (co-mimetic pair), showed intermediate values of transgression of the mean but almost null range transgression. In the latter case, this means that few of the hybrids display phenotypes outside the parental ranges, yet the majority of the hybrids deviate from a pure intermediate. Such deviation is missed by our range index, which supposes isotropy (equal variations in all directions), but it is captured by the second axis of the bg-PCA and by observing, in the index of transgression strength, a transgressive effect associated with some morphological directions. The fact that some directions of morphological transgression are occupied by a group of hybrids means that, although, individually, few hybrids represent new phenotypes, some hybrid morphologies are produced at a higher probability. This possibly reveals genetic interactions due to hybridization and a non additive/non-linear underlying genetic architecture and may also represent the most available and heritable directions of morphological variation.

Hybrid novelty and genetic architecture of wing morphology

The traits considered in this study rely on different genetic mechanisms, which likely explain the differences in transgressive effect. Patterns in hybrids between the non-mimetic Panamanian species, H. melpomene rosina and H. cydno chioneus, are totally new because they recombine discrete elements of parental colour pattern (such as black "shutters" or large colour patches) whose presence is controlled by unlinked loci of major effect in parental species (Martin et al.;Mallet, 1989;[START_REF] Gilbert | Adaptive novelty through introgression in Heliconius wing patterns: evidence for shared genetic "tool box" from synthetic hybrid zones and a theory of diversification[END_REF][START_REF] Naisbit | Mimicry: developmental genes that contribute to speciation[END_REF][START_REF] Jiggins | A genetic linkage map of the mimetic butterfly Heliconius melpomene[END_REF][START_REF] Joron | A conserved supergene locus controls colour pattern diversity in Heliconius Butterflies[END_REF][START_REF] Martin | Diversification of complex butterfly wing patterns by repeated regulatory evolution of a Wnt ligand[END_REF][START_REF] Huber | Conservatism and novelty in the genetic architecture of adaptation in Heliconius butterflies[END_REF]. Dominance at those loci is strong and several loci are found in linkage, leading to the segregation within the back-crosses of phenotypes that are either F1-like or similar to one of the parental phenotypes [START_REF] Naisbit | Mimicry: developmental genes that contribute to speciation[END_REF]. For example, the presence of a red forewing patch is associated with the B locus that is expressed in H. melpomene and hybrids carrying the dominant B allele, while no red patch is observed in H. cydno and backcrosses with bb genotype.

In contrast, the red patch of H. melpomene amaryllis, H. timareta thelxinoe and their hybrids are controlled by the same B locus, introgressed between the species [START_REF] Heliconius | Butterfly genome reveals promiscuous exchange of mimicry adaptations among species[END_REF]Pardo-Diaz C. et al., 2012), and no effects of dominance on hybrid variation is expected. Yet, our projection of hybrid phenotype in the morphospace shows that the position and the morphology of colour patches vary between species and in hybrids. This suggests that subtle variations in colour pattern involve additional small-effect quantitative loci, which may affect the accuracy of mimicry, and are expected to produce transgressive phenotypes in hybrids. Such subtle variation in the position and morphology of the colour patch is also observed in the non-mimetic pair in addition to the presence/absence of colour patches. This subtle variation has been partly attributed to epistatic interactions between heterospecific colour pattern major-effect genes [START_REF] Naisbit | Mimicry: developmental genes that contribute to speciation[END_REF] and may also involve small-effect loci. For instance, in H. melpomene, the morphology of the red forewing patch has been associated with several QTLs on different chromosomes and distinct from the major switch genes [START_REF] Baxter | Butterfly speciation and the distribution of gene effect sizes fixed during adaptation[END_REF].

The genetic basis of wing shape is completely unknown in Heliconius. Generally, shape is a complex trait involving multiple small-effect loci, e.g., in Drosophila [START_REF] Pitchers | The power of a multivariate approach to genome-wide association studies: an example with Drosophila melanogaster wing shape[END_REF].

Here, in both sets of crosses, between co-mimics or non-mimics, hybrid wing size and shape are largely intermediate between parental phenotypes. Backcross individuals are generally closer to the phenotype of the parental species representing the majority of the ancestors (BT is closer to H. t. thelxinoe than to H. melpomene amaryllis, for instance). This suggests an additive genetic basis for wing shape, possibly with complementary additive alleles underlying transgression and stronger effects of certain alleles determining the preferential directions explored in hybrids.

Transgressive effect and parental similarity

When transgression is due to the complementary action of multiple additive alleles that are present in parental lines, it is predicted that transgression is higher in phenotypically similar taxa [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF]. This is because a polygenic trait under stabilizing selection involves multiple loci with opposing effects (within species), but not necessarily the same alleles and the same directional effect in two phenotypically similar species. In hybrids between phenotypically similar parents, the different effects can thus sum for several genes, resulting in high transgression. In contrast, dissimilar parental phenotypes, evolving under directional or divergent selection, are expected to accumulate alleles driving phenotypic variation consistently in opposing directions, which then cancel each other in hybrids [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF][START_REF] Rieseberg | The genetic architecture necessary for transgressive segregation is common in both natural and domesticated populations[END_REF]Albertson & Kocher, 2005;Stelkens & Seehausen, 2009).

This prediction is rather appropriate to highly-polygenic traits involving multiple small-effect loci [START_REF] Rieseberg | Transgressive segregation, adaptation and speciation[END_REF][START_REF] Stelkens | Hybridization facilitates evolutionary rescue[END_REF], such as wing shape or size. Our results confirm that the mimetic parents are more similar in shape and size than non-mimetic species, suggesting that mimicry goes beyond colour pattern and might include shape and size, and perhaps therefore flight behaviour. Our prediction was that transgression for size and shape should be higher in crosses between co-mimics than between non-mimics. Results provided mixed support for this. On the one hand, in both pairs of species, size was generally not transgressive, and shape displayed low or null transgression when considering indexes based on the range. On the other hand, wing shape does display a higher transgression of the mean in the co-mimetic pair than in the non-mimetic pair. This could confirm the polygenic nature of wing shape and would deserve further investigation.

In contrast, for wing pattern, transgression was higher in hybrids between non-mimics than between co-mimics. This is probably not surprising given that pattern involves few largeeffect loci operating in different compartments and with dominance relationships, and hybrids between non-mimetic parents were expected generally to be very different from the parents.

Altogether, comparing the species pairs and separating the different morphological components of the visual wing phenotype, our results confirms that hybrid transgression is a complex phenomenon which is difficult to predict due to the variation in genetic architectures underlying those different traits.

Hybrid wing morphology and fitness

Wing morphology is a key component of adaptation in Heliconius butterflies because of Müllerian mimicry (Mallet & Gilbert, 1995). Mimicry relies on the fact that predators learn to associate unpalatability with a previously encountered warning signal. Thus, resemblance between chemically-defended species coexisting in the same habitat contributes to enhancing the protection conferred by the local warning signal [START_REF] Müller | Ituna and Thyridia: A remarkable case of mimicry in butterflies[END_REF][START_REF] Turner | Butterfly mimicry -the genetical evolution of an adaptation[END_REF]. Exotic or novel patterns generally suffer from higher attack rates [START_REF] Langham | Specialized avian predators repeatedly attack novel color morphs of Heliconius butterflies[END_REF][START_REF] Arias | Crossing fitness valleys: empirical estimation of a fitness landscape associated with polymorphic mimicry[END_REF].

For instance, F1 hybrids between H. melpomene rosina and H. cydno chioneus, which combine red and white on the forewing and form a novel pattern, suffer from higher predation than parental patterns, and this contributes to reproductive isolation [START_REF] Merrill | Disruptive ecological selection on a mating cue[END_REF].

However, we note that in certain cases, patterns not mimicking any known species in their local community may persist in the wild. This is the case, for instance, for the Colombian species, Heliconius heurippa, whose very prominent red and yellow forewing patch is not found in any coexisting species. This pattern is very similar to the one observed here in hybrids between H. cydno and H. melpomene, and is an example of pattern novelty produced via hybridization and introgression [START_REF] Mavarez | Speciation by hybridization in Heliconius butterflies[END_REF]Pardo-Diaz C. et al., 2012).

Because the same "toolbox" of major-effect colour pattern loci is re-used in the different Heliconius species, hybrid phenotypes may in fact match co-existing warning patterns and benefit from their protection [START_REF] Gilbert | Adaptive novelty through introgression in Heliconius wing patterns: evidence for shared genetic "tool box" from synthetic hybrid zones and a theory of diversification[END_REF][START_REF] Enciso-Romero | Evolution of novel mimicry rings facilitated by adaptive introgression in tropical butterflies[END_REF]. A fraction of the back-crosses between the non-mimetic Panamanian species present an overall pattern that roughly mimics one of the parents [START_REF] Naisbit | Mimicry: developmental genes that contribute to speciation[END_REF]. Yet, subtle variation in the morphology or position of the colour patches makes them imperfect mimics. Similarly, hybrids between the Peruvian co-mimics are very similar at first sight with both of the parents and within the overall parental range, but they slightly deviate from the parental patterns. As for wing shape and size, subtle phenotypic differences possibly translate into differences in fitness. In H. timareta thelxinoe, within-species geographic variation tend to follow the most abundant warning pattern in the local community, mimicking the abundant H. erato favorinus/H. melpomene amaryllis at mid-altitude and another species, H. telesiphe, at high altitude [START_REF] Mérot | Refining mimicry: phenotypic variation tracks the local optimum[END_REF]. H. telesiphe has narrower colour bands and a Z-shaped FW pattern which are not readily mimicked by the hybrids between H. t. thelxinoe and H. m. amaryllis, but the hybrids are more similar to H. m. amaryllis Therefore, depending on the surrounding community, hybrids are either better or worse mimics than their parents.

Depending on the extent of selection on wing morphology, possibly linked to the degree of generalization by bird predators or the local community of mimics, hybrid phenotypes may thus be disadvantageous or advantageous. If disadvantageous, it may reinforce species isolation. If neutral or advantageous, for instance if the hybrid phenotype includes discrete elements of co-occurring warning pattern or subtle improvements of mimicry accuracy, variations displayed in hybrids can spread morphological possibilities in parental species and contribute to phenotypic evolution. In Heliconius, this mechanism of novelty by hybridization is thought to be more important than novelty by mutation, since a large portion of pattern diversity observed in Heliconius (at least in the cydno-melpomene subclade) can be explained by interspecific introgression and a re-shuffling of the different loci regulating each subelements of the pattern [START_REF] Gilbert | Adaptive novelty through introgression in Heliconius wing patterns: evidence for shared genetic "tool box" from synthetic hybrid zones and a theory of diversification[END_REF][START_REF] Wallbank | Evolutionary Novelty in a Butterfly Wing Pattern through Enhancer Shuffling[END_REF][START_REF] Enciso-Romero | Evolution of novel mimicry rings facilitated by adaptive introgression in tropical butterflies[END_REF].

Nothing is known, however, about other components of wing morphology involved in mimetic accuracy, such as subtle variations of the pattern or wing shape and size, which, given the extent of interspecific genome sharing, may also have partly evolved with input from hybridization.

Conclusion

Overall, our study reveals that the occurrence and the extent of transgression depend on the genetic architecture of the trait. Completely novel phenotypes were only found for the wing colour pattern of non-mimetic hybrids, due to the shuffling of alleles at unlinked large-effect genes, switching on and off the components of colour pattern. In contrast, for traits with more likely polygenic and quantitative variation, such as wing shape or wing colour pattern in the co-mimetic pair, hybrid phenotypes were generally within the range of parental phenotypes.

However, in that case, the average hybrid morphology was distinct from a truly intermediate phenotype and transgressive along particular directions in the morphological space, perhaps revealing the underlying bias in the genetic architecture of the traits. In the mimetic species studied here, mimicry has been achieved by the introgression of at least a few large-effect loci controlling colour pattern loci. Whether more loci, affecting for instance the accuracy of mimicry, have also flowed among species and were positively selected is yet unknown. Here we focused on two pairs of butterfly species through the prism of mimicry, because selective factors associated with mimicry are well understood. Yet, any functional trait shared between species still exchanging gene flow also evolves in a multi-species context, with rare hybrids offering transgressive phenotypes and new combination of traits to selection. This raises the question to what extent occasional hybridization intervenes in evolution, by increasing standing variation, either with subtle changes or largely novel variation.

BMa

0.41** 0.45 (0.13) 0.08 (0.04) BMr 0.22** 0.22 (0.03) 0.00 (0.00) F1 TMa 0.27** 0.30 (0.07) 0.05 (0.03) F1 CMr 0.48** 0.48 (0.02) 0.00 (0.00) pattern BT 0.42** 0.46 (0.11) 0.08 (0.05) BC 0.21** 0.21 (0.04) 0.00 (0.00) BMa 0.27** 0.38 (0.12) 0.14 (0.07) BMr 0.16** 0.20 (0.05) 0.07 (0.04) F1 TMa 0.20** 0.27 (0.07) 0.09 (0.04) F1 CMr 0.22** 0.25 (0.05) 0.04 (0.02) Transgression of the mean (Tm) shape BT 0.44** 0.48 (0.11) 0.12 (0.06) BC 0.16** 0.22 (0.06) 0.07 (0.03) BMa 0.10** 0.12 (0.07) 0.00 (0.00) BMr 0.83** 0.81 (0.06) 0.43 (0.10) F1 TMa 0.05** 0.15 (0.06) 0.00 (0.00) F1 CMr 1.00

1.00 (0.00) 1.00 (0.01) pattern BT 0.08** 0.12 (0.08) 0.00 (0.00) BC 0.71** 0.73 (0.08) 0.15 (0.14) BMa 0.01** 0.02 (0.02) 0.00 (0.00) BMr 0.01** 0.01 (0.02) 0.00 (0.00) F1TMa 0.00 0.01 (0.02) 0.00 (0.00) F1 CMr 0.03** 0.05 (0.03) 0.00 (0.00)

Range Transgression (Tr)

shape BT 0.00 0.01 (0.03) 0.00 (0.00) BC 0.02** 0.04 (0.04) 0.00 (0.00) BMa 0.09** 0.06 (0.03) 0.00 (0.00) BMr 0.51** 0.34 (0.04) 0.00 (0.00) F1 TMa 0.02** 0.06 (0.02) 0.00 (0.00) F1 CMr 0.40** 0.32 (0.03) 0.00 (0.00) pattern BT 0.01* 0.03 (0.01) 0.00 (0.00) BC 0.45** 0.36 (0.01) 0.00 (0.00) BMa 0.06** 0.05 (0.03) 0.01 (0.01) BMr 0.04** 0.03 (0.01) 0.00 (0.00) F1 TMa 0.01** 0.04 (0.02) 0.00 (0.00) F1 CMr 0.02** 0.04 (0.02) 0.00 (0.00) Transgression Strength (Ts) shape BT 0.01** 0.03 (0.02) 0.01 (0.01) BC 0.06** 0.04 (0.02) 0.00 (0.00) All groups were tested together with a pairwise t-test corrected following Benjamini and Hochberg (1995).

  et al. (2009), followed by[START_REF] Husemann | Body shape differences in a pair of closely related Malawi cichlids and their hybrids: Effects of genetic variation, phenotypic plasticity, and transgressive segregation[END_REF] andHolzman and Hulsey (2017) which call this measure "distance expansion", proposed to evaluate the degree of transgression by evaluating to what extent hybrid range exceeds the parental range. In multivariate data, range was calculated as the difference between the minimal and the maximal value along each PC, and transgression strength was first evaluated on each PC separately followingStelkens et al. (2009):𝑇𝑠 (𝑃𝐶 𝑖 ) = (𝑡𝑜𝑡𝑎𝑙 𝑟𝑎𝑛𝑔𝑒 𝑎𝑙𝑜𝑛𝑔 𝑃𝐶 𝑖 -𝑝𝑎𝑟𝑒𝑛𝑡𝑎𝑙 𝑟𝑎𝑛𝑔𝑒 𝑎𝑙𝑜𝑛𝑔 𝑃𝐶 𝑖 )𝑝𝑎𝑟𝑒𝑛𝑡𝑎𝑙 𝑟𝑎𝑛𝑔𝑒 𝑎𝑙𝑜𝑛𝑔 𝑃𝐶 𝑖A global index of the strength of multidimensional transgression (Ts) was then estimated as the sum of these values, weighted by the variance explained by each PC:𝑇𝑠 = ∑ 𝑖 (𝑡𝑜𝑡𝑎𝑙 𝑟𝑎𝑛𝑔𝑒 𝑎𝑙𝑜𝑛𝑔 𝑃𝐶 𝑖 -𝑝𝑎𝑟𝑒𝑛𝑡𝑎𝑙 𝑟𝑎𝑛𝑔𝑒 𝑎𝑙𝑜𝑛𝑔 𝑃𝐶 𝑖 )𝑝𝑎𝑟𝑒𝑛𝑡𝑎𝑙 𝑟𝑎𝑛𝑔𝑒 𝑎𝑙𝑜𝑛𝑔 𝑃𝐶 𝑖 * 𝑉𝑎𝑟 (𝑃𝐶 𝑖 )

H

  . c. chioneus; BMr, backcross with H. m. rosina. (B-D) Illustrations of alternative measures of hybrid transgression. (B) Transgression of the mean, measured as the departure of the mean hybrid phenotype from that of the two parental phenotypes, relative to the distance between the parents (based on Renaud et al. (2012); Renaud et al. (2017)). (C) Range transgression, in which transgressive hybrids are those for which the distance from the parental phenotypes exceeds the parental range (based on Parsons et al. (2011)). (D) Transgression strength along principal component axes, measured as the extent to which hybrid range extends beyond the range of parental phenotypes (based on Stelkens et al. (2009) and Husemann et al. (2017)).

Figure 2 .

 2 Figure 2. Wing size distribution of hybrids and parents.

Figure 3 .

 3 Figure 3. Variation of hybrid wing shape and pattern.

Figure S3 .

 S3 Figure S3. PCA on forewing and hindwing shape PCA on wing shape for parents and hybrids between H. cydno chioneus and H. melpomene rosina on forewing pattern (A, C) and the hindwing pattern (B, D), between H. timareta thelxinoe and H. melpomene on forewing pattern (E, G) and the hindwing pattern (F, H). Open circles represent females and filled dots represent males. Shape variation is illustrated next to each axis, where broken lines represent minimum negative values of the axis, and full lines represent maximum values.

Table 1 . Transgression in hybrids.

 1 Each index (left column) is calculated for wing shape and wing pattern for the F1 and backcrosses from the non-mimetic (H. t. thelxinoe [T] and H. m. amaryllis [Ma]) and mimetic (H. c. chioneus [C] and H. m. rosina [Mr]) crosses. Indices based on 30 randomly-chosen or simulated hybrids were based on 1000 bootstraps. Stars indicate the probability that the observed values of the transgression indices were drawn from the simulated distribution of theoretically non-transgressive hybrids; ** p < 0.001 and * p < 0.01.

	Index	Mean	Mean	Index	Mean	Mean
	based	index (sd),	index (sd),	based	index (sd),	index (sd),
	on all samples	30 randomly-chosen	30 simulated hybrids	on all samples	30 randomly-chosen	30 simulated hybrids
		samples			samples

Figure 1. Study species and methods.

  Phylogenetic relationship between the parental species, map of repartition, parental phenotypes and design of the crosses. Between Peruvian mimetic parents: F1 T/Ma , first generation hybrids between H. t. thelxinoe and H. m. amaryllis; BT, backcross with H. t. thelxinoe; BMa, backcross with H. m. amaryllis. Between Panamanian non-mimetic parents: F1 C/Mr , first generation hybrids between H. c. chioneus and H. m. rosina; BC, backcross with

(A)

Table S1 . Shape and pattern variation between parental species

 S1

		ANOVA/MANOVA	Mean of individual		Mahalanobis
				distances		distances
							co-	non-
		co-mimics	non-mimics	co-mimics non-mimics		mimics	mimics
		Ma-T	Mr-C	Ma-T	Mr-C	ratio	Ma-T	Mr-C
		F1,156=41,	F1,129=90,
	size	P<0.001	P<0.001	0.8	1.05	1.3
		F14,143=32,	F14,116=62,
		Pillai=0.76	Pillai=0.88
	shape	P<0.001	P<0.001	0.063	0.076	1.2	11.1	28.8
		F10,139=80,	F10,99=3734,
		Pillai=0.85	Pillai=0.99
	pattern	P<0.001	P<0.001	71.9	224.1	3.1	28.7	437

Table S2 . Shape and pattern variation between hybrids and parents LDA

 S2 (CV) indicates the cross-validation rate when doing a linear discriminant analysis between each hybrid group and each parents

				SHAPE						PATTERN
			LDA (CV)	MANOVA			LDA (CV)	MANOVA
			Ma T	Ma	T			Ma	T	Ma	T
				F1,20=16.5;	F1,20=5.4;					F1,20=10.02;	F1,20=6.4;
		BT	93% 77%	p<0.001	p<0.001		BT	87%	81%	p<0.001	p<0.001
	FW	BMa 79% 95%	F1,20=9.0; p<0.001	F1,20=26.8; p<0.001	FW	BMa	77%	98%	F1,20=6.8; p<0.001	F1,20=55.8; p<0.001
				F1,20=11.4;	F1,20=7.1;					F1,20=6.5;	F1,20=10.1;
		F1T/Ma 93% 83%	p<0.001	p<0.001		F1T/Ma 82%	85%	p<0.001	p<0.001
			Mr C	Mr	C			Mr	C	Mr	C
				F1,20=19.7;	F1,20=7.7;					F1,20=257.6;	F1,20=17.7;
		BC	89% 78%	p<0.001	p<0.001		BC	98%	91%	p<0.001	p<0.001
	FW	BMr 87% 89%	F1,20=11.3; p<0.001	F1,20=17.8; p<0.001	FW	BMr	88% 100%	F1,20=14.4; p<0.001	F1,20=2071; p<0.001
				F1,20=7.7;	F1,20=6.5;					F1,20=248.9;	F1,20=1134;
		F1C/Mr 84% 85%	p<0.001	p<0.001		F1C/Mr 99% 100%	p<0.001	p<0.001
			Ma T	Ma	T			Ma	T	Ma	T
				F1,20=11.6;	F1,20=4.5;					F1,20=104;	F1,20=9.9;
		BT	88% 75%	p<0.001	p<0.001		BT	84%	87%	p<0.001	p<0.001
	HW	BMa 74% 93%	F1,20=6.6; p<0.001	F1,20=22.1; p<0.001	HW	BMa	71%	95%	F1,20=4.3; p<0.001	F1,20=37.5; p<0.001
				F1,20=6.2;	F1,20=5.9;					F1,20=3.8;	F1,20=7.3;
		F1T/Ma 79% 79%	p<0.001	p<0.001		F1T/Ma 75%	81%	p<0.001	p<0.001
			Mr C	Mr	C			Mr	C	Mr	C
				F1,20=42.2;	F1,20=9.1;					F1,20=3994;	F1,20=10.1;
		BC	96% 82%	p<0.001	p<0.001		BC	100%	83%	p<0.001	p<0.001
	HW	BMr 82% 93%	F1,20=12.3; p<0.001	F1,20=24.5; p<0.001	HW	BMr	90% 100%	F1,20=23.5; p<0.001	F1,20=279; p<0.001
				F1,20=7.5;	F1,20=14.2;					F1,20=888;	F1,20=22.4;
		F1C/Mr 86% 91%	p<0.001	p<0.001		F1C/Mr 100%	90%	p<0.001	p<0.001

Journal of Evolutionary Biology

Supplementary materials

Hybridization and transgressive exploration of wing morphology in Heliconius butterflies Table S3. Transgression strength along each PC direction %var stands for the variation explained by the given PC, Ts Pci is the difference between total range and parental range, expressed as a fraction of parental range along each PC. Ts is the global index and Max the maximum Ts Pci. Table S4. Mean transgression controlled by the ratio of parental distances Since the index Tm is standardized by the distance between parents, this might artificially bias index values towards higher values when species are more similar. For instance, we noted that, in the null model, simulated hybrids have a higher transgression of the mean in the co-mimetic pair than in the non-mimetic pair. Therefore, we also compared the transgression value when correcting from the ratio of parental distances (Table S1).


Ivan V Zhukov

Alexey S Kiryutin

Ziqing Wang

Milan Zachrdla

Alexandra V Yurkovskaya

Konstantin L Ivanov
email: ivanov@tomo.nsc.ru

Fabien Ferrage
email: fabien.ferrage@ens.psl.eu

Surprising absence of strong homonuclear coupling at low magnetic field explored by two-field NMR spectroscopy

Strong coupling of nuclear spins, which is achieved when their scalar coupling is greater than or comparable to the difference in their Larmor precession frequencies in an external magnetic field, gives rise to efficient coherent longitudinal polarization transfer. The strong-coupling regime can be achieved when the external magnetic field is sufficiently low, as is reduced proportional to the field strength. In the present work, however, we demonstrate that in heteronuclear spin systems these simple arguments may not hold, since heteronuclear spin-spin interactions alter the value. The experimental method that we use is two-field NMR (Nuclear Magnetic Resonance), exploiting sample shuttling between a high field, at which NMR spectra are acquired, and low field, where strong couplings are expected, at which NMR pulses can be applied to affect the spin dynamics. By using this technique, we generate zero-quantum spin coherences by means of non-adiabatic passage through a level anticrossing and study their evolution at low field. Such zero-quantum coherences mediate the polarization transfer under strong coupling conditions. Experiments performed with an 13 C labelled amino acid clearly show that the coherent polarization transfer at low field is pronounced in the 13 C-spin subsystem under proton decoupling. However, in the absence of proton decoupling, polarization transfer by coherent processes is dramatically reduced, demonstrating that heteronuclear spin-spin interactions suppress the strong coupling regime even when the external field is low. A theoretical model is presented, which can model the reported experimental results.

I. Introduction

The topological and conformational information provided by scalar couplings lies at the foundation of the analytical power of NMR spectroscopy [START_REF] Ernst | Principles of Nuclear Magnetic Resonance in One and Two Dimensions[END_REF][START_REF] Keeler | Understanding NMR Spectroscopy[END_REF][START_REF] Levitt | Spin Dynamics : Basics of Nuclear Magnetic Resonance-2nd Ed[END_REF][START_REF] Cavanagh | Protein NMR Spectroscopy : Principles and Practice[END_REF]. The strong coupling case is encountered when scalar coupling constants are not negligible with respect to the difference of resonance frequency between the coupled spins [START_REF] Keeler | Understanding NMR Spectroscopy[END_REF]. Understanding strong scalar couplings and their spectral signature was essential when NMR was introduced for chemical analysis, which was typically performed at magnetic fields considered today as low [START_REF] Bodenhausen | Proton-Coupled Carbon-13 J Spectra in the Presence of Strong Coupling[END_REF][START_REF] Pfändler | Strong Coupling Effects in Z-Filtered Two-Dimensional NMR Correlation Spectra[END_REF]. Modern high-field NMR is widely based on the exploitation of weak scalar couplings, so that strong scalar couplings have remained a nuisance, in particular in aromatic spin systems (Vallurupalli et al., 2007;[START_REF] Foroozandeh | Ultrahigh-Resolution NMR Spectroscopy[END_REF]. Recently, the development and availability of benchtop NMR spectrometers operating at low or moderate magnetic fields [START_REF] Grootveld | Progress in Low-Field Benchtop NMR Spectroscopy in Chemical and Biochemical Analysis[END_REF], has revived the interest in the understanding of strong scalar couplings in conventional NMR.

Contrarily to conventional NMR, NMR at near-zero or ultralow magnetic fields (ZULF-NMR), explores the benefits of NMR in the strong scalar-coupling regime. At such magnetic fields, typically smaller than 1 T, scalar coupling interactions dominate all Zeeman interaction and dictate the eigenstates of spin systems and transition energies obtained in spectra [START_REF] Ledbetter | Near-Zero-Field Nuclear Magnetic Resonance[END_REF][START_REF] Tayler | Invited Review Article: Instrumentation for Nuclear Magnetic Resonance in Zero and Ultralow Magnetic Field[END_REF][START_REF] Blanchard | Zero-to Ultralow-Field NMR[END_REF]. However, for homonuclear couplings, the transition between the weak-and strongcoupling regimes occurs in a range of magnetic fields, where the Zeeman interaction is still dominant [START_REF] Miesel | Coherence Transfer During Field-Cycling NMR Experiments[END_REF][START_REF] Ivanov | Coherent Transfer of Hyperpolarization in Coupled Spin Systems at Variable Magnetic Field[END_REF][START_REF] Ivanov | The Role of Level Anti-Crossings in Nuclear Spin Hyperpolarization[END_REF][START_REF] Appelt | Paths from Weak to Strong Coupling in NMR[END_REF]Türschmann et al., 2014). This transition between weak and strong couplings can be investigated by varying the magnetic field applied to the sample on a high-field magnet, which is usually performed by moving the sample through the stray field with a shuttle system (Roberts and Redfield, 2004a, b;Redfield, 2012;[START_REF] Wagner | High-Resolution Magnetic Relaxation Dispersion Measurements of Solute Spin Probes Using a Dual-Magnet System[END_REF][START_REF] Bryant | Nuclear Magnetic Resonance and Spin Relaxation in Biological Systems[END_REF][START_REF] Goddard | The Magnetic Field and Temperature Dependences of Proton Spin-Lattice Relaxation in Proteins[END_REF]Chou et al., 2016;Chou et al., 2017;Charlier et al., 2013;Cousin et al., 2016a;Cousin et al., 2016b;[START_REF] Zhukov | Field-Cycling NMR Experiments in Ultra-Wide Magnetic Field Range: Relaxation and Coherent Polarization Transfer[END_REF][START_REF] Kiryutin | A Fast Field-Cycling Device for High-Resolution NMR: Design and Application to Spin Relaxation and Hyperpolarization Experiments[END_REF]. These studies have highlighted the effects of level anti-crossings (LACs) [START_REF] Miesel | Coherence Transfer During Field-Cycling NMR Experiments[END_REF][START_REF] Ivanov | The Role of Level Anti-Crossings in Nuclear Spin Hyperpolarization[END_REF]. When the passage through a LAC is slow, the transition is adiabatic and the population of eigenstates is smoothly converted to the new eigenstates. When the transition is fast, coherences can be generated between the new eigenstates and time-oscillations of the population of high-field eigenstates can be observed (Pravdivtsev et al., 2013;Kiryutin et al., 2013). As usual, non-adiabatic variation, which gives rise to excitation of coherences, means that the adiabatic eigenstates of the spin system change with time fast as compared to the rate of internal evolution of the system. Specifically, for each pair of adiabatic states, and , the parameter is much greater than 1 (here is the energy difference between the states, measured in angular frequency units). When , switching is adiabatic and populations follows the time-dependent eigenstates. This phenomenon has been observed on a variety of homonuclear spin systems. Heteronuclear scalar couplings have been shown to alter LACs in homonuclear spin systems [START_REF] Korchak | High Resolution NMR Study of T 1 Magnetic Relaxation Dispersion. III. Influence of Spin 1/2 Hetero-Nuclei on Spin Relaxation and Polarization Transfer among Strongly Coupled Protons[END_REF]; yet, the properties of such heteronuclear couplings on LACs are not fully understood, in particular, in spin systems with extensive networks of homo-and heteronuclear scalar couplings.

Here, we investigate the effect of heteronuclear scalar couplings on LACs in a spin system typical of biomolecular NMR, a uniformly carbon-13 labeled amino acid (leucine), which combines extensive networks of homo-and heteronuclear scalar couplings. Essentially, we exploit the ability to apply composite pulse decoupling on our two-field NMR spectrometer (Cousin et al., 2016a) to switch on and off heteronuclear scalar couplings at low magnetic field. We demonstrate that heteronuclear scalar couplings alter LACs by sustaining the weak-coupling regime in a carbon-13 homonuclear spin system. Composite pulse decoupling at low magnetic field restores the strong scalar coupling regime in the carbon-13 nuclei of the isopropyl group of leucine at 0.33 T. Our results identify how heteronuclear couplings alter homonuclear couplings at low magnetic fields, which could be exploited in low-field NMR methodology and may be considered in further developments of total correlation spectroscopy (TOCSY) [START_REF] Braunschweiler | Coherence Transfer by Isotropic Mixing -Application to Proton Correlation Spectroscopy[END_REF] mixing sequences in high-field NMR.

II. Methods

A. Sample preparation

Experiments have been performed using the following sample: 76 mM 99% enriched 13 C, 15 N labeled Lleucine (Leu) in 90% H 2 O 10% D 2 O solution. 13 C, 15 N enriched L-leucine were purchased from Sigma-Aldrich and used as it stands. 13 C-NMR spectrum of the labelled Leu molecule is shown in Figure 1. We also show separately the signals of the individual carbon nuclei. Broadband proton decoupling was used to simplify the spectrum. Here, we will focus on a three-spin system, formed by the C γ and two C δ nuclei of the isopropyl moiety. We will study polarization transfer in this subsystem upon fast switch of the external magnetic field obtained by a transfer of the sample though the stray field of a high-field NMR magnet.

B. Field-cycling NMR experiments

NMR experiments were performed on a two-field NMR spectrometer (Cousin et al., 2016a) with fast sample shuttling (Charlier et al., 2013). The high field 14.1 T is the detection field of a 600 MHz NMR spectrometer while the low field is 0.33 T corresponding to 14 MHz 1 H Larmor frequency. The magnetic field in the low-field centre is sufficiently homogeneous (inhomogeneities of the order of 10 ppm) so that radiofrequency (RF) pulses can be applied by using a triple-resonance NMR probe, as described previously (Cousin et al., 2016a).

Field-cycling NMR experiments were run according to the pulse sequences depicted in Figure 2. First, a non-equilibrium state is generated at by applying a selective  pulse to the C δ2 nucleus (shaped RE-BURP pulse (Geen and Freeman, 1991), the pulse duration was 46.4 ms, the peak RF-field amplitude was adjusted to cover ca. 100 Hz bandwidth around the center of C δ2 signal). RE-BURP shaped pulse was used since it is less sensitive to the initial nuclear magnetization state than I-BURP (Geen and Freeman, 1991) and has narrow excitation profile. To improve the selectivity of the pulse, simultaneous proton decoupling was used, which reduces multiplet overlap in the carbon-13 NMR spectrum. Following this preparation, the sample was shuttled from the high-field center to the low-field center with a duration 107 ms. The field jump is fast enough to be non-adiabatic and it is aimed to excite a spin coherence. Subsequently, the coherence evolves at during a variable time period . The shuttle transfer back to the high-field center leads to a second field jump with a duration 94 ms. This second non-adiabatic field jump to converts the coherence into a population difference. Detection is performed after a /2 pulse on the carbon-13 channel in the presence of proton decoupling. We perform two types of experiments, in which the carbon spin coherence (zero-quantum coherence, ZQC) evolves at in the absence (see Figure 2A) and in the presence (see Figure 2B) of proton composite-pulse decoupling. Decoupling at has been performed using composite pulse decoupling pulse with the WALTZ decoupling with MLEV-64 supercycle (Shaka et al., 1983;[START_REF] Levitt | Supercycles for Broad-Band Heteronuclear Decoupling[END_REF] at low field on the proton RF-channel (operating at 14 MHz corresponding to the proton NMR frequency at 0.33 T). Composite pulse decoupling is used because of the rather high field inhomogeneity at 0.33 T, which is of the order of 10 ppm: under such conditions continuous-wave decoupling would require more power, potentially giving rise to sample heating. The -dependence of polarization is expected to be oscillatory, due to the coherent polarization exchange within the expectedly strongly coupled system of the C γ and two C δ carbon-13 nuclei.

III. Theory

A. Polarization transfer in a 3-spin system

In this subsection, we provide a theoretical description of the field-cycling NMR experiments. First, we present the analytical treatment of polarization transfer among two nuclei of the same kind, here spin and spin (e.g. two carbon-13 nuclei), in the presence of a third spin , which can be a heteronucleus (e.g. here a proton). This is the minimal system allowing us to detail the effect of a heteronucleus on polarization transfer among strongly coupled spins. We assume that spins and are in strong coupling conditions, meaning that the difference, , in their Zeeman interaction frequencies with the external field is smaller than or comparable to the scalar-coupling constant, , between them. When the strong coupling regime is achieved, the zero-quantum part of the scalar coupling, given by the operator , becomes active, giving rise to flips and flops of spins and . The couplings to the third spin S, and , are assumed to be unequal (otherwise coupling to the proton give rise to an identical shift of the NMR frequencies of spins 1 and 2 and does not modify the eigenstates of this subsystem). The Hamiltonian of the spin system can be written as follows (in units):

(1)

Here , and are the spin operators;

, and stand for the NMR frequencies of the corresponding nuclei (with being the corresponding gyromagnetic ratios and being the chemical shifts). We assume that the heteronucleus is coupled weakly to I spins due to the large difference in their NMR frequencies, i.e., , and keep only the secular part of the heteronuclear coupling Hamiltonian.

In the present case, the nuclear magnetic number, , of spin is a "good quantum number", which is conserved because commutes with the Hamiltonian. For this reason, it is possible to find the solution for the spin dynamics of spins and for two separate cases, which corresponds to the two different values of being and , i.e., spin is in the "spin-up" state or "spin-down" state. In each case, the Hamiltonian of the carbon subsystem is as follows:

(2)

Hence, in the Hamiltonian given by eq. ( 1) we replace the operator by the value, which is .

Hence, the value is modified and it depends on the value:

(3)

The eigenstates of the subsystem of spin 1 and spin 2 are (4)

Here the "mixing angle" is given by the values of and :

. When approaches zero, the mixing angle goes to meaning that the eigenstates become singlet and triplet states: the spins are strongly coupled. When is much greater than the coupling, the eigenstates are obviously the Zeeman states.

Even in this simple system, it is clear that the condition is not sufficient to guarantee strong coupling of the two carbons. Indeed, when is greater than and the carbon spins become weakly coupled in the two sub-ensembles, corresponding to .

How do heteronuclear couplings affect polarization transfer in the carbon system? We assume that at one of the spins has polarization and the other spin is not polarized, 0. Hereafter, it is convenient to use normalization 1. The state of the spin system is then given by the density operator (5) As shown previously [START_REF] Miesel | Coherence Transfer During Field-Cycling NMR Experiments[END_REF], in the two-spin system of and , in the absence of coupling to any other spin, the polarization evolves with time as follows:

(6) where and the oscillation frequency is the frequency of the ZQC between the eigenstates and . Hence, coherent exchange of polarization is taking place. As becomes smaller the frequency of the oscillations decreases, but the amplitude increases: at we obtain and complete exchange is possible when . In the presence of scalar couplings to the third spin , here a proton ( and are carbon-13 nuclei), the expressions should be modified: the evolution should be calculated for each specific spin state of the proton, and , and sum of the two curves should be taken. We obtain at the following expression:

(7)

where the evolution frequencies are equal to .

The time dependence of the expectation value for the longitudinal polarizations of spins and is presented in Figure 3 in the presence and the absence of scalar couplings to a heteronucleus. In the absence of heteronuclear coupling the two strongly coupled spins (the strong coupling condition is fulfilled since ) almost completely exchange polarizations. The polarization transfer is of a coherent nature and the frequency of the oscillations is close to the scalar-coupling constant . In the presence of different heteronuclear scalar couplings to the third spin , the time-evolution changes considerably. The two spins are no longer in the regime of strong coupling, since . The efficiency of polarization transfer is reduced and complete exchange of polarization is no longer possible. The time dependence also becomes more complex: instead of a single frequency found in the previous case, here two frequencies appear: and . Hence, when couplings to heteronuclei are present, the condition does not guarantee that the homonuclei are in the strong-coupling regime. These results show that the interaction with a heteronucleus clearly alters polarization transfer in strongly coupled networks. Consequently, we expect strong effects of heteronuclear interactions on polarization transfers in systems with several heteronuclei. Notably, we anticipate that polarization transfer among strongly coupled carbon spins will be dramatically different in the presence of proton decoupling, which effectively removes proton-carbon spin-spin interactions.

B. Spin dynamics simulations

In addition to this simple model, we carried out numerical simulations in a realistic multi-spin system: the isopropyl group of carbon-13 labeled leucine. This spin system contains three carbon-13 nuclei , , and : the C γ carbon-13 and the two C  carbon-13 nuclei. In addition, the spin system includes seven protons : each C  nucleus is coupled to the three protons of the methyl group, and the C γ carbon-13 nucleus is coupled to one proton. We model the effects of fast field variation and coherent spin dynamics at low field. We consider two cases, namely, polarization transfer in the presence and in the absence of proton decoupling.

The simulation method is as follows. The band-selective inversion pulse on spin generates the initial density operator for the three-spin system:

(8)

Hence, we generate a population difference for the states , and : the first state is overpopulated, while the other two states are underpopulated. The three-spin system under study, C γ , C δ1 and C δ2 , has a LAC at 1.1 T, see Figure 4. Upon passage through a LAC during the field jump due to the sample shuttle transfer, the population difference is expected to be converted into a coherence between the states, which have the LAC: these adiabatic states correspond to the and states at high fields. To calculate the actual spin state at we solve numerically the Liouville-von Neumann equation for the spin density operator (9)

The Hamiltonian of the spin system at a magnetic field is as follows:

(10) Here and are the carbon and proton gyromagnetic ratios, and are the chemical shifts of the -th carbon and -th proton, is the scalar coupling constant between the -th and -th carbon, is the scalar coupling constant between the -th and -th proton, is the scalar coupling constant between the -th carbon and -th proton, and are the spin operator of the -th carbon and -th proton. Given the range of magnetic fields considered here, heteronuclear scalar couplings are considered to be weak. The precise values of the calculation parameters are given in Table 1. Since the magnetic field changes with time, the Hamiltonian is also time-dependent. In the calculation we consider three carbons and seven protons (six protons in three CH 3 -groups and the -proton). Using this Hamiltonian we evaluate the density operator after the first field jump, . The Liouville-von Neumann equation is integrated using 1 ms time increments and assuming that for each step the Hamiltonian is constant, similarly to simulations carried out for relaxation experiments [START_REF] Bolik-Coulon | Theoretical and Computational Framework for the Analysis of the Relaxation Properties of Arbitrary Spin Systems. Application to High-Resolution Relaxometry[END_REF]. In the calculation, we ignore relaxation effects, since the dimensionality of the relaxation superoperator is too big for the multi-spin system considered here and our focus is on coherent effects.

At

the density operator evolves under a constant Hamiltonian, at the end of the evolution period it becomes as follows:

(11)

The field jump is simulated numerically in the same way as the first field jump (the time interval is split into many small steps). Finally, knowing the density operator at , we evaluate the NMR signals of the nuclei of interest as the expectation values of their -magnetization .

The method used for modelling the experiments with decoupling at is different. After evaluating the density operator we trace out the proton degree of freedom and define the density operator of the carbon subsystem as , with the argument that proton polarization is destroyed by decoupling. The partial trace procedure implies that when is a protoncarbon density operator (in the notation of spin states stand for the proton states and stand for the carbon states), the elements of the carbon density operator are:

. One should note that proton two-spin operators may contain a zero-quantum component, which would withstand proton decoupling. Consideration of effects of such coherences is beyond the scope of this work: we expect this to only lead to small perturbations of the observed behavior. Then we introduce the Hamiltonian of the carbon subsystem (12) Using this Hamiltonian, we evaluate the density operator of the 13 C spins at the end of the evolution period as follows (13) The final step in evaluating the ZQC evolution is introducing the carbon-proton density operator. This is done by multiplying and the density operator of non-polarized protons (as decoupling removes any proton spin order). Hence (14) where is a unity matrix. The final step of the calculation, the field jump , is modelled in the same way as in the previous case.

Finally, we would like to comment on the dependence, which was used in calculation. The distance dependence of the magnetic field is precisely known but the precise is not known. We modelled this dependence assuming that motion goes with a constant speed (in experiments, constantspeed motion is achieved after a 5-10 ms lag delay for acceleration). Non-ideal agreement between theory and experiment can be attributed to the fact that the precise dependence is not known (our previous works (Pravdivtsev et al., 2013;Kiryutin et al., 2013) show that the knowledge of is required for modeling).

IV. Results and discussion

The experimental τ-dependences of the measured spin polarization are shown in Figure 5. One can see that without decoupling no coherent behavior is found: polarization simply decays due to relaxation and no coherent oscillations are visible (Figure 5A). In the presence of proton decoupling the situation is drastically different: coherent oscillations are clearly observed, which mediate polarization exchange between the C δ1 and C δ2 nuclei. We attribute such polarization exchange to the ZQC, which is generated by passage through the LAC. The coherence gives rise to exchange of the populations of the two states, which experience the LAC. These levels are correlated with the and high-field states. Hence, polarization transfer gives rise to population exchange of the states (initially overpopulated state) and

(initially underpopulated state). As a result, the state of the first spin, C γ , does not change, but the other two spins, C δ1 and C δ2 , exchange polarizations. With the available speed and range of the field-cycling, other coherences are not excited, i.e., non-adiabatic variation of the Hamiltonian is achieved only for the pairs of levels that have the LAC in between and , i.e., only the LAC shown in Figure 4 contributes to spin mixing. The C γ spin never shows any oscillatory polarization transfer, which is an indication that the specific LAC is responsible for the observed effect. In conclusion, a zero-quantum coherence of the two carbon-13 nuclei C δ1 and C δ2 is excited by fast magnetic field jump between 14.1 T and 0.33 T.

The oscillatory behavior does not show up in the absence of proton decoupling. There are two reasons for that. First, the multiple proton-carbon couplings give rise to a set of ZQC frequencies, instead of a unique frequency in the presence of decoupling. Second and more importantly, proton-carbon-13 couplings prevent the carbon subsystem from reaching the strong-coupling regime. Thus, the amplitude of coherent evolutions is drastically reduced (see Eq. 7) and becomes negligible (Fig. 5.A). As a result, in experiments without decoupling the ZQC decays because of inhomogeneous broadening of the ZQC evolution frequency, i.e. relaxation. We would like to stress that the ZQC of interest is excited by the field jump, which is identical for experiments with and without proton decoupling at low field. However, the ZQC does not reveal itself and does not give rise to efficient polarization transfer in the experiment without decoupling.

These considerations are confirmed by theoretical modeling (Figure 6). In the presence of carbonproton couplings coherent oscillations are hardly observed: only fast oscillations of very small amplitude can be seen in the simulated curves. By contrast, in the absence of the proton-carbon couplings, i.e., when decoupling is used, coherent evolutions become manifest with slower oscillations of larger amplitude. The results of numerical modeling are in good agreement with the experimental data. As relaxation effects are not taken into account in simulations, to ease comparison we subtracted the slowly relaxing background from the experimental time traces. In addition, we rescaled all calculated traces with the same factor; then the starting polarization values were adjusted individually to achieve the best agreement with the experimental data. Such a data treatment becomes necessary because relaxation is active not only during spin mixing at the field, but also during the field jumps. The agreement between the experimental data and simulation in Figure 6 is not ideal, possibly because some small long-range scalar couplings are not included in the simulation but most likely because the field switching profile is not known exactly: previous studies of the spin dynamics in field-cycling NMR experiments (Pravdivtsev et al., 2013;Kiryutin et al., 2013) suggest that using the precise profile is crucial for simulating coherent polarization transfer phenomena. The absence of strong-coupling regime, in spite of scalar coupling constants larger than the difference in Larmor frequencies is somewhat counterintuitive but clearly explained when taking into account the effect of large heteronuclear scalar couplings (Eqs. 2-4). In the present case, the effect is even more pronounced since the two  carbon-13 nuclei of leucine are coupled to no less than 3 protons each, further splitting resonance frequencies in the absence of proton decoupling. A conventional way to present the weak coupling regime consists in stating that the part of the scalar coupling Hamiltonian (Eq. 1) that is proportional to a zero-quantum product operator is non-secular in the frame of the Zeeman interactions of the two coupled spins, which is true if the scalar coupling constant is much smaller than the difference in Larmor frequencies of the two spins. Here, the weak-coupling regime is extended because this zero-quantum part can be considered non-secular in the interaction frame of the heteronuclear scalar couplings (note that the perturbative treatment is allowed to the extent that the heteronuclear coupling constants are much larger than the homonuclear coupling).

A particular consequence of the observation we report here can be relevant for experiments where the strong scalar-coupling regime is created by radio-frequency irradiation: isotropic mixing for TOCSY [START_REF] Braunschweiler | Coherence Transfer by Isotropic Mixing -Application to Proton Correlation Spectroscopy[END_REF]. We have recently introduced a two-field TOCSY experiment where isotropic mixing is carried out at 0.33 T and chemical shift evolutions occur at high field [START_REF] Kadeřávek | Full Correlations across Broad NMR Spectra by Two-Field Total Correlation Spectroscopy[END_REF], which makes broadband carbon-13 TOCSY straightforward. This study included a control experiment where no radio-frequency pulses were applied at low field (see Figure 3.b in reference [START_REF] Kadeřávek | Full Correlations across Broad NMR Spectra by Two-Field Total Correlation Spectroscopy[END_REF]. Intuitively, one would have expected cross-peaks to be observed for carbon-13 nuclei in strongly-coupled networks at 0.33 T. Some cross peaks could indeed be observed within the aliphatic carbon region of leucine and in the aromatic ring of phenylalanine. The current investigation suggests that strong scalar couplings between carbon-13 nuclei are less prevalent than expected at 0.33 T. The observed cross-peaks were possibly due to cross relaxation and not necessarily coherent evolution under strong scalar couplings. Conventional TOCSY experiments might also be altered by the effect of large heteronuclear scalar couplings. In this case, isotropic mixing sequences have been optimized on isolated pairs of two coupled spins [START_REF] Kadkhodaie | Broad-Band Homonuclear Cross Polarization Using Flip-Flop Spectroscopy[END_REF], excluding the effects of scalar couplings to heteronuclei or as heteronuclear decoupling sequences that happen to be efficient at isotropic mixing [START_REF] Rucker | Broad-Band Homonuclear Cross Polarization in 2D NMR Using DIPSI-2[END_REF][START_REF] Shaka | Iterative Schemes for Bilinear Operators -Application to Spin Decoupling[END_REF]. Although isotropic mixing sequences decouple heteronuclear scalar couplings, optimizing simultaneously for homo-and heteronuclear scalar coupling operators may improve homonuclear coherence transfers. Such effects of couplings to heteronuclei are of relevance for abundant nuclei such as protons or 13C spins in uniformly 13C-labelled molecules.

V. Conclusions

In this work, we present a study of coherent polarization transfer in a system of (strongly) coupled 13 C nuclei. Spin coherences are zero-quantum coherences, which are generated by a fast non-adiabatic magnetic field jump. Such coherences are excited most efficiently when the system goes through a LAC during the field switch. Here we indeed pass through a LAC in a system of three coupled 13 C spins and investigate the spin dynamics at low fields, where strong couplings of the carbon spins are expected.

We can clearly demonstrate that the polarization transfer in the carbon-13 spin subsystem is strongly affected by spin-spin interactions with the protons in the molecule. In this situation, the role of these interactions can be determined by comparing the experiments with and without proton decoupling at low fields. When decoupling is used, we observe coherent polarization exchange between two of the three carbons: such a behavior is typical when the spin coherences are excited upon non-adiabatic passage through a specific LAC. In the absence of decoupling, i.e., when heteronuclear interactions are present, we cannot observe such a behavior: polarization transfer is very inefficient and coherent phenomena are not found. We attribute this to the fact that relatively strong proton-carbon couplings (i) drive the carbon system away from the strong coupling condition and (ii) give rise to a set of evolution frequencies instead of a unique ZQC frequency. These considerations are supported by an analytical model of a three-spin system and numerical simulations in a multi-spin system.

Our results are of importance for analyzing polarization transfer phenomena at low magnetic fields and for interpreting NMR data obtained under apparent strong coupling conditions. Under such conditions heteronuclear spin-spin interactions might disturb "strong coupling" of homonuclei and substantially alter spin dynamics. Similar effects also often arise in dynamic nuclear polarization, where the difference in the electron-nuclear couplings for nuclei located at different distances from the electron hampers nuclear spin diffusion, giving rise to the spin diffusion barrier around the electron spin [START_REF] Ramanathan | Dynamic Nuclear Polarization and Spin Diffusion in Nonconducting Solids[END_REF].

Figure

  Figure 1. Structure of 13 C, 15 N L-leucine (A) and 150.9 MHz 13 C-NMR spectrum (B) under broadband 1 H decoupling. Signal of each carbon nuclei is also shown separately (C)-(F). The multiplet structure in the spectrum is due to 13 C-

Figure 2 .

 2 Figure 2. Experimental protocols of field-cycling NMR experiments without 1 H decoupling at the low field (A) and

Figure 3 .

 3 Figure 3. Polarization transfer among two strongly coupled nuclei (A) in the absence and (B) in the presence (bottom) of a heteronucleus. Here, we present the time dependence of (black solid lines) and (red dashed lines), normalized to the initial value of . The density operator at time 0 is . Parameters of the simulation were 10 Hz, 30 Hz, and (A) 0 Hz and (B) 100 Hz.

Figure 4 .

 4 Figure 4. (A) Energy levels of the spin system at variable magnetic field strength in the absence of scalar coupling with protons. Levels are assigned at high field, where the spin system is weakly coupled. (B) Energy levels, corresponding to the ααβ and αβα states at high field, have a LAC at 1.1 T, which is responsible for generation of the zero quantum coherences. To visualize the energy levels better, in the right panels we have subtracted the large Zeeman energy from the actual energy and show the energy difference. The calculation is done using parameters listed inTable 1 and neglecting carbon-proton couplings.

Figure 5 .

 5 Figure 5. Observed -dependence of the polarizations of carbon-13 nuclei C γ , C δ and C δ measured (A) without 1 H

Figure 6 .

 6 Figure 6. Calculated -dependence of polarization (lines) overlaid with the observed time traces (points) obtained (A) without 1 H decoupling, and (B) with 1 H decoupling. The slowly relaxing background (compare with the data shown in Figure 5) has been subtracted from the time traces, to enable comparison between theory and simulations. Observed NMR intensities are normalized to intensities in 150.9 MHz (14.1 T) 13 C spectra at thermal equilibrium. We use the subtraction procedure because relaxation effects were not taken into account in the calculation; consequently, we are unable to consider polarization decay due to relaxation at and during the field variation. To enable comparison of the experiment and calculation results, the amplitude of oscillations in polarization transfer traces were scaled with the same factor, then the starting polarization values were adjusted individually to give best agreement with experimental data.

Table 1 .

 1 Parameters used for energy calculations. Proton-carbon direct scalar coupling values marked by asterisk has been used in numerical simulations polarization transfer:

	Chemical shifts
	C γ	24.14 ppm
	C δ1	22.05 ppm
	C δ2	20.92 ppm
	Scalar couplings
	J(C γ -C δ1 )	35 Hz
	J(C γ -C δ2 )	35.4 Hz
	J(C δ1 -C δ2 )	0 Hz
	J(C γ -H γ )*	127.4 Hz
	J(C δ1 -H δ1 )*	124.8 Hz
	J(C δ2 -H δ2 )*	124.8 Hz

Acknowledgements

This work has been supported by the Russian Foundation for Basic Research (grants Nos. 19-29-10028 and 19-33-90251). I.V.Z. acknowledges support from the French embassy in the Russian Federation, in the framework of the Ostrogradsky fellowship (project No. 933824A).

Code availability

Scripts used to run numerical modeling of polarization transfer in carbon-13 subsystem of leucine in the presence and in the absence of proton decoupling are provided in supplement.