BUSCADOR RECOLECTA

[Aim] There is little consensus about the relative roles of biotic versus abiotic factors in setting limits to species distributions or in generating geographical patterns of species richness. However, despite the probable importance of host availability in governing the distribution and diversity of consumers, few studies have simultaneously tested the effects of resource distribution and diversity on consumer ranges and richness patterns., [Location] Sierra de Guadarrama, central Spain., [Methods] We examined the effects of biotic resources, consumer attributes and climate on the ranges and species richness patterns of 43 specialist butterflies at 40 sites over a 1800-m elevational gradient. Evidence for resource use was based on comprehensive field records of oviposition and larval feeding on host plants., [Results] We show that limitation by either biotic interactions with resources (the distributions and parts eaten of the larval host plants) or intrinsic dispersal ability was stronger at upper than lower elevational range limits for butterflies. Both resource and consumer richness followed a unimodal, humped pattern over the elevational gradient, but host plant richness peaked 300 m lower than butterfly richness. In addition, whereas changes in butterfly species richness were roughly symmetrical around peak richness over the gradient studied, the host plants showed markedly lower species richness at high elevations (> 1750 m). Butterfly species richness increased with host plant resource diversity and relative humidity, with a steeper response to host plant richness in cooler sites (at higher elevations)., [Main conclusions] The results demonstrate the role of bottom-up control by resource availability in limiting consumer distributions and richness. Importantly, resource limitation had increasing relevance towards the coolest parts of environmental gradients and those poorest in resource species, with potential consequences for ecological responses to environmental change., The research was funded by Universidad Rey Juan Carlos/Comunidad de Madrid (URJC-CM-2006-CET-0592), the Spanish Ministry of Economy and Competitiveness (REN2002-12853-E/GLO, CGL2005-06820/BOS, CGL2008-04950/BOS, CGL2011-30259, CGL2013-48277-P and CGL2014-57784-P), the British Ecological Society and the Royal Society. Research permits were provided by Comunidad de Madrid, Parque Nacional de Guadarrama, Parque Natural de Peñalara, Parque Regional de la Cuenca Alta del Manzanares, Parque Regional del Curso Medio del Río Guadarrama, Patrimonio Nacional and Ayuntamiento de Cercedilla., Peer reviewed

1. Phenological change is the most widely documented biological impact of climate change, but shows marked variation in magnitude among populations and species. Thus, quantifying the environmental factors and organismal differences driving this intra- and interspecific variability in phenology is vital to understand and forecast the ecological consequences of climate change.
2. Here, we test intra- and interspecific differences for a set of butterfly species in the organismal sensitivity of flight phenology and its dependence on environmental factors, using as our model system an elevation gradient in a Mediterranean mountain range where temperature and relative humidity vary substantially over space and time.
3. We use field-collected meteorological data, and butterfly counts for 20 univoltine species over 14 years, to test the relative effects on phenology of temperature and relative humidity, the sensitivity of phenology to spatial and temporal variation in temperature and whether ecological traits account for inter-specific variation
in sensitivity.
4. For all species, temperature in the months immediately preceding adult emergence had the strongest relationship with phenology. All species appeared earlier in warmer years, with those flying earlier in the season showing the greatest sensitivity to annual (temporal) variation in temperature. However, only a minority
of species showed evidence of plastic, space-for-time responses to temperature. Instead, most species showed strong evidence that phenology was more sensitive to temporal than spatial variation in temperature.
5. Our results support the dominant influence of temperature on phenology, even in Mediterranean environments suffering summer drought. They also suggest that accurate forecasts of species' phenological shifts could require the isolation of spatial from temporal components of temperature variation, because the sensitivity of populations and species may differ across these two dimensions. The factors driving synchronisation of phenology over space merit particular research in the context of climate change, given their potential to expose populations simultaneously to environmental extremes., British Ecological Society; Ministerio de Economía y Competitividad. Grant Numbers: CGL2005‐06820/BOS, CGL2008‐04950/BOS, CGL2011‐30259, CGL2014‐57784‐P, REN2002‐12853‐E/GLO; Universidad Rey Juan Carlos / Comunidad de Madrid. Grant Numbers: PEJD‐2017‐PRE/AMB‐4075, URJC‐CM‐2006‐CET‐0592; Royal Society. Grant Number: JP10052., Peer reviewed

Parasitism is a key factor in the population dynamics of many herbivorous insects, although its impact on host populations varies widely, for instance, along latitudinal and altitudinal gradients. Understanding the sources of geographical variation in host–parasitoid interactions is crucial for reliably predicting the future success of the interacting species under a context of global change.
Here, we examine larval parasitism in the butterfly Aglais urticae in south-west Europe, where it is a mountain specialist. Larval nests were sampled over 2 years along altitudinal gradients in three Iberian mountain ranges, including the Sierra Nevada, home to its southernmost European population. Additional data on nettle condition and adult butterflies were obtained in the study areas.
These data sources were used to investigate whether or not differences in parasitism rates are related to the geographical position and phenology of the host, and to the availability of the host plants.
Phenological differences in the host populations between regions were related to the severity of summer drought and the corresponding differences in host plant availability. At the trailing-edge of its distribution, the butterfly's breeding season was restricted to the end of winter and spring, while in its northern Iberian range the season was prolonged until mid-summer. Although parasitism was an important source of mortality in all regions, parasitism rates and parasitoid richness were highest in the north and lowest in the south. Moreover, within a region, there was a notable increase in parasitism rates over time, which probably led to selection against an additional late summer host generation in northern regions. Conversely, the shorter breeding season in Sierra Nevada resulted in a loss of synchrony between the host and one important late season parasitoid, Sturmia bella, which may partly explain the high density of this butterfly species at the trailing-edge of its range.
Our results support the key role of host phenology in accounting for differences in parasitism rates between populations. They also provide insights into how climate through host plant availability affects host phenology and, ultimately, the impact of parasitism on host populations., Sierra Nevada fieldwork was partly supported by the Global Change Observatory of Sierra Nevada that is funded by the Andalusian Government (Ministry of Agriculture, Livestock, Fisheries and Sustainable Development) with the collaboration of the University of Granada. Funding for this research was provided by Ministerio de Ciencia e Innovación, R + D Programa Nacional, Proyecto I + D + I (grant CGL2014-57784-P to D.G.)., Peer reviewed

3 tables. -- Stefanescu_et_al_JAE_data_butterfly_counts.csv: overwintering adult, fresh adult and larval nest counts of Aglais urticae for the three study regions. Counts were made fortnightly (15 visits, conditions permitting). -- Stefanescu_et_al_JAE_data_nettle_phenology.csv: nettle height (cm) and quality (1, 2, 3, 4) for the three study regions. Data were recorded fortnightly (16 visits per site, conditions permitting). -- Stefanescu_et_al_JAE_data_larval_parasitism.csv: larval nests collected for estimating parasitism rates for the three study regions. -- All models were run in R (R Core Team, 2018)., Parasitism is a key factor in the population dynamics of many herbivorous insects, although its impact on host populations varies widely, for instance, along latitudinal and altitudinal gradients. Understanding the sources of geographical variation in host-parasitoid interactions is crucial for reliably predicting the future success of the interacting species under a context of global change. Here, we examine larval parasitism in the butterfly Aglais urticae in south-west Europe, where it is a mountain specialist. Larval nests were sampled over two years along altitudinal gradients in three Iberian mountain ranges, including the Sierra Nevada, home to its southernmost European population. Additional data on nettle condition and adult butterflies were obtained in the study areas. These data sources were used to investigate whether or not differences in parasitism rates are related to the geographical position and phenology of the host, and to the availability of the host plants. Phenological differences in the host populations between regions were related to the severity of summer drought and the corresponding differences in host plant availability. At the trailing-edge of its distribution, the butterfly’s breeding season was restricted to the end of winter and spring, while in its northern Iberian range the season was prolonged until mid-summer. Although parasitism was an important source of mortality in all regions, parasitism rates and parasitoid richness were highest in the north and lowest in the south. Moreover, within a region, there was a notable increase in parasitism rates over time, which probably led to selection against an additional late-summer host generation in northern regions. Conversely, the shorter breeding season in Sierra Nevada resulted in a loss of synchrony between the host and one important late-season parasitoid, Sturmia bella, which may partly explain the high density of this butterfly species at the trailing-edge of its range. Our results support the key role of host phenology in accounting for differences in parasitism rates between populations. They also provide insights into how climate through host plant availability affects host phenology and, ultimately, the impact of parasitism on host populations., [Study system] We studied the complex of larval parasitoids of Aglais urticae in three regions encompassing the whole of its latitudinal range in the Iberian Peninsula: the Pyrenees in the north of the Iberian Peninsula, the Sierra de Guadarrama in central Spain, and Sierra Nevada in southern Spain. Sampling sites were established along an altitudinal gradient in each region that covered most of the altitudinal range in which this species breeds. Nine sites were sampled for parasitoids in the Pyrenees at 1,127–2,560 m a.s.l., seven sites in the Sierra de Guadarrama at 1,150–2,004 m a.s.l., and six sites in Sierra Nevada at 975–2,532 m a.s.l. As part of a larger project aimed at investigating various aspects of the ecology of A. urticae, additional sites were surveyed in each region. The information gathered at these additional sites was used in this work to improve knowledge of the phenology of this butterfly., [Field sampling] To study the phenology of A. urticae adults, we used 500-m transects on which butterflies were counted every two weeks from March to September (a total of 15 sampling visits), following the standard methodology of the Pollard walks (Pollard & Yates, 1993). Butterflies were classified either as overwintered or freshly emerged based on wing colouration (i.e. dull or brightly coloured, respectively). Transects were walked at 15, 24 and 14 sites in 2016, and at 16, 24 and 20 sites in 2017, in the Pyrenees, Sierra de Guadarrama and Sierra Nevada, respectively. To study the phenology of A. urticae larvae, we counted all larval nests found in four (exceptionally, just two and three at two sites) focal U. dioica patches in a subsample of the sites used for adult counts in each region. The focal patches were randomly selected along the butterfly transects and, if not available, at other accessible sites that were as close as possible to the transect route. The focal patches were visited every two weeks from March to September, whenever possible during the same visits as for the adult transect counts. To study larval parasitism, larval nests detected at focal U. dioica patches (see above) were marked and, if larvae were in the third or later instars, 20 individuals were collected to assess parasitism. Because the total number of larvae per nest was sometimes less than 20, the overall average number of larvae per sample (± SD) was 16.2 ± 6.9. Moreover, given that the opportunistic parasitoid, Cotesia vestalis, is known to parasitize first instar larvae of the small tortoiseshell and to emerge mainly from the second instar (Audusseau et al., 2021), in 2017 we also collected 8 samples of five second instar larvae in all three regions (3 samples in the Pyrenees, 2 in the Sierra de Guadarrama and 3 in Sierra Nevada). We did not assess pupal parasitism, even if it may be important (Pyörnilä,1977; Shaw et al., 2009) because pupae are difficult to locate in the field, thereby precluding any reliable estimates of mortality. Larvae were reared indoors in transparent plastic containers (155 x 105 x 45 mm) in groups of up to five individuals, all from the same sample. To avoid possible contamination, larvae were always reared with nettle leaves collected from their original nettle patch; if not available, nettles were harvested from sites where A. urticae and its closest congener, A. io, were absent, since some common parasitoids (e.g. the tachinids Sturmia bella and Pales pavida) lay microtype eggs on nettle leaves that can infect caterpillars if they eat these leaves. When a larva or pupa (in the case of larva-pupal parasitoids) produced a parasitoid, we recorded the stage at which the host was killed and kept the parasitoid individually in a vial until the adult emerged. Adults were preserved in pure ethanol for identification (Ichneumonoidea by MRS, Tachinidae by DH). Although hatching success was generally poor for most tachinids, careful inspection of puparia allowed for correct identification in almost all cases. To examine host plant availability, we recorded the growing condition (i.e. quality level and height) of nettles over the season at a subsample of sites used for larval nest counts in each region. At each visit, two stems were randomly selected from each nettle patch. Their height was measured (in cm) and they were given a categorical value from 1 (worst quality) to 4 (best quality) in which (1) corresponds to already dry or withered plants, with senescent leaves; (2) to flowering plants and plants with green but not fresh leaves; (3) to old plants in which regrowth leaves were beginning to become visible (a common situation at the end of summer after rain or after herbivory); and (4) to vigorous plants, with fresh leaves. Category ranking was based on previous work showing how nettles in these various phenological stages differently affect larval growth rates, pupal and adult weights in the small tortoiseshell (Pullin, 1987) and the map butterfly, Araschnia levana (Mevi-Schütz & Erhardt, 2005)., [Host phenology] A combination of the standardized adult and larval count data was used to define the phenology of the species. GAM models were fitted to the adult (overwintering and fresh butterflies separately) and larval nest counts, which allowed us to extract the Julian day corresponding to each peak of abundance in a given region and season. GAM models were built using the package mgcv in R (Wood, 2011). In these models we used pooled data from 2016 and 2017 to increase the sample size and to improve the overall phenological picture in each of the study regions. To investigate the potential altitudinal delay in larval phenology, we regressed separately the timing of larval nest appearance against site altitude for each year and region. The timing of larval nest appearance was summarised as the weighted mean appearance date., [Host plant phenology] We tested for differences in the phenology of nettles between regions using GAMM models, in which either nettle condition or height were the response variables and altitude, region, year, visit number (i.e. the timing of the season, used as the smoothing term) and the interactions of region with both visit number and altitude were the predictors. In these models, each individual stem was used as a data point and 'nettle patch' was entered as a random factor., [Impact of parasitism on host populations] To test whether or not the number of parasitoid species was comparable between regions (because sample sizes differed greatly between regions, see below), we computed the most common nonparametric estimators of species richness for each region and year separately (based on all recorded parasitoid species and genera) using the SpadeR package in R (Chao & Chiu, 2016). To assess the impact of parasitoids on host populations, we calculated the parasitism rate for each larval nest as the number of larvae killed as a result of parasitism, after discounting those that died for unknown reasons (i.e. our calculations were always based on effective larval samples). To avoid biases resulting from low sample sizes, the parasitism rate was calculated for effective samples of ≥ 5 larvae. We obtained very similar results (not shown) when models of parasitism rate were built following a more restrictive criterion of effective samples of ≥ 10 larva. The parasitism rate was modelled with generalized linear models (GLMs) using a binomial distribution and logit link function, with region, Julian day (date of nest collection) and altitude as predictors. However, because altitude and Julian day were highly correlated (r=0.76), only models with just one of these two variables were retained in the end. Models were built separately for 2016 and 2017 because the sampling sites in Andorra differed slightly between the two years. All possible models were built using the package lme4 in R (Bates et al., 2015); the best models were selected using the MuMIn package in R (Barton, 2015), with model selection being based on the Akaike Information Criterion (AIC). Models that differed by < 2 points from the lowest AIC (∆AIC < 2) were considered the top-ranked models (statistically equivalent to the best model of the set)., [Phenological overlap between the host and its main parasitoids] The overlap (i.e. temporal co-occurrence) between the host and its two main parasitoids, Pelatachina tibialis and Sturmia bella, was estimated using the Overlap Parasitoid-Host index (OPH), as described by Audusseau et al. (2020), which is bouinbded from zero to one. This index was calculated for all possible combinations of site and year in each region. The maximum value of 1 is obtained in the hypothetical case when all individuals recorded in a given season, both of the parasitoid and the host, are concentrated in the same sampling event k. The opposite situation occurs (index value equaling to zero) when in all available samples one of the interacting species is always missing.

To understand which factors explain the degree of overlap between the parasitoid and the host, we used GLM models with a binomial distribution and a logit-link function. For each parasitoid species and year we applied a model in which OPH was the dependent variable and site altitude and region were the variable predictors., Ministerio de Ciencia e Innovación, R+D Programa Nacional, Proyecto I+D+I , Award: CGL2014-57784-P, Stefanescu_et_al_JAE_data_butterfly_counts.csv; Stefanescu_et_al_JAE_data_nettle_phenology.csv; Stefanescu_et_al_JAE_data_larval_parasitism.csv, Peer reviewed

Glacial cycles lead to periodic population interbreeding and isolation in warm-adapted species, which impact genetic structure and evolution. However, the effects of these processes on highly mobile and more cold-tolerant species are not well understood. This study aims to shed light on the phylogeographic history of Aglais urticae, a butterfly species with considerable dispersal ability, and a wide Palearctic distribution reaching the Arctic. Through the analysis of genomic data, four main genetic lineages are identified: European, Sierra Nevada, Sicily/Calabria/Peloponnese, and Eastern. The results indicate that the Sardo-Corsican endemic taxon ichnusa is a distinct species. The split between the relict lineages in southern Europe and the main European lineage is estimated to have happened 400¿450 thousand years ago, with admixture observed during the Quaternary glacial cycles, and still ongoing, albeit to a much smaller extent. These results suggest that these lineages may be better treated as subspecific parapatric taxa. Ecological niche modelling supported the existence of both Mediterranean and extra-Mediterranean refugia during the glacial periods, with the main one located on the Atlantic coast. Nevertheless, gene flow between populations was possible, indicating that both differentiation and admixture have acted continuously across glacial cycles in this cold-tolerant butterfly, generally balancing each other but producing differentiated lineages in the southern peninsulas. We conclude that the population dynamics and the processes shaping the population genetic structure of cold-adapted species during the Quaternary ice ages may be different than those classically accepted for warm-adapted species., This research was supported by grants PID2019-107078GB-I00 and PID2022-139689NB-I00 (funded by Ministerio de Ciencia e Innovación (MCIN)/Agencia Estatal de Investigación (AEI)/10.13039/501100011033 and ERDF A way of making Europe) and by grant 2021-SGR-00420 from Departament de Recerca i Universitats de la Generalitat de Catalunya to RV, by grant CGL2014-57784-P funded by MCIN/AEI/10.13039/501100011033 to DG, by PRE2020-094870 funded by MCIN/AEI/10.13039/501100011033 and ‘European Social Fund (ESF) Investing in your future’ to VM, and by BES-2017-080641 funded by MCIN/AEI/10.13039/501100011033 and by ‘ESF Investing in your future’ to JCH. L.D. acknowledges the support of NBFC to University of Florence, Department of Biology, funded by the Italian Ministry of University and Research, PNRR, Missione 4 Componente 2, ‘Dalla ricerca all'impresa’, Investimento 1.4, Project CN00000033 and by several projects with the Direttiva Biodiversità of the Italian Environmental Ministry. G.T. was supported by grant PID2020-117739GA-100 funded by Ministerio de Ciencia e Innovación (MCIN)/Agencia Estatal de Investigación (AEI)/10.13039/501100011033 and by grant 2021-SGR-01334 from Departament de Recerca i Universitats de la Generalitat de Catalunya., Abstract
1 INTRODUCTION
2 MATERIALS AND METHODS
2.1 Sampling
2.2 ddRADseq library preparation
2.3 Data processing
2.4 Population genetic structure
2.5 Phylogenetic inference and gene flow analyses
2.6 Genetic diversity analyses
2.7 Demographic history
2.8 Ecological niche modelling (ENM)
2.9 Wolbachia infection analysis
3 RESULTS
3.1 Principal component analysis
3.2 Genetic structuring and hybridization between lineages
3.3 Phylogenetic inference
3.4 Gene flow analyses
3.5 Genetic diversity analyses
3.6 Demographic history
3.7 Ecological niche modelling
3.8 Wolbachia infection analysis
4 DISCUSSION
4.1 Genomic data clarify previous mitochondrial results
4.2 Evidence of admixture
4.3 Genetic diversity
4.4 The opposed forces of differentiation by isolation and admixture by gene flow
4.5 A complex scenario of differentiation with gene flow
AUTHOR CONTRIBUTIONS
ACKNOWLEDGEMENTS
CONFLICT OF INTEREST STATEMENT, Peer reviewed

Au_filtered_PCA.vcf - used as input for PCA, Au_STRUCTURE.str - used as input for STRUCTURE analysis, Au_IQTREE.phy - used as input for phylogenetic inference in IQTREE, Au_PhyloNetworks.loci - used as input to create gene trees for subsequent use in PhyloNetworks analysis, Au_filtered_TreeMix.vcf - used as input for TreeMix analysis, Au_TASSEL_western.vcf; Au_TASSEL_eastern.vcf; Au_TASSEL_central.vcf - used as input for genetic diversity analysis in TASSEL, Au_filtered_dadi.vcf - used to produce SFS file for dadi analysis, Au_EuSicily_hybridsim; Au_EuSierra_hybridsim - used to simulate hybrids with HYBRIDLAB, Peer reviewed

Au_filtered_dadi.vcf - used to produce SFS file for dadi analysis, Au_EuSicily_hybridsim; Au_EuSierra_hybridsim - used to simulate hybrids with HYBRIDLAB, Peer reviewed

Au_filtered_PCA.vcf - used as input for PCA, Au_STRUCTURE.str - used as input for STRUCTURE analysis, Au_IQTREE.phy - used as input for phylogenetic inference in IQTREE, Au_PhyloNetworks.loci - used as input to create gene trees for subsequent use in PhyloNetworks analysis, Au_filtered_TreeMix.vcf - used as input for TreeMix analysis, Au_TASSEL_western.vcf; Au_TASSEL_eastern.vcf; Au_TASSEL_central.vcf - used as input for genetic diversity analysis in TASSEL, Peer reviewed