Biogeography and ecology of geographically distant populations of sibling Cryptocephalus leaf beetles

Different populations of two closely related species, Cryptocephalus ﬂ avipes and C. bameuli , from western (Alps, Apennines and Pyrenees) and central Europe (Poland, Ukraine and Pannonia) were analysed. On the basis of DNA sequences from two genes, cox1 and ef1- α , distinctiveness of both species was con ﬁ rmed. Nevertheless, possible hybrids were identi ﬁ ed in Carpathian mountains. We found a signi ﬁ cant genetic differentiation among populations of C. ﬂ avipes and C. bameuli from distant regions but a high genetic similarity between populations of C. bameuli from north and south of the Carpathians. Demographic estimates suggest a past population expansion in the case of C. bameuli and a recent one for C. ﬂ avipes , possibly occurred during Pleistocene and Holocene, respectively. Distribution modelling showed that C. ﬂ avipes is typically present in the mountain systems, whereas C. bameuli is associated with hilly areas of central and eastern Europe. Based on the present data, Last Glacial Maximum refugia of both species were located in the Alpine region and Black Sea coasts, but on different elevations. The characterization of the insect diet, through a DNA metabarcoding approach targeting the trnL plant intron, demonstrated a signi ﬁ cant differentiation of food preferences between the two species, as well as between geographic populations within the species.


Introduction
Cryptocephalus Müller (Chrysomelidae) is one of the most species-rich genera in the order of Coleoptera, and because of that, it represents an excellent model for studying ecological and evolutionary processes, in particular speciation. Recent studies focused on leaf beetles have revealed interesting evolutionary features of Cryptocephalus, as the presence of highly divergent lineages (i.e. Cryptocephalus coryli (L.), C. decemmaculatus (L.), C. nitidulus (F.) and C. barii Burl.) (Piper & Compton 2003;Brunetti et al. 2019); cases of hybridization and mtDNA paraphyly in the only two species complexes studied so far (i.e. C. sericeus and C. flavipes species complexes; Gómez-Zurita et al. 2012; Montagna et al. 2016a). The latter finding was explained by the authors as the result of recent cladogenetic events and by episodes of hybridization in areas of secondary contact among the closely related species.
The C. flavipes complex represents an interesting case study for population genetics and ecology (e.g. the host plants preferences) of sibling species (Lopatin & Nesterova 2002). C. flavipes complex includes eight taxa with highly similar external morphologies; the species have medium size (2.5-5.0 mm) and are distinguishable by slight differences in the shape of the spermatheca and by patterns of yellow stripes on the pronotum and elytra margins (Duhaldeborde 1999;Warchałowski 2010). Recently, the shape of the metaepisternum has been demonstrated to be a highly effective character in differentiating among the species of the complex (Montagna et al. 2016a). At present, the following species are listed within the complex: C. signatifrons Suffr. (distributed in Europe), C. turcicus Suffr. (present in the Mediterranean Basin and the Near East), C. peyroni Mars. (Near East), C. bameuli Duha. (Europe to West Asia), C. flavipes F. (Europe to West Asia), and C. alborzensis Rapi. (Near East). Recently, integrating molecular and morphological data, C. quadripustulatus Gyll. (distributed in Europe) was assigned to the C. flavipes group and the existence of a new taxon from Turkey was demonstrated adopting an integrated taxonomic approach (Montagna et al. 2016a), which it was formally described as C. cilicus by Duhaldeborde (2018).
Cryptocephalus flavipes and C. bameuli inhabit similar environments such as the xeric grasslands and the scrublands. However, the localization of known populations of the two species suggests that C. flavipes is more associated with dry scrublands and woodland edges of mountain areas, whereas C. bameuli inhabits mainly steppes and xeric grasslands. Because of their similar phenotypes and ecology, in association with the hybridisation hypothesis (Montagna et al. 2016a), these species can be considered excellent subjects to investigate evolutionary processes that led to their separation, as the pastpresent population diversity and demography, host plants adaptation and microbiota associations, as postulated in the case of the C. marginellus species complex (Montagna et al. 2015b).
Using molecular data and species distribution modelling analyses, three main goals have been addressed in the present study: i) estimate the genetic diversity and distinctiveness among geographically distant European populations of C. flavipes -C. bameuli; ii) infer their demographic history during the Quaternary climatic fluctuations, in order to evaluate the responses of the two species and their potential glacial refugia; iii) identify the range of host plants exploited by the species, through a DNA metabarcoding approach.

Sampling
We used the collection of specimens of C. flavipes and C. bameuli from the majority of their ranges in Europe, previously used in the work of Montagna et al. (2016a) (Supplementary Table I,  Supplementary File I). See for details about the preservation and identification of specimens in Montagna et al. (2016a).
In order to perform the subsequent analyses, the collected beetles were clustered into the following four major groups according to the collection localities: i) C. flavipes from Alps and Apennines; ii) C. bameuli from the Alps (population from the Pyrenees was excluded due to low sampling and genetic monomorphism of individuals), iii) C. bameuli from the north of the Carpathians (Poland and Ukraine); and iv) C. bameuli from the south of the Carpathians -Pannonia. Such geographic grouping follow current knowledge on the phylogeography of beetles related with steppes and other xeric habitats of Europe (see Kajtoch et al. 2016) and is justified by phylogenetic trees presented in the preceding study (Montagna et al. 2016a).
In the niche modelling analyses we used geographic coordinates of the specimens collected during 2008-2014 field trips; in addition, the geographic coordinates of the specimens from Central, Eastern and South-Eastern Europe preserved in museum collections (Nature Museum of Institute of Systematics and Evolution of Animals PAS, Museum and Institute of Zoology PAS, Upper Silesian Museum in Bytom) were included in the analyses. The analyzed dataset includes specimens collected in 25 localities for C. flavipes, and in 32 localities for C. bameuli (Supplementary  Table II).

DNA extraction and PCR amplification
DNA isolation and sequencing of two genes (mitochondrial Cytochrome Oxidase I; cox1 and one nuclear Elongation Factor 1 alpha; ef1-α) were described in the earlier work of Montagna et al. (2016a).
In order to characterize the host plants of the two species, DNA isolates from the whole organism were used to amplify trnL(UAA) plant intron using A49325/B49863 primers (Taberlet et al. 1991), successfully used for leaf beetles (e.g. Montagna et al. 2013;De la Cadena et al. 2017). The choice of this marker was also dictated by previous development of barcode database for xeric plants of Central Europe , which was crucial for plant taxa assignment (see further chapter for details). PCR amplicons were obtained from the four geographical areas previously defined; the PCR amplicons were then mixed according to the region of collection and normalized, tagged and run on Illumina MiSeq 2 × 150 bp read lengths (using a MiSeq Reagent Kit v2).

Population genetic analysis
Standard genetic indices such as haplotype diversity (h), nucleotide diversity (π), and the number of private haplotypes (Np) for populations were computed with DnaSP v.5 (Librado & Rozas 2009). Ambiguous nucleotides in ef1-α gene sequences were resolved through the use of the Unphased/ Genotype Data option of DnaSP v.5 (Librado & Rozas 2009). The HAPAR (Wang & Xu 2003) algorithm was used to infer the haplotype sequences. These haplotypes were also inspected by eye and compared to other ef1-α sequences. Haplotype networks of cox1 and ef1-α sequences were constructed using the Median-Joining network method (Bandelt et al. 1999) in the Network 4.6.1.0 software (http:// www.fluxus-engineering.com/).

Demographic and migration analyses
A mismatch distribution (MD) (Rogers & Harpending 1992) on cox1 sequences was calculated in ARLEQUIN 3.5 (Excoffier & Lischer 2010) in order to examine demographic history, and specifically to test for historical (temporal) expansions of both species populations (in the case of C. bameuli, separately for each regional group of populations, as previously defined). The probable expansion time was estimated by the upper and lower limits of the τ parameter, adopting the Coleopteran evolutionary mean rate of 0.0177 (SD = 0.0019) for the mitochondrial cox1 (Papadopoulou et al. 2010) and using one generation per year.

Niche modelling
For environmental niche modelling (ENM) was used the MaxEnt algorithm (MaxEnt Model v3.3.3 k available at www.cs.princeton.edu/~scha pire/maxent; Phillips et al. 2006), with the following settings: number of iterations = 500, convergence threshold = 10 −5 , regularization multiplier = 1; maximum number of background points = 10 3 . The run was repeated 10 times, with cross-validation as a method for evaluation, and adding samples to the background to encompass all combinations of environmental data present in the data points. We used a mask in order to exclude certain areas from background points available for sampling that are beyond a known range of a species (Merow et al. 2013), confining the extent to 12°W, 66°N, 60°E, 30°N (excluding Iceland and areas beyond the Polar Circle). Area under the curve (AUC) of the receiver operating characteristic (ROC) plot was used as a measure of accuracy of the models produced.
For contemporary distribution models, we used 19 bioclimatic variables from the WorldClim data set (Hijmans et al. 2005), generated at 2.5-arc-minute resolution. The model was projected onto LGM data, reconstructed and made available by University of Miami -RSMAS (CCSM4), with the same variables and resolution. All climatic layers were downloaded from the website www.worldclim.org. In order to determinate the influence of the variables in the model were calculated the MaxEnt's heuristic estimates: the percentage contribution and the permutation importance, as well as the jackknife test performed for each of the variables separately (Phillips & Dudík 2008).

Host plant identification
The reads of trnL intron were identified by megablast analyses (E value cut-off = 1e-20; Altschul et al. 1990) and by comparison with the accurate trnL dataset obtained for dry grassland plants from central and eastern Europe . The methods used for this purpose (including thresholds of considered values) follow procedures adopted in Kajtoch et al. (2015). A read was considered to have a unique hit only if a single hit was reported or when the bitscore of the second-best hit was not greater than 0.95 x the bitscore of the best hit. When this condition was not met, then all hits (species) with the bitscores >0.95 x the bitscore of the best hit were considered as matching the read equally well. The obtained trnL illumina reads were identified at genus and species level (when possible) and then arranged in the operational taxonomic units (OTU) table in order to perform the comparative analyses. Diversity indices calculation and the following analyses (exceptions are specified) were performed using the software R package (R Project 3.0.2; http://cran.r-project.org/) package vegan (Dixon 2003). In order to test dissimilarities between the host plant communities associated with the four groups of specimens (i.e. C. bameuli from (i) the Alps, (ii) area north of the Carpathians, (iii) area south of the Carpathians; and (iv). C. flavipes from the Alps-Apennines), the dissimilarity matrices, obtained as previously described in Montagna et al. (2015aMontagna et al. ( , 2016b, were subjected to a nonparametric one-way analysis of similarity (ANOSIM; Clarke 1993) as in a previous work (Montagna et al. 2015b). To investigate the OTUs shared among the four analyzed groups, an analysis of commonality was performed and visualized through Venn diagrams using the g-plots package in R.

Population genetics and demographic history
The mitochondrial genetic diversity of the populations of C. flavipes and C. bameuli was found to be moderate or high, generally higher in the latter species (Table I); this observation could be the result of the unequal sample size. However, even in the case of an identical sample size (N = 9 for both), as for the Alpine populations, C. bameuli resulted to be more diverse than C. flavipes (mean Hdiv = 0.83-0.90 and 0.22, respectively). The most diverse populations of C. bameuli were localized in the Alps (mean Hdiv = 0.83-0.90), and in the north of the Carpathians (southern Poland) (mean Hdiv = 0.90). The populations from Pannonia (south of the Carpathians) (mean Hdiv = 0.60-0.70) and Ukraine (northeast of the Carpathians) (mean Hdiv = 0.53) were less diverse. No genetic diversity was found among the specimens from Pyrenees. On the other hand, with respect to the nuclear ef1-α gene, both species were found moderately diverse (mean Hdiv = 0.40-0.70), but some populations resulted monomorphic. In the single Carpathian locality was determined hybrid population of Cryptocephalus beetles. The inclusion of that population increases the values of genetic indices of all populations of C. flavipes, erroneously suggesting that this species has greater genetic diversity than C. bameuli. The Carpathian population was included in C. flavipes as its individuals have morphology resembling C. flavipes, despite their hybrid origin.
Haplotype networks constructed using cox1 and ef1-α haplotypes supported the distinctiveness of both the two species (Figure 1). In addition, the networks showed a higher diversity for C. bameuli populations, at least with respect to the cox1 marker, confirming previous results. There were interesting findings concerning the position of three specimens from the Carpathians that were morphologically identified as C. flavipes (identifiers F-KP1, F-KP2 and F-KP3). In the cox1 network, a single individual from the Carpathians (F-KP2) represented a distinct isolated lineage compared to other C. flavipes haplotypes from Italian Alps, whereas the haplotypes of F-KP1 and F-KP3 clustered Table I. Basic genetic diversity indices of cox-I and ef1-α markers calculated for C. flavipes and C. bameuli. nsample size, Vvariable sites, Hnumhaplotype number, Hdivhaplotype diversity with SE, лnucleotide diversity with SE, Snumber of segregating sites, Hprivnumber of private haplotypes.   Table III, Supplementary  Figure 1).
When tested in isolation (jackknife test), the variables of the highest gain were isothermality (0.972) for C. bameuli and precipitation of the wettest quarter (1.186) for C. flavipes. During the LGM, the extension of the C. bameuli predicted niche is broader, then that inferred for C. flavipes (Figure 2a of the predicted suitability ranges showed that under the current climatic conditions, C. bameuli areas of suitability (higher than 0.25 and higher than 0.5) are, respectively, about 2.3 and 1.8 times bigger than those of C. flavipes (Supplementary Table V). The inferred niche models for LGM showed similar proportions of the suitability areas. At Figure 2(c,d) are shown on maps the models of the current distribution of the species.
In terms of the current geographic distribution, C. flavipes resulted more confined to the mountain environments ( Figure 2c) while C. bameuli found to be widespread in various elevations (Figure 2d). The latter finds suitable living conditions also at lower elevations from western Europe to the Middle East and has its northern limit in the south of Scandinavian Peninsula (Figure 2d); the mountains of the northern Balkan Peninsula appeared to be suitable for C. bameuli but not for C. flavipes. Under LGM conditions, the inferred niches were less extended; suitable areas for C. bameuli were Jura Mountains, Po Valley, and the coast of the Black Sea, whereas C. flavipes was restricted to the Alpine region and to northern Anatolia ( Figure  2a,b).

Host plants
The analyses performed on the trnL Illumina reads confirmed the polyphagy of C. flavipes and C. bameuli. Indeed, they revealed 56 plant taxa (SD ± 11.1) on average associated with the three analyzed groups of C. bameuli (i.e. specimens collected from the Alps, north and south of Carpathians), while, a total of 58 plant taxa present in the diet of C. flavipes from the Alps-Apennines ( Figure 3). Interestingly, no significant differences were present in the diet composition of the three analysed groups of C. bameuli (ANOSIM analysis), while some differences were recovered when the diet composition of C. flavipes from the Alps-Apennines was included in the analysis (p-value = 0.042).
Regarding the taxonomical aspect of diet composition, several plant genera and species, mainly belonging to Rosaceae, Fabaceae, and Betulaceae, were identified in the diet of the regional populations of both species (Figure 3a). At the family level, no differences were observed among the groups of C. bameuli, whereas in the diet of C. flavipes were not present members of Asparagaceae, which instead resulted abundant in all groups of C. bameuli. Plant taxa shared among the three groups of C. bameuli and between C. flavipes and C. bameuli from the Alps were investigated using Venn diagrams (Figure 3b). As expected, the three groups of C. bameuli shared a core diet, encompassing 35 taxa of plants (e.g. Onobrychis viciifolia, Prunus spinosa, Rosa canina). Moreover, unique plant taxa were recovered in all of them, one in the populations from the Alps-Apennines, two and ten, respectively, in those inhabiting the areas in north and in south of the Carpathians. In detail, the analysed specimens of C. bameuli from Alps-Apennines mostly fed on plants of the genus Lotus and Onobrychis but also Rosa and Corylus; however, populations from both north and south of the Carpathians fed primarily on Crataegus, Prunus, Rosa, and on some herbs belonging to Rosaceae as well as on Anthericum, Asparagus, Genista and Onobrychis. The comparison between C. flavipes and C. bameuli inhabiting the same geographical area (i.e. Alps-Apennines) led to similar results, with 34 plant taxa (including, e.g. Corylus, Crataegus, Lotus and Onobrychis) shared among them out of a total of 68 (58 associated with C. flavipes and 44 with C. bameuli).

Different evolutionary histories
Population genetic analyses, niche modelling and the estimation of the demographic changes in the populations of C. flavipes and C. bameuli show that these two sibling species experienced rather different demographic histories and that these differences are not only species-specific but also geographically related. Interesting is that for both species similar sets of bioclimatic variables found to be crucial, that isothermality and precipitation. The same bioclimatic variables were also mostly responsible in niche modelling of other Cryptocephalus (e.g. C. barii, Brunetti et al. 2019), suggesting that the species of this genus respond to a similar set of climatic factors, however not necessarily in the same direction. So far, it was assumed that C. flavipes was widespread across the temperate zone of Europe and the Near East (Löbl & Smetana 2010); however, in some areas, its presence should be verified due to possible misidentifications of the two species. Our study does not allow for a definitive estimation of the C. flavipes range, but niche modelling suggests that it is restricted to the mountain systems across Europe; contrary to the much wider range reported by the Catalogue of Palaearctic Coleoptera (Löbl & Smetana 2010), our analyses showed that the species should be absent in the Balkans and in Anatolia. It is noteworthy that we found C. flavipes in a single locality in the Carpathians, despite the extensive collecting campaigns performed from 2008 to 2014 (~30 collecting campaigns). Even if Cryptocephalus species seem to persist at very low population densities (Piper 2002), our findings suggest that most of the C. flavipes specimens collected across central and eastern Europe and deposited in museum collections should be revisited. On the other hand, C. flavipes is widespread and abundant in the Alps and the Apennines. Demographic data estimated for the analysed Italian populations of this species suggest a possible expansion just in the middle of the Holocene period. The most probable explanation for the obtained pattern is that C. flavipes survived to the less favourable conditions of Pleistocene glaciations in refugia located approximately in the same areas where it is present now; this result is consistent with the predicted species distribution during LGM. The finding of phylogenetically divergent lineages in the Carpathians suggests the presence of a LGM refugium within the Carpathian Basin, which has been recently confirmed as an important extra-Mediterranean refugium (see Stewart & Lister 2001;Schmitt & Varga 2012). Niche modelling failed to locate this refugium, but it could be the consequence of a bias induced by the relatively low number of sampled sites in Central Europe. Conversely, genetic data and niche modelling analyses on C. bameuli suggested that this species experienced quite the opposite historical events. The fact that the subalpine populations of this species proved to be the most genetically diverse suggests the role of refugium of the Alps, especially of their western slope; this statement is also supported by the niche modelling. A further refugium could have existed in the Pyrenees; however, the population from this area proved to be genetically monomorphic. The observed genetic variability might be an effect of the low sample size (N = 5) or, most likely, of the high level of inbreeding of the population in which the specimens were collected. The western population of C. bameuli could have experienced an expansion between 170 Kya and 36 Kya. On the other hand, C. bameuli populations could have expanded to central and eastern Europe during the LGM or soon after it; this is consistent with the predicted species distribution during LGM (restricted to areas around the Alps and the Black Sea) and with present distribution (high probability of occurrence across hilly areas of Central, South-Eastern and Eastern Europe). The differences in the time of C. bameuli population expansions are not easily explained, but they might be linked to the different environmental and climatic conditions present in the subalpine areas and around the Carpathians. The glaciated areas during LGM were more extensive in central and eastern than in western (especially southwestern) Europe. The current south Alpine piedmont was generally ice-free, whereas an ice sheet covered most of the areas north of the Carpathians (Ehlers & Gibbard 2004). This assertion, however, does not explain the similar and relatively late expansion of Pannonian populations, but the observed patterns for populations north and south of the Carpathians suggest that C. bameuli could have expanded from some eastern refugia (e.g. Black Sea coast, as it was predicted in niche modelling) and colonized both sides of the Carpathians in the same time. The estimates of the time of expansion for the populations of C. bameuli strongly suggest that this species is much more coldtolerant than C. flavipes. Climate warming could lead to the further expansion of C. flavipes and to the retreat of C. bameuli, at least from the southernmost localities. The latter scenario was also recently predicted for cold-adapted mountain C. barii in the Alps (Brunetti et al. 2019).
The different stories experienced by the two species suggest that despite their similar morphologies, they are probably adapted to different climatic and environmental conditions (including some differences in food preferences, as seen below). We hypothesize that C. flavipes might be associated with mountain Mediterranean xeric habitats on limestone, whereas C. bameuli is probably more related to continental steppes. Differently to the majority of steppic beetles (see Kajtoch et al. 2016 for summary), the populations of C. bameuli from the inner and outer side of the Carpathians express a high level of genetic similarity. The haplotypes in common among C. bameuli populations living north and south of the Carpathians could be a consequence of (i) a recent expansion of this species in central and eastern Europe, as supported by demographic analyses and niche modelling and (ii) current high migration rates. The limited dispersal abilities of Cryptocephalus beetles (e.g. as demonstrated for Cryptocephalus decemmaculatus and Cryptocephalus nitidulus; Piper & Compton 2010) support the first scenario. On the other hand, shared haplotypes were detected only in adjacent areas on both sides of the Carpathians (Moravia and southern Poland), suggesting a limited bidirectional dispersal of individuals in current times. However, both explanations are not exclusive since the recent species expansion across the Carpathians could have been continued via, e.g., the Moravia Gate, a large pass between the Carpathians and the Sudetes known as one of the major migration routes for xeric species in central Europe (Mazur 2001).

Possible hybridization
The only population of central Europe in which were present specimens morphologically identified as C. flavipes was located in the Carpathians (Pieniny Mountains), well-known local refugium with the presence of xeric habitats. Interestingly, the only three examined specimens of this population expressed mixed genotypes. Despite a C. flavipes phenotype, only one specimen harboured a C. flavipesrelated mtDNA, whereas the other two had C. bameulirelated mtDNA. Moreover, two of these specimens showed intermediate variants of their nuclear DNA (ef1-α gene). These findings strongly suggest that this isolated population of C. flavipes had a substantial gene flow with C. bameuli of central Europe. It has been hypothesized that when one of the two closely related species becomes rare, it could mate and hybridize with its more abundant relative (Montagna et al. 2016a). Moreover, it has been postulated that hybridization is not a rare phenomenon among Cryptocephalus beetles (e.g. the Cryptocephalus sericeus complex; Gómez-Zurita et al. 2012). The alternative scenario, which is in our opinion much less plausible, is that specimens from the Carpathians harbour ancestral polymorphisms and are suffering from incomplete lineage sorting. Further studies using larger samples, including presumable hybrid populations and using next-generation sequencing technologies, are necessary to be able to disentangle among these scenarios.

Host plants
Until now, C. flavipes has been considered a polyphagous species since specimens were observed and collected on different plant species, mainly deciduous trees and bushes (e.g. as Corylus, Crataegus, Quercus, Crataegus, Salix; Müller 1949Müller -1953Burlini 1955;Fogato & Leonardi 1980;Burakowski et al. 1990); however, no conclusive evidence have been provided. With regard to the recently described species C. bameuli, the host plants were unknown and the existence of possible variation in feeding preferences, with respect to the different geographic areas, has never been investigated.
Metabarcoding analysis on the gut content of C. bameuli and C. flavipes confirmed their polyphagy; the most commonly exploited host plants belong to the families of Betulaceae, Fabaceae and Rosaceae. The same analysis showed that the three groups of C. bameuli share a core diet, but some differences were observed when populations from different regions were compared, especially between the Alpine and central-eastern European populations. Such slight regional differences in the diet of polyphagous beetles have also been reported for the weevil Centricnemus leucogrammus Germ. (Kajtoch 2014), and it may reflect some local feeding adaptations to different sets of plants, even if the core diet of the species is maintained. Considering the overall populations of both species, significant differences in diet composition were observed. It is noteworthy that the Alpine populations of C. bameuli and C. flavipes, inhabiting the same habitats and geographical area, share some host plants but also possess a higher number of unique hosts. Further field or lab-based observations could clarify the exact feeding preferences of these beetles.

Conclusions
The presented results provide evidence that these two sibling and sympatric taxa have experienced different histories along their evolutionary path, which may reflect their different adaptations to climatic and environmental conditions. Moreover, both species have differentiated their ecological niche (e.g. living in different xeric habitats at different elevations) and, possibly, have undergone ecological speciation through switches in their host plants. Further studies, including denser sampling as well as the use of different types of genetic markers, could help in estimating the magnitude of hybridization (frequency and geographic span) and in identifying the biological and ecological adaptations that have differentiated these sibling beetles.