Non-parallel morphological divergence following colonization of a new host plant

Adaptation to new ecological niches is known to spur population diversification and may lead to speciation if gene flow is ceased. While adaptation to the same ecological niche is expected to be parallel, it is more difficult to predict whether selection against maladaptive hybridization in secondary sympatry results in parallel divergence also in traits that are not directly related to the ecological niches. Such parallelisms in response to selection for reproductive isolation can be identified through estimating parallelism in reproductive character displacement across different zones of secondary contact. Here, we use a host shift in the phytophagous peacock fly Tephritis conura, with both host races represented in two geographically separate areas East and West of the Baltic Sea to investigate convergence in morphological adaptations. We asked (i) if there are consistent morphological adaptations to a host plant shift and (ii) if the response to secondary sympatry with the alternate host race is parallel across contact zones. We found surprisingly low and variable, albeit significant, divergence between host races. Only one trait, the length of the female ovipositor, which serves an important function in the interaction with the hosts, was consistently different between host races. Instead, co-existence with the other host race significantly affected the degree of morphological divergence, but the divergence was largely driven by different traits in different contact zones. Thus, local stochastic fixation or reinforcement could generate trait divergence, and additional evidence is needed to conclude whether divergence is locally adaptive.


Introduction
Adaptation to novel niches through ecological selection has been identified as a major driver of speciation (Schluter 2008(Schluter , 2009Nosil 2012). Classic adaptive radiations including the Darwin Finches (Weiner 1994;Loo et al. 2019), the Hawaiian tarweeds and silverswords (Baldwin et al. 2003) and the recurrent evolution of limnetic-benthic stickleback pairs (Schluter 1993;Bay et al. 2017) are driven by adaptations to novel diets and habitats. Speciation driven by expansion of novel dietary niches has the potential to be rapid, like in Rhagoletis flies, where populations of hawthorn-feeding flies have adapted to use the introduced apple (Malus domesticus) as a host during the last 150 years (Bush 1969;Feder et al. 1988;Linn et al. 2004;Meyers et al. 2020). It remains an empirical challenge to resolve to what extent morphological adaptations are parallel and predictable under parallel selection regimes (Bolnick et al. 2018). Few studies have explicitly compared parallelism in multiple traits in such scenarios, but a recent study of Bahamas mosquitofish (Gambusia hubbsi) suggests that only a few traits are highly predictable across parallel predation environments (Langerhans 2018).
In addition to the need to better understand parallelism within certain environments, we further need to understand parallelism in response to co-existence with diverged lineages, host races or incipient species, as evolutionary divergence may be affected by the local presence or absence of such groups (Amarasekare 2003;Calabrese and Pfennig 2020). Reducing gene flow from populations adapted to other niches is crucial for ecological specialization to result in diversification, as ecological divergence does not necessarily lead to speciation (Taylor et al. 2006;Bolnick 2011;Lackey and Boughman 2017). Divergence and specialization may be reversed in the absence of reproductive isolation (Taylor et al. 2006;Lackey and Boughman 2017). Hence, if reproductive isolation between the incipient species is incomplete, secondary contact may break up important ecological adaptations. This is particularly evident when different populations of the diverging species are inhabiting discrete niches, like in the case of insect host races. Close congeners inhabiting different niches may be poorly adapted to the alternative environment (Rundle and Nosil 2005;Nilsson et al. 2017;Cronemberger et al. 2020;Martin et al. 2020), which may lead to reinforcement (Servedio and Noor 2003) and reproductive character displacement, with greater differences manifesting in traits involved in mate choice (Hinojosa et al. 2020). The expectations for reinforcement driven character displacement are, however, complicated by the frequency of non-parallel evolutionary responses at the morphological (Langerhans and Riesch 2013;Stuart and Losos 2013) and genetic levels (Stuart et al. 2017), meaning it remains a challenge to disentangle to which extent reinforcement selection results in parallel patterns of character displacement.
To clarify expectations for trait divergence in the case of parallel selection, both from the environment and from co-existence with a close congener, we outline three extreme scenarios of how populations may diverge in Fig. 1. These are expected under (a) parallel ecological selection (Johannesson et al. 2010), (b) parallel reproductive character displacement in response to selection against maladaptive gene flow with the other host race (Pfennig and Pfennig 2009), or (c) introgression (Wallbank et al. 2016). While a parallel response can be the result of parallel ecological selection acting on the study traits (Fig. 1a), parallel reproductive character displacement (Fig. 1b) or reduced divergence due to introgression (Fig. 1c) it could also arise for traits that are not strongly coupled with host plant use. While these are only three out of many possible scenarios (see also Supplementary Material 1), they are chosen based on expectations of strong parallel selection (Diegisser et al. 2008), expectations of reinforcement selection to reduce gene flow in sympatry due to the discrete genetic clusters host races maintain in sympatry (Pfennig and Pfennig 2009), and because of indications of some gene flow in genomic data (Ortega et al. [in prep.]) in the study system (see below).
Reinforcement selection is selection to reduce maladaptive mating in the presence of incipient populations. To understand whether reinforcement selection results in parallel morphological divergence to the same extent as ecological adaptations to a novel niche, it is necessary to compare patterns of divergence in allopatry and sympatry across geographic regions. Such studies would test the hypothesis that parallel evolutionary responses to reinforcement is less common than parallel ecological adaptations. The rationale for this is that we could expect traits under sexual selection to diverge in a more arbitrary fashion than adaptations of ecological traits, which are expected to diverge in parallel across populations of similar niches (Fig. 1a). In particular, studies asking how co-existence with congeners affects trait distributions, e.g. through parallel reproductive character displacement (Fig. 1b), may provide insight into the nature of the selection pressures exerted by co-existence of close congeners. Here, we address how morphological traits evolve and diversify in allopatric and sympatric areas using two host races of the fly Tephritis conura. This species has recently undergone a host shift resulting in host races specializing on different Cirsium thistle species as larval food plants (Diegisser et al. 2006a(Diegisser et al. , b, 2007(Diegisser et al. , 2008. Macroevolutionary analyses reveal that host plant driven speciation, as in the case of T. conura, is one of the main factors explaining the tremendous diversity of phytophagous insects, as many speciation events can be attributed to host plant shifts (Berlocher Fig. 1 Hypothetical predictions for morphological divergence under different scenarios. Three hypothetical scenarios expected under extreme parallel ecological selection, extremely parallel character displacement in response to reinforcement selection, and introgression to an extent that populations converge. These are only three of many possible patterns of divergence (see Supplementary Material 1). Bars represent mean values of morphological traits, with one population being represented by each bar. While the white areas of the graph depict allopatric populations, the grey area includes co-existing populations. 'N. allo.' stands for northern allopatry, 'Symp.' stands for sympatry and 'S. allo.' stands for southern allopatry. Two hypothetical transects are depicted side by side. A: Scenario 1, parallel ecological selection between host races. This scenario would be expected e.g. for ecologically important traits where natural selection imposed by the environment provided by the host differ. B: Scenario 2, parallel reproductive character displacement. This pattern could arise e.g. if a character under sexual selection would be reinforced by maladaptive hybridization in a parallel manner in both transects. C: Scenario 3, hybridization in sympatric areas is causing introgression resulting in reduced divergence in these areas. Whereas these explanations are not the only possible explanations, they were chosen as we expect some degree of parallelism to result from the adaptation to a similar environment and potentially some gene flow in the sympatric regions. Alternative explanations to the pattern in C would be clinal host plant independent variation, for instance 1 3 and Feder 2002;Dres and Mallet 2002;Nylin et al. 2014). Moreover, diversification rates are elevated in herbivorous insects compared to their non-herbivorous relatives (Mitter et al. 1988;Farrell 1998;Wiens et al. 2015). For an insect specialist like T. conura, the host plant provides a discrete environment imposing a multidimensional selection pressure. Hence, insect host plant adaptations are predicted to affect multiple traits involved in phenological matching, female host preference and larval performance (Matsubayashi et al. 2010). The putatively strong selection pressures involved in host plant adaptation provide a strong prediction of parallelism in host-use related traits of phytophagous insects at early stages of the speciation process (Nosil et al. 2002;Meyers et al. 2020). However, the prediction is less straightforward for responses to sympatry, because the traits involved in putative reinforcement are not necessarily associated with host use.
Tephritis conura provides an excellent system for studying parallelism in divergence in response to ecological adaptation and co-existence with congeners, as it enables comparisons of trait variation across replicated sympatric and allopatric settings. The two recently established host races inhabit two geographically isolated sympatric zones with adjacent allopatric populations, one East and one West of the Baltic Sea (Fig. 2a). Tephritis conura females of the two host races are ovipositing into buds of either C. heterophyllum (the ancestral host; Fig. 2c) or C. oleraceum (the novel host; Fig. 2d), and larvae show host specific performance and survival (Diegisser et al. 2008). Hence, the system enables investigating response to sympatry across two geographically separated contact zones. Importantly, the strong ecological selection brought on by the host plants (Diegisser et al. 2007(Diegisser et al. , 2008 led us to expect a degree of parallelism between separated populations. Moreover, the two host races are consistently genetically differentiated from each other, both in allopatry and sympatry, with host races being more diverged from each other than populations of the same host race from different geographic regions (Diegisser et al. 2006a, b;Ortega et al. [in prep.]). The highly consistent genetic divergence (Diegisser et al. 2006a), the documented strongly divergent ecological selection pressures acting on larval ability to process the host plant tissue (Diegisser et al. 2008), the divergent phenologies (Romstock-Volkl 1997), and the differences in host plant preference (Diegisser et al. 2008) all suggest that these host races could potentially be forming new species. Based on these ecological and genetic differences we predict that the host races should be morphologically divergent, and this divergence to be parallel. Furthermore, as larval survival is strongly reduced in the alternate host (Diegisser et al. 2008), we predict hybridization to be maladaptive, potentially resulting in reinforcement and reproductive character displacement between the host specialists in sympatry. Potentially, mate choice is based on traits under divergent ecological selection (Thibert-Plante and Gavrilets 2013) leading to parallel character displacement. In contrast, traits not involved in ecological adaptation to the host plants may follow a scenario similar to a mutation order scenario, where traits not necessarily involved in reproductive isolation arbitrarily diverge depending on the order of novel mutations and stochastic fixation (Mendelson et al. 2014). Here, we specifically test the hypotheses (i) that there are consistent morphological responses to a host plant shift across the two transects separated by the Baltic Sea and (ii) that we find parallel responses to sympatry in the two contact zones if ecologically selected traits are used to discriminate against the alternative host race to avoid maladaptive gene flow. As we lack replication within the Eastern and Western transects we will, however, not be able to differentiate between transect specific character displacement following a mutation order scenario and genetic drift.

Study species and sampling
The dipteran T. conura infest several species of the thistle genus Cirsium (Asteraceae) (Romstock-Volkl 1997). Adult T. conura oviposit into thistle buds during early summer, wherein the adolescent flies remain during larval development and pupation. Tephritis Principal component analyses are used for illustrating how the host races of T. conura differ in morphology using an unguided approach. Both morphological traits including including body length, ovipositor length, wing length, wing width and melanisation ratio as well as wing shape traits in form of six relative warps were included. E: For female flies, where the ovipositor also was included in the analysis. F: For male flies conura infesting C. heterophyllum, the melancholy thistle, have recently colonized and adapted to C. oleraceum, the cabbage thistle. Throughout the manuscript, we refer to flies infesting C. heterophyllum as CH-flies, while the flies infesting C. oleraceum are referred to as CO-flies for simplicity. Haplotype analyses suggest that a peripatric host shift took place in the Alps during the last ice age (Diegisser et al. 2006b). The two host races are largely reproductively isolated, but there is evidence of some small amounts of gene flow (Diegisser et al. 2006b). The CO-flies have adapted ecologically to the smaller bud sizes, as reflected in their significantly shorter ovipositor to body size ratio (Diegisser et al. 2007). Another variable character is wing pigmentation. Wing patterns in tephritid flies have been suggested to be under sexual selection, as males perform dances with their wings to attract females (Sivinski and Pereira 2005). Tephritid males attempting to initiate copulation situate themselves in front of female flies and posture with their wings (Sivinski 2000; but see Briceno and Eberhard 2017). While there is no direct evidence to support that wing traits affect sexual selection in T. conura we included it within a range of traits with the expectation that they could to be under different selective regimes, to investigate to which extent different traits diverge in parallel ( Fig. 1a,b). As we know that the ovipositor is a key trait involved in host plant adaptation (Diegisser et al. 2007), we focus on the findings from female T. conura in the manuscript, while the results for males are more briefly referred and readily available in the Supplementary Material.
We used a parallel sampling design to examine phenotypic adaptation to a novel host plant and effects of co-existence between host races, by sampling each host race both in sympatry and allopatry on each side of the Baltic (Fig. 2a). We collected thistle buds infested by T. conura larvae/pupae and allowed the adults to eclose in a common environment. CO-fly larvae were sampled in Germany and Lithuania (allopatric areas), and both host races were collected in sympatric areas in southern Sweden and Estonia. Allopatric CH-fly larvae were sampled in central-Sweden and Finland (Fig. 2a, Table S1). All sampling took place during June and July 2018. We define sympatric as presence of both thistle species in a region and allopatric as only one thistle species being present in this study (Fig. 2a). The sampling scheme is designed to enable examining to what extent patterns of phenotypic divergence are explained by host plant adaptations and by parallel adaptations to co-existence with the other host race. It does, however, not, allow us to differentiate transect specific character displacement from genetic drift, since differentiating these processes would require more than one sympatric population in each transect to be sampled.

Morphological measurements
Tephritis conura adults eclosed from field-collected thistle buds in a common lab environment (see Supplementary material S2). One male and one female per bud were euthanized by freezing a few days after emergence, and subsequently included in the morphological analysis. For each individual, we took magnified photographs using a Celestron 44,308 USB microscope. We photographed a lateral image of the fly body after removal of the wings and a dorsal image of the right wing on a transparent background to allow better visibility of the wing veins. Body length and ovipositor length were measured digitally from lateral photographs (Fig. S1). We placed 14 landmarks, adapted from Pieterse et al. (2017), digitally on the dorsal wings (Fig. 2d, S2) for geometric morphometrics (Zelditch 2004). We added a landmark; number 15, to reflect the high variance in the proximal area on the wing. Digitization was performed in TPSDig2 v2.31(Rohlf 2015) and we used TPSUtil v1.76 (Rohlf 2015) for file handling. Wing melanisation was measured with an automated script without any user queries developed in MATLAB (Matlab 2017). As the wings do not display any chromaticity, analysis is based on the red spectral band only. The script extracted the intensity of the red spectral band for each pixel, and performed white calibration by division with the statistical mode corresponding to the white background of the images. Subsequently, images were inverted to represent absorbance rather than reflectance, so the melanised wing would extend from the white background (set to zero). The script identified the wing and separated melanised areas from non-melanised areas, and the size of these areas were divided to estimate the fraction of the wing that was melanised.

Statistical analysis
We used PAST3 v3.20 (Hammer et al. 2001) to apply a Procrustes fit to the landmark data to align and scale the wings (Fig. 2d). To produce relative warps (i.e. principal components of shape) to compare shape between groups, a wing-shape principal component analysis (PCA) was performed with the Procrustes fitted data using PAST3 v3.20 (Hammer et al. 2001) (Fig. S3). Based on the variance explained by the eigenvalues (Fig. S4) and the broken stick criterion (Jackson 1993), six relative warps, principal components of shape, which jointly explain 68% of the wing shape variance were identified (Fig. S5). These relative warps were included in subsequent analyses of phenotypic divergence to represent wing shapes. All subsequent statistical analysis were performed in the statistical software R (R Core Team 2019).
We quantified five morphological traits (body length, ovipositor length, wing length, wing width and melanisation ratio) in addition to wing shape (represented by relative warps produced from landmark analysis) for 583 flies from eight populations (see Supplemental  Table 1 for population sizes). As exploratory analyses, aimed to provide illustrations of how the morphology of host races differed in multivariate space, we performed one full Principal Component Analysis (PCA) per sex. We used the variables fly body length (mm), wing length (mm), wing width (mm), wing melanisation ratio (%) and relative warps 1-6 reflecting wing shape in the analyses for both sexes, but also included ovipositor length (mm) in the analysis of females.
To formally test if the two host races were significantly differentiated, we applied two multivariate analyses of variance (MANOVA), with all variables measured included as response variables (body length, ovipositor length, wing width, wind length, wing melanisation and PC1-6 of wing shapes). The first one addressed if the differences between host races were affected by co-existence with the other host race as we included both host race, co-existence and their interaction as factors. Second, to test if the patterns of morphological adaptation were parallel in the Eastern and the Western transects, and explicitly address if host plant affected morphology in the same way in these replicates, we used a MANOVA with host race, geographical setting and their interaction as factors. These MANOVAs were both performed separately on females and males as this enabled including the biologically important trait ovipositor length (Diegisser et al. 2007) in analyses of females. This division of sexes also reduced the multicollinearity of explanatory factors to below recommended values (Hair 2010). To identify if specific traits are divergent between host races, geographic contexts and transects, we followed up on these MANOVAs with post hoc analyses consisting in ANOVAs per each individual trait. As we could expect to see strong correlations between size traits, we also performed analyses of the residuals of ovipositor length corrected for body length. Finally, we also performed vector analyses capturing the difference in how host races have evolved in trait space to address if the host races have adapted to co-existence and geography in a parallel manner. By comparing P-matrix vector length differences and angles, it is possible to infer the magnitude (corresponding to the length of the vector) and direction of multivariate trait divergence (Oke et al. 2017). The length of a divergence vector reflects the magnitude of divergence across all traits, while the angle between two divergence vectors reflect similarity in patterns of divergence across traits (a small angle indicates that traits have diverged in a similar way). We performed analyses of whether the host races were differentially affected by co-existence, reflecting how much they have diverged and whether the same traits had diverged in parallel in the sympatric population compared to the allopatric population for both host races. We also examined if the magnitude or angle of divergence between host races differed between the Eastern and Western transect by using host race and transect as factors in a vector analysis. These two types comparisons were performed separately for males and females.
To further investigate parallelism in differentiation of the host races in the two transects, we applied linear discriminant function analyses (LDAs) separately on the data from the Eastern and Western transects using population as a grouping factor. These analyses were performed seperately on males and females, and reveal which groups are most differentiated and what traits distinguish these groups. To test if the patterns of divergence differed significantly between transects, we performed each LDA 10,000 times for each sex using the bootstrap R package 'boot' (Canty and Ripley 2020), and used the confidence intervals to assess if the loadings differed between analyses. To assess discreteness of the clusters, we also estimated the percentage correctly classified individuals per group using confusion matrices. To further test whether the host race differences found in one transect also result in correct classification of flies in the other transect, we used the trait loadings of the traits that discriminate the Eastern populations to classify the Western individuals and vice versa. Finally, we assessed the proportions of divergence that is shared among host races and unique among the populations, respectively. To this end, we used one nested MANOVA per sex, including size and shape variables with host race and transect as independent factors following Langerhans and DeWitt (2004).
Interestingly, host race differences in females depend on co-existence, as illustrated by a significant interaction between these factors, driven by that CO-fly females become more morphologically similar in sympatry compared to CH-fly females (Pillai's trace = 0.155, F 11, 271 = 4.52, p < 0.001; Table 1; Fig. 3 and S6). Similarly, there was a significant interaction between co-existence and host race for males, with CO-fly males becoming more similar to CH-fly males in sympatry (Pillai's trace = 0.141, F 10, 275 = 4.5, p < 0.001; Table 1; Fig. S7-S10). The differences between host races were also affected by geography, as illustrated by the significant interaction between host race and transect (Pillai's trace = 0.19, F 11, 271 = 5.78, p < 0.001; Table 1; Fig. 3 and S6), illustrating that females of the two host races were more diverged in the Western than the Eastern transect. This pattern holds true also for males (Pillai's trace = 0.16, F 10, 275 = 5.37, p < 0.001; Table 1; Fig. S7-S10). Only ovipositor length and wing shape (Table S2) differed significantly between females of the two host races in the post hoc tests for the analysis including coexistence, while body length, ovipositor length, wing length and wing shape (Table S2) for the analysis including geography. We also found that coexistence affected host races differently for body length, wing width and wing shape (Table S2) and that geography affected host races differently for body length, ovipositor length and wing shape (Table S2). All other results and findings for males are reported in Supplemental Table 2. The findings for ovipositor remained the same also when corrected for body length (Fig. S11).
Finally, we find that sympatric and allopatric populations show approximately equal differences between females of the two host races, as very similar vector lengths are retrieved when applying a vector analysis (sympatric: 0.54% difference; allopatric: 0.65% difference). The traits that differ between sympatric and allopatric populations are not the same, as the angle between the host race vectors is 67.1°. This mostly reflects that body length Table 1 Models of host race divergence (A for females, C for males) including coexistence and its interaction with host race and (B for females, D for males) including transect and its interaction with host race MANOVA models to investigate the effect on multivariate morphological divergence in five morphological and six wing shape traits (A, C) using host race, co-existence and their interaction and (B, D) using host race, transect and their interaction. These analyses were performed separately for females (A, B) and males (C, D) to include the biologically important ovipositor for females differences were larger in sympatric populations. When examining if divergence between host races are affected by geographic settings for females, we find small differences in the transect specific overall magnitude of divergence between host races as vector lengths were similar across (East: 0.36% difference; West: 1.03% difference). Which traits that have diverged differs between transects though, as the angles of the vectors differ by 59.45°. This difference mainly reflects larger body length differences between the host races in the Eastern transect compared to the Western transect. The findings for males were overall similar but showed a larger angle of divergence between allopatric and sympatric populations (Table S3). As host race affected morphology differently in the two transects, we further tested if the major axis of divergence separated host races in both transects, and if the same traits separated groups using one Linear Discriminant Analysis (LDA) for each sex. This analytical approach revealed that the importance of host race for population separation differed between transects. In the Western LDA, host races separated along the first discriminant axis and the sympatric and allopatric populations separated along the second discriminant axis for both females (Fig. 4) and males (Fig S12). In the Eastern LDA, the first discriminant instead divided the two CO-fly populations, and the second discriminant axis divided the two CH-fly host race populations for females (Fig. 4). Males show a similar result for CO-flies, but sympatric CH-flies land intermediately to other populations and allopatric CH-flies separate on the second axis (Fig S12). Based on the confusion matrix, 76.13% and 72.5% of the CH females and males and 81.94% and 75.29% of the CO females and males, respectively, were correctly classified in the Eastern transect, whereas 80.94% and 74.76% of the CH females and males and 78.41% and 76.43% of the CO females and males, were correctly classified in the Western transect (Table S4). To further test whether the host race differences found in one transect also result in correct classification of flies in the other transect, we used the trait differences discriminating the Eastern populations to classify the Western individuals and vice versa. We find that the trait loadings from the Eastern LDA enable us to classify 28.78% and 31.78% of Western females and males to the correct host race, whereas the trait loadings from the Western LDA correctly classified 32.88% of the Eastern females and 35.22% of the Eastern males (Table S4). These analyses LDAs illustrating differences in how the host races group along the first two linear discriminant axes in Western (a) and Eastern (b) flies. c: The morphological traits loading on the two first discriminant axes for Eastern and Western flies. Colors illustrate how much the standard error diverges from zero based on 100 000 bootstrap replications. Loading that surpass zero are depicted in red colors whereas loadings significantly lower than zero are colored in blue. Plots and analyses are based on female fly morphology. Findings based on males are reported and illustrated in Supplementary Fig. 12 also revealed that the host races from transects East and West of the Baltic had diverged in different sets of characters, as bootstrap loadings from the two LDAs show that different characters loaded on the first discriminant axes (Table S5). For the Eastern transect, body length and wing warp 6 had positive loadings, and wing size, ovipositor length and wing warp 1 had negative loadings on LD1 in females (for males see Table S5). Wing warp 5 loaded positively on LD1 for both transects, but the loading was substantially stronger for Western populations. In the Western transect, LD1 had positive loadings for wing width and shape. Females also differed in which traits that loaded strongly on LD2. In the Eastern transect, wing shape, body length and wing width had negative loadings, but wing length positive loadings on LD2, while for the Western transect body length, ovipositor length and wing lengths loaded positively (see Table S5 for males).
To assess how much of the female divergence was shared among all populations of one host race compared to divergence unique to one transect, we also estimated Wilk's partial η 2 . While a high share of the partial variance was explained by shared divergence between host races (F 11.267 = 8.42, 25.75%) and transect specific patterns of host race divergence (F 11.267 = 5.86, 19.45%), divergence between transects explained the highest percentage of partial variance (F 11.267 = 16.73, 40.81%. For full results of females and males, see Supplemental table S6.

Discussion
Contrary to previous studies that found clear differences between the T. conura host races (Diegisser et al. 2007), we found low and variable morphological differentiation. The multivariate differentiation between females of the two host races was driven by subtle differences across many traits, and the traits that were divergent between host races differed depending on whether they co-existed with the other host race or not, as well as between the two transects. Previous work on the T. conura host races has reported consistent divergence for loci inferred to be involved in host plant adaptation (Peptidase D and Hexokinase; see (Diegisser et al. 2006b), and whole genome data (Ortega et al. [in prep.]) support the presence of discrete genetic clusters for each host race. Moreover, the flies show very poor performance on the alternative host plant (Diegisser et al. 2008). In light of this, the utterly moderate and variable host race divergence in morphology was unexpected. However, our findings are in line with the observation that parallelism in fitness is typically higher than parallelism in both phenotypic divergence and genetic divergence (Bolnick et al. 2018). Patterns of parallel divergence only in traits under very strong ecological selection is consistent with the findings of parallelism in response to different predation regimes in Bahamas mosquitofish (Langerhans 2018), where only a few traits show highly predictable patterns of diversification. Hence, at early stages of diversification driven by ecological adaptation, parallelism may be high only for traits that are strongly coupled to the ecological factor that adaptively diversify the incipient species.
We find that one out of five female traits, ovipositor length, match a parallel divergence scenario. Hence, the degree of parallelism among host races in morphological divergence was lower than expected. Ovipositor length was significantly diverged between host races, but the strength of divergence differed between sympatric and allopatric populations as well as between the Eastern and Western populations. In Western populations, the sympatric CH-fly population has intermediate ovipositor lengths, whereas CO-flies and Eastern CH-flies show parallel ecological divergence between host races both in sympatry and 1 3 allopatry, corresponding to the expectations under parallel ecological selection (Fig. 1a). The intermediate measurements for traits found in the Western CH-flies could suggest selection for shorter ovipositors in sympatry or reflect introgression (Wallbank et al. 2016;Llaurens et al. 2021). We find host race divergence in body size, with clear host race differences in body length in the Western transect, while Eastern flies have a shuffled distribution of body sizes. Size typically varies with temperature in insects (Atkinson 1994; but see Shelomi 2012). While the parallel sampling design would partly correct for temperature effects, temperature could still differ between the Eastern and Western transects. Alternative explanations to the lack of parallelism include consistent divergent selection only on a few traits (Lande 1979) or strong genetic covariances of some traits preventing divergence (Arnold 1992). Given the ecological importance of an ovipositor with the optimal length for the specific host plant, other traits could be constrained by the balancing selection on the ovipositor optima. Lastly, different levels of gene flow have previously been found to interact with (non-)parallel responses to similar environments (Stuart et al. 2017), and studying gene flow between host races in the T. conura system would be an interesting future line of inquiry that could shed light on whether there is evidence for asymmetric gene flow.
When we investigated effects of co-existence and geographic origin and their separate interactions with host race divergence, we found that both co-existence with the other host race and whether the population originates from East or West of the Baltic Sea affected host race differences (Fig. 3). Hence, the low host race divergence may partly be explained by interacting effects of co-existence with the other host race, and with non-parallel patterns of divergence in the two transects, which could indicate a mutation order scenario (Mendelson et al. 2014) or genetic drift. If reinforcement due to maladaptive hybridization would have been a prominent force acting on the T. conura host races, we would expect traits to be more divergent between host races in sympatry compared to allopatry, i.e. for reproductive character displacement to arise (Comeault et al. 2016;Calabrese and Pfennig 2020;Kyogoku and Wheatcroft 2020;Fig. 1b). This was generally not the case, and we found no evidence for parallel reproductive character displacement. Instead, the most common pattern was that traits in sympatric populations are more similar between host races than traits in allopatric populations. Contrary to our expectations, based on the fact that tephritid flies perform a mating ritual that includes elaborate wing movements from both sexes during courtship (Sivinski 2000; but see Briceno and Eberhard 2017), this pattern of convergence in sympatry was strongest for wing morphology traits. This is the opposite of what would be predicted under a reinforcement scenario but consistent with introgression ( Fig. 1c), as found by Ortega et al. [in prep.], or effects of a shared environment (Llaurens et al. 2021). An intriguing future direction would be to include further populations of each transect, to test patterns of reproductive character displacement with a more extensive sampling scheme.
Wing size was most strongly affected by co-existence. Interestingly, wing length was overall significantly shorter in sympatry than in allopatry, consistent on both sides of the Baltic Sea. Although these differences are small, we could speculate that the disproportionally shorter wing length in relation to body size could be consistent with selection for shorter dispersal distances (Claramunt et al. 2012), potentially to avoid dispersal to the alternative host plant which could lead to maladaptive introgression, though formal tests of wing length and dispersal would be required to test this. We found patterns that could be consistent with trait convergence in the Eastern transect where sympatric CH-flies have broader and sympatric CO-flies narrower wings, but this pattern was not parallel. Therefore, while this observation could be consistent with a mutation order scenario where the difference may be important but the specific trait or direction of the difference is arbitrary (Mendelson et al. 2014), our experimental design does not enable us to differentiate this from genetic drift. An alternative explanation to that phenotypes are more similar in sympatry is that sympatric flies of the two host races are exposed to a shared environment, selecting for similar morphological adaptations (Grenier and Greenberg 2005). Given the weak parallelism in morphological traits in allopatry, the degree to which host race differences could be explained by co-existence with the other host race is expected to be low. Still, small and variable differences in sexually selected traits could result in amplified and parallel patterns of divergence in these traits (Blankers et al. 2019). For example, if the size or direction of the G-matrix (the genetic covariance among traits) differ consistently between host races (Arnold et al. 2008), then the lines of least resistance along which populations can evolve (Schluter 1996), or the autonomy of how traits can evolve in relation to each other (Bolstad et al. 2014), may differ between host races. Prior knowledge of traits under sexual selection in this system would be valuable to guide hypotheses and studies on this subject. The strong selection against using the wrong host plant (Diegisser et al. 2008) and consistent genetic differentiation between host races Diegisser et al. 2006a, b;Ortega et al. [in prep.]) suggest that there could be reinforcement selection, further underlining the need to investigate whether there is assortative sexual selection at work in this system and which traits that are affected by it.
Contrary to studies of parallelism in sexually selected traits in ecotypes in e.g. sticklebacks (Boughman et al. 2005), we find mixed evidence for parallelism. In the Western transect, shared host race divergence constitutes the main share of variation as reflected by the host races being separated by the main discriminant axis. Conversely, host race explains less of the variance in the Eastern transect, where the first discriminant axis separates the two CO-fly populations. The traits that differ between host races also differs between our two geographic replicates. The differences in patterns between the transects East and West of the Baltic Sea could have several additional or alternative explanations related to demographic history (Tregenza et al. 2000), population size (MacPherson and Nuismer 2017) or the extent to which the host plant races co-exist locally (Rundle and Schluter 1998). Possibly, one contact zone may be older than the other and populations in older sympatry would have had more time for reproductive character displacement to develop (Johannesson et al. 2020). Alternatively, if the Eastern transect has a higher proportion of suitable thistle habitats, this could have increased both within-and between host race connectivity and potentially gene flow (Servedio and Noor 2003). Genetic data and detailed analyses of introgression should be used to resolve whether selection against hybridization could be expected.
Morphological differentiation does not always strongly reflect even crucial ecological adaptations. For instance, cultural evolution contributes to reproductive isolation in Cassia crossbills, Loxia sinesciuris (Porter and Benkman 2019) and Rhagoletis pomonella have adapted their phenologies to host fruits ripening at different times of the year (Filchak et al. 2000). Another potential explanation to the low consistency of host race divergence may be the traits included in this analysis. Females were more diverged than males likely as ovipositor length was differentiated, consistent with previous findings (Diegisser et al. 2007). These findings are similar to those of Jourdan et al. (2016) where divergence in female mosquitofishes (Gambusia) was more parallel than male divergence. Potentially, the other traits measured are not important enough for host plant adaptation to result in strongly parallel divergence. Importantly, both maternal effects (Mousseau and Fox 1998;Bernardo 2015) and environment driven phenotypic plasticity (Moczek 2010;Clissold and Simpson 2015), or host plant quality (Egan and Ott 2007) could affect our findings, especially for size traits (Steiger 2013). If these effects differ between populations, they could potentially reduce parallelism stemming from genetic differences.
Traits that have been shown to differ strongly between host races include female ovipositor length (Diegisser et al. 2007) and the larval ability to survive on the different host plant species (Diegisser et al. 2008). A performance experiment showed low viability of T. conura larvae adapted to C. heterophyllum when reared on C. oleraceum (Diegisser et al. 2008). Thus, the larval ability to process plant tissue is likely under strong selection and expected to show similar patterns across populations. Furthermore, host plant preference may have a potential to act as a magic trait (Gavrilets 2004;Thibert-Plante and Gavrilets 2013) in T. conura, separating the habitats of the populations and simultaneously providing reproductive isolation as these flies mate on their host plant (Diegisser et al. 2007). Finally, we cannot rule out that other, unmeasured traits are important for sexual selection, as pheromones have been suggested to play a role in tephritid mate choice (Roriz et al. 2019), as well as overexpression of antioxidants, which have been shown to increase male performance under certain conditions (Teets et al. 2019). Regardless of the exact selection pressures acting on our set of study traits, our findings are important, because they show a context dependence for host race adaptation. Parallelism was found only in ecologically strongly selected traits, whereas we found non-parallel changes for many other traits. These insights should guide the design and interpretation when studying ecologically driven divergence.
In conclusion, our work suggests that morphological responses to niche shifts can be highly context dependent. Co-existence with closely related congeners and demographic origin may affect easily measured morphological characters, potentially masking underlying parallelism in traits important for adaptation to specific niches.