Stasis in the adaptive landscape of two ecologically similar 2 and sympatric damselfly species

Abstract


Introduction
began to be recognized as descriptors of the dynamics of adaptive zone-a collection of adaptive 86 landscapes inhabited by populations and species that share ecological mode of life (Simpson, 1944; 87 Hansen et al., 2008) rather than the dynamics of species evolving on a fixed adaptive zone. Thus, 88 the adaptive landscape estimated from present-day populations and the adaptive zone estimated 89 from a dataset that contains multiple species can be conceptually related to study the link between 90 contemporary selection and macroevolutionary dynamics. However, this idea is currently 91 considered a metaphor with limited bearings on empirical research, because many previous studies 92 have shown that phenotypic selection often fluctuates over short time scales (Kinnison & Hendry,93 2001; Grant & Grant, 2002;Gosden & Svensson, 2008;Siepielski et al., 2017) and that stabilizing 94 selection is not necessarily a common form of selection in nature (Kingsolver & Diamond, 2011).

95
Such findings are at odds with the core premise of the adaptive zones-the assumption that adaptive 96 phenotypes are stable over hundreds of millions of years. Therefore, to construct an operational 97 research program of the adaptive landscape "as a conceptual bridge between micro-and 98 macroevolution" (Arnold et al., 2001), we need to elucidate how dynamic adaptive landscapes are 99 at the population level, and if this does or does not result in a stable adaptive zone at 100 macroevolutionary time scale.

102
The stability of ecological niche has long been a major hypothesis to explain the persistence of 103 adaptive zone. In the classic work on the evolution of high-crown cheek teeth (hypsodonty) in 104 horses, Simpson, 1944 argued that the evolution of hypsodonty in the early to middle Miocene The insect order Odonata (damselflies and dragonflies) is a popular system for studies of 119 phenotypic selection because a rich background knowledge of natural history (Corbet, 1999)  We have three objectives in this study. First, we estimated univariate and multivariate selection on 139 two important fitness-related traits -body size and wing shape, and compared these within and 140 between species. Second, we quantified demographic and mating system parameters across 141 multiple populations of these two species. These demographic and mating system parameters 142 include the opportunity for sexual selection (Arnold & Wade, 1984), mating frequencies, male and 143 female densities, operational sex ratio (OSR) and the frequency of male-mimicking female color 144 morphs ("androchrome females"). These mating system parameters provide information about the 145 local sexual selection regimes, male-male competition (opportunity for sexual selection, OSR, 146 male density) and the intensity of male mating harassment and sexual conflict (mating frequencies 147 and the frequency of androchrome female morphs see Gosden & Svensson, 148 2009; Takahashi et al., 2014;Blow et al., 2021). Finally, by combining information from these 149 different sources about the ecological and social causes of selection, we aim to understand and 150 explain the stability and dynamics of the adaptive zones in these pond damselflies. waterbody, chase females that approach the waterbody, and compete for the opportunity for 166 mating. A male that successfully clasps a female will form the "tandem position", followed by the 167 formation of "a mating wheel" (Fig. 1D) during which egg fertilization takes place (Corbet, 1999 were captured using hand-held nets while slowly walking around waterbodies. Upon capture of 175 individuals flying without a partner, we examined sex and kept males and females separately in 176 small net cages. Individuals found as either a tandem or as a copulating couple were kept in plastic 177 cups. We visited these populations between the hour of 08.00 and 13.00 in all partially or fully 178 sunny days with temperature >15°C in May, June and July of 2020 and 2021. At each visit, we 179 sampled between 20 and 30 minutes, and between 3-5 people participated in the sampling. The 180 captured damselflies were kept in cooling bags to protect them from overheating, and were brought  cyathigerum. This dataset was used to estimate selection gradients presented further in this study.

188
In a complementary field study, we quantified social organization and mating systems of I. elegans 189 and E. cyathigerum. This was part of a community survey where we captured individuals of 190 damselflies and dragonflies regardless of the species identity along predefined transects, with the 191 aim of quantifying species composition, operational sex-ratio and density of local Odonata fauna.

192
At each visit, 3-5 people sampled between 10-20 minutes, who identified and recorded all captured 193 individuals in terms of species, sex, female color morph, sexual maturity and mating status (either 194 captured as single or copula/tandem). After recording this information, animals were released. We  severe in short-lived insects (Nishida, 1989) such as pond damselflies which typically live for only 217 a few days as reproducing adults (Corbet, 1999).   wing. Then, we tested the ability of ML-Morph to landmark wings accurately using shape.tester 245 function, where we compared automated landmarking by ML-morph to manual landmarks, which 246 responded to 99.2% precision in the landmarking using automation in a test set build from the 247 training set. Finally, we applied shape.predictor function to landmark the remaining image set that 248 have not been landmarked. All landmarked images were later checked manually for any errors, 249 and inconsistent landmarks across all images were corrected. We removed 357 images (178 for I. 250 elegans and 179 for E. cyathigerum) from the dataset before analysis because they contained an 251 injured/broken wing where position of 1 or more landmarks could not be deduced accurately.

252
Examination of erroneous images were done by one observer (AG).  coefficients were estimated for both species and sexes for all the traits using linear mixed effect 305 models with locale as the random factor. The partial regression coefficients of models that includes 306 only linear (i.e. unsquared) term were used to estimate b, while g was estimated as the partial 307 regression coefficients of models that includes both unsquared and squared terms. The quadratic 308 regression coefficients and associated standard errors were multiplied by two (Stinchcombe et al.,309 2008). Note that these operations are not necessary for quadratic terms between different traits (i.e. Analyses, selection gradients evaluated in these traits are on the mean-standardized scale without additional transformations. For comparison, we also estimated the variance-standardized selection 321 gradients by dividing them with standard deviation of the trait within species. We also performed 322 bivariate model to evaluate correlational selection between LD1 and size-PC1. We analyzed 323 forewings and hindwings separately in these analyses because LD1 of fore-and hindwings are 324 highly correlated (results not shown).

325
Estimation of selection gradients and statistical comparisons of adaptive landscape between 326 species were made using random mixed effect models implemented in lme4 package (Bates et al.,327 2015). We constructed three sets of models, which all included sampling location as the random 328 effect. First, for each sex and species separately, we modeled fitness component as the response 329 variable and one of examined traits (size-PC1, forewing shape LD1, hindwing shape LD1) and its 330 squared terms as the fixed explanatory variables. Second, we constructed two bivariate models, 331 for each sex and species separately, with the same model specifications as above but include either 332 size-PC1 and forewing shape LD1 or size-PC1 and hindwing shape LD1 and their interaction and 333 squared terms as the fixed explanatory variables. Finally, we constructed two models for each sex 334 separately, with the fitness component as the response variable, size-PC1 and forewing shape LD1 335 or size-PC1 and hindwing shape LD1 with their squared terms, species (either I. elegans or E.

365
Dimension reduction of body size and exploration of wing morphospace: A scatterplot of the 366 first two principal components of five body size measurements (Fig. S3) revealed that the first PC 367 axis explained the vast majority (93.0%) of the total variation in log of body size measurements.

368
Based on this, we hereafter use this axis (size-PC1) as our measure of body size. A scatterplot of 369 the first two principal components of the entire wing shape data (including both species, sex, and 370 both fore-and hindwings) is shown in Figure S4.  elegans. Subsequent PC axes up to PC5 explained 8.9%, 6.3%, and 5.9% of total variation, 378 respectively.

379
A linear discriminant function analysis (LDA) with species as a grouping factor revealed that the 380 axis that most effectively separates wings of the two species (wing shape LD1) explains 33.3% of 381 the total variation in wing shape and represents the variation in the width of the wing, that occurs 382 together with the stretch of two landmarks (LM9 and LM16) at the center of wing (Fig. 2). LD1 is 383 correlated with shape-PC1 (r 2 ; forewing: 62.5%, hindwing: 68.6%) and shape-PC2 (r 2 ; forewing: 384 86.9%, hindwing: 89.0%), but not with shape PC3, PC4, and PC5 (r 2 < 5%, Fig. S5). Within two 385 subsets of data that include either forewings or hindwings, LD1 explained 47.1% (forewing) and 386 57.4% (hindwing) of the phenotypic variation within each of the subset, meaning that LD1 is the 387 primary axis of phenotypic variation in these traits. In addition, LD1 captures variation that are 388 most relevant for the goal of our study to compare adaptive landscape between I. elegans and E. 389 cyathigerum. Moreover, given that our working hypothesis is that the two species reside at 390 different parts of a shared adaptive landscape, studying LD1 will make our results conservative 391 because we examine similarity of the adaptive landscape in the most divergent axis of wing 392 morphospace. For these reasons, the rest of our analyses will use LD1 as a representative wing 393 shape trait.

395
Selection gradients: Selection gradients evaluated at the species level are presented in Table 1. 396 Across both species, we found that selection favors females with large body sizes. Estimates of  Selection on wing shape LD1 was relatively weak across both sexes in both species (Fig. 3, Table   403 1). Estimates of b showed that the selection typically favors narrow and elongated wings in I. 404 elegans while broad and round wings are favored in E. cyathigerum. Considering the mean shape 405 differences between the two species (Fig. 2), these estimates suggest that, if all else being equal, 406 directional selection would drive the convergent evolution of wing shape. The quadratic selection 407 gradients (g) suggested stabilizing selection (i.e., negative γ) on wing shape LD1 in all trait-sex-408 species combination except for forewing of E. cyathigerum, although none of these negative estimates were significantly different from zero (Table 1). Therefore, wing shape appear to be 410 under weak stabilizing selection in both I. elegans and E. cyathigerum, and the result of b that 411 favors intermediate wing shape of the two species is consistent with the hypothesis that the two 412 species are evolving on different sides of a shared adaptive peak. For both size PC1 and wing 413 shape LD1, variance-standardized gradients showed qualitatively equivalent results (Table 1). 414 Results of correlational selections between wing shape and body size are presented in 415 Supplementary Material (Table S1, Fig. S6, S7).   By comparing linear (b) and quadratic (g) selection gradients at these population between I. 431 elegans and E. cyathigerum, we found no statistical evidence that the local fitness surfaces are 432 different between the two species (Fig. 4, Table S6). b on body size was far more variable than b 433 on wing shape (variance of b: male body size = 2.36, female body size = 10.1, male forewing 434 shape LD1 = 0.062, female forewing shape LD1 = 0.19, male hindwing shape LD1 = 0.021, female 435 hindwing shape LD1 = 0.027, unit: (relative fitness/loge(mm)) 2 for body size, (relative 436 fitness/centroid size) 2 for wing shape), which globally fits with selection estimates at the species 437 level that the adaptive landscape of wing shape was similar between the two species whereas that 438 of body size showed more variation. Because g was poorly estimated in most populations due to 439 small sample sizes (Table S2-S6), rest of analyses will be performed on b.

440
The summary of models relating b with social organization indices estimated at corresponding 441 location and year are shown in Figure 5, Table S7 and Table S8. We found that the strength of 442 directional selection on male body size is correlated with operational sex ratio (slope ± SE = 0.61 443 ± 0.15, p < 0.01) and opportunity for sexual selection (slope ± SE = 12.43 ± 4.95, p = 0.025), and 444 that a non-significant trend exists with respect to androchrome frequency (slope ± SE = -4,43 ± 445 2.57, p = 0.11). These relationships (Fig. 5) (Table S8 and Figure  S8). Although these relationships were somewhat weaker than those on b on body size, the way 454 by which social organization indices was related to the strength of selection was remarkably 455 consistent between I. elegans and E. cyathigerum, suggesting that the adaptive landscape of body 456 size in these two species are governed by similar ecological processes.   One question that we aimed to address here is whether the adaptive landscape remains stable over 478 long evolutionary time periods as has often been assumed in macroevolutionary models of 479 phenotypic evolution (Simpson, 1944;Hansen, 1997;Estes & Arnold, 2007;Uyeda et al., 2011). 480 Estimates of linear and quadratic selection on wing shape in I. elegans and E. cyathigerum revealed 481 that wing shape generally seems to experience weak stabilizing selection (Fig. 3, Table 1). There 482 was no evidence of divergent selection on wing shape between these two species, suggesting that 483 the general form in their adaptive landscapes of wing shape is relatively stable. This is quite 484 remarkable, considering that these two species belong to genera that diverged about 12 million shape in I. elegans were all negative in sign (Table 1). Based on the difference in wing shape 492 between E. cyathigerum and I. elegans (Fig. 2), these estimates can be interpreted as convergent 493 selection favoring an intermediate wing shape of these two species.

494
By plotting the population-specific selection gradients of wing shape together with the global 495 selective surfaces of E. cyathigerum and I. elegans (Fig. 6), we found that selection often favors a wing shape intermediate in between these two species, with a possible exception in hindwing of 497 females where b seems to push populations to species' specific adaptive peaks (Fig. 6D). However, 498 the apparent peak of E. cyathigerum reflects one population with small sample size (Höje A6, 499   Table S4) and the strength of selection is the weakest in hindwing among the three traits examined. 500 Therefore, our results are largely consistent with the idea that wing shape of E. cyathigerum and 501 I. elegans are evolving under a shared adaptive zone. It is conceivable that the intermediate wing 502 phenotype represents a global adaptive optimum determined by some fundamental mechanical and 503 physiological demands of flight, that would constitute the primary optimum (Hansen, 1997). Each 504 one of these two species may occupy a fraction of this adaptive zone according to their ecological 505 mode of life (Simpson, 1944) that are similar but also exhibit subtle differences (Fig. S9,   506 Supplementary Results). Our analyses illustrate that the adaptive zone can take stabilizing form of 507 selection even when individual selection estimates at any given place and time may give an 508 impression of directional selection (Table S2-  Compared to selection on wing shape, selection on body size was considerably more variable.

515
Females in both species experienced positive directional selection towards large body size 516 although this fecundity selection was stronger in magnitude in E. cyathigerum than in I. elegans 517 (Fig. 3, Table 1, Table 2). Thus, the fitness benefits in terms of fecundity increased faster with 518 female body size in E. cyathigerum than in I. elegans. Stronger sexual conflict and more intense 519 male mating harassment in females of I. elegans compared to E. cyathigerum can potentially 520 explain divergence in selection on female body size. Male mating harassment in I. elegans has 521 been shown to be mainly directed towards high-fecundity females (Gosden & Svensson, 2009). In 522 other insects such as Drosophila, large females have intrinsically higher fecundity but they also 523 suffer more from male mating attempts that reduce these intrinsic fitness advantages of large 524 females (Long et al., 2009;Chenoweth et al., 2015). In contrast to the situation in I. elegans, in E. 525 cyathigerum, female fecundity benefits of large size might be less affected by male mating 526 harassment because male mating harassment appear to be lower in this E. cyathigerum compared 527 to I. elegans based on several mating system indices (Supplementary Results, Table S9, Fig. S9). males and females typically co-occupy the same microhabitat and sex ratios are more even in this 537 species (Fig. S9). Sexual habitat segregation might therefore further release females of E. 538 cyathigerum from the cost of male mating harassment, enabling them to enjoy greater fecundity 539 advantages of large body size compared to I. elegans. 540 In males, we found divergent patterns of selection on body size between the two species (Fig. 3A). 541 Specifically, I. elegans is characterized by a weak (non-significant) stabilizing selection while E. 542 cyathigerum experienced disruptive selection (Fig. 3A). These results indicate that sexual selection 543 on male body size is highly context-dependent in these damselflies (see considerable variation in selection. However, we found no consistent difference in selection 556 evaluated at the local population level between I. elegans and E. cyathigerum (Fig. 4). This result 557 corroborates the similarity in these two species' adaptive landscapes of wing shape at the species 558 level. These data further suggest that the adaptive zone of wing shape has been stable and is largely 559 shared between these two species, and might have been so for millions of years. However, it 560 somewhat contradicts our finding that selection on body size evaluated at the species level seems 561 to be different between the two species. Below, we propose that ecology can reconcile this 562 discrepancy.

564
By investigating local sexual selection regimes on male body size, we found that selection on male 565 body size became stronger when the operational sex ratio (OSR) was male-biased (Fig. 5A), when 566 androchrome females were rare (Fig. 5B) and when the opportunity for selection (Arnold & Wade,567 1984) was high (Fig. 5C). In insects that engage scrambling male-male competition, OSR and the 568 opportunity for selection are important descriptors of mating systems, and have previously shown 569 to be correlated with sexual selection regimes in water striders (Arnqvist, 1992)   We also found that the local social environments affected the strength of sexual selection in I. 581 elegans and E. cyathigerum in a similar fashion, both in terms of body size (Fig. 5) and wing shape 582 (Fig. S8). These results suggest that the adaptive landscapes in males of these two species are   Figure 3: Univariate sexual selection (males, A-C) and fecundity selection (females, D-F) on body size (PC1) and wing shape (LD1) in E. cyathigerum and I. elegans. Fitness functions were visualized using cubic splines. Fitness data on the Y-axis were relativized at each species. The gray shaded regions around each spline represents the 95% confidence intervals. A) Selection on body size (PC1) in males, B) Selection on male forewing shape (LD1), C) Selection on male hindwing shape (LD1), D) Selection on body size (PC1) in females, E) Selection on female forewing shape (LD1) F) Selection on female forewing wing shape (LD1). The estimates for the coefficients of selection are presented in Table 1.  . Boxplots (barline in the boxplot represent 50 th percentile and the intervals represent the range between 1.5 times above 75 th percentile and below 25 th percentile) comparing linear selection gradients (b) between E. cyathigerum and I. elegans evaluated at 17 populations. b of A) male body size, B) male forewing shape, C) male hindwing shape, D) female body size, E) female forewing shape, F) female hindwing shape are compared. The estimates for the mean and standard errors of these parameters are as well as the statistical analysis for comparison between the two species (none are statistically different between the species) is presented in Table S6.  Table S2-S5). Note that one population (I. elegans from Hoje A7 in 2021, see Table S2 for estimates of this population) was removed from plots for better visualization, while the regressions are from analyses using the whole dataset. Opportunity for Selection C Enallagma cyathigerum Ischnura elegans Figure 6. Comparison of species-specific fitness surface with directional selection gradients at each locality (3 populations in E. cyathigerum and 14 populations in I. elegans). Sexual selection (males, A-B) and fecundity selection (females, C-D) on wing shape (LD1, both fore-and hindwings) are plotted. Fitness functions of each species (dashed lines) were visualized using cubic splines. Fitness data on the Y-axis were relativized at each location. Directional selection gradients (solid lines) were estimated using a regression model of relative fitness against LD1 at each location. Estimates of these models are presented in Table S2-S5. A) Selection on male forewing shape (LD1), B) Selection on male hindwing shape (LD1), C) Selection on female forewing shape (LD1) D) Selection on female forewing wing shape (LD1).