The evolution of climate tolerance in conifer‐feeding aphids in relation to their host's climatic niche

Abstract Climate adaptation has major consequences in the evolution and ecology of all living organisms. Though phytophagous insects are an important component of Earth's biodiversity, there are few studies investigating the evolution of their climatic preferences. This lack of research is probably because their evolutionary ecology is thought to be primarily driven by their interactions with their host plants. Here, we use a robust phylogenetic framework and species‐level distribution data for the conifer‐feeding aphid genus Cinara to investigate the role of climatic adaptation in the diversity and distribution patterns of these host‐specialized insects. Insect climate niches were reconstructed at a macroevolutionary scale, highlighting that climate niche tolerance is evolutionarily labile, with closely related species exhibiting strong climatic disparities. This result may suggest repeated climate niche differentiation during the evolutionary diversification of Cinara. Alternatively, it may merely reflect the use of host plants that occur in disparate climatic zones, and thus, in reality the aphid species' fundamental climate niches may actually be similar but broad. Comparisons of the aphids' current climate niches with those of their hosts show that most Cinara species occupy the full range of the climatic tolerance exhibited by their set of host plants, corroborating the hypothesis that the observed disparity in Cinara species' climate niches can simply mirror that of their hosts. However, 29% of the studied species only occupy a subset of their hosts' climatic zone, suggesting that some aphid species do indeed have their own climatic limitations. Our results suggest that in host‐specialized phytophagous insects, host associations cannot always adequately describe insect niches and abiotic factors must be taken into account.

reconstructed at a macroevolutionary scale, highlighting that climate niche tolerance is evolutionarily labile, with closely related species exhibiting strong climatic disparities. This result may suggest repeated climate niche differentiation during the evolutionary diversification of Cinara. Alternatively, it may merely reflect the use of host plants that occur in disparate climatic zones, and thus, in reality the aphid species' fundamental climate niches may actually be similar but broad. Comparisons of the aphids' current climate niches with those of their hosts show that most Cinara species occupy the full range of the climatic tolerance exhibited by their set of host plants, corroborating the hypothesis that the observed disparity in Cinara species' climate niches can simply mirror that of their hosts. However, 29% of the studied species only occupy a subset of their hosts' climatic zone, suggesting that some aphid species do indeed have their own climatic limitations. Our results suggest that in host-specialized phytophagous insects, host associations cannot always adequately describe insect niches and abiotic factors must be taken into account.

K E Y W O R D S
climate, insect-plant interactions, niche equivalency tests, niche evolution, phylogeny, phytophagous insects

| INTRODUC TI ON
The niche theory is a central concept in ecology and evolution (Duran & Pie, 2015;Pearman, Guisan, Broennimann, & Randin, 2008;Pearman et al., 2014;Petitpierre et al., 2012;Wiens et al., 2010), with the set of ecological conditions where species occur (Hutchinson, 1957) considered one of the key factors that drive diversification and biogeographic patterns on Earth. With more than 500,000 known species, phytophagous insects represent nearly a quarter of the whole terrestrial macroscopic biodiversity (Daly, Doyen, & Purcell, 1998;Southwood & Norton, 1973;Strong, Lawton, & Southwood, 1984). However, the factors that shape the geographic distribution of phytophagous insects are still poorly understood, limiting our ability to model the effects of climate change on their future distributions; there is also an additional complexity of combining the insects' and host plants' climate niches (Barredo et al., 2015).
Most phytophagous insect species are host-specific, feeding on one or a few plant species (Jaenike, 1990), and thus, their ecology and life history are strongly dependent on their interaction with their host plants. Consequently, the study of phytophagous insect niches is often restricted to the traits characterizing their association with their host plants: the range of plant species they feed on (Kaplan & Denno, 2007), the plant organs on which they develop (Condon & Steck, 1997;Cook, Rokas, Pagel, & Stone, 2002), or the range of natural enemies encountered on these plants (Singer & Stireman, 2005). This traditional host-centered conception of the ecological niche of phytophagous insects implies that their distribution is primarily constrained by the geographic range of their hosts, whereas other environmental abiotic factors, such as geographic or climatic barriers, are rarely taken into account (Futuyma & Moreno, 1988;Jaenike, 1990). Similarly, insect diversification and biogeographic history are also often interpreted as the consequence of host shifts (von Dohlen, Kurosu, & Aoki, 2002;Janz, Nylin, & Wahlberg, 2006;Nyman, Linder, Peña, Malm, & Wahlberg, 2012).
In the last decades, the study of the environmental conditions relevant to large-scale ecological and geographic properties of species has been conducted in numerous organisms such as amphibians (Bonetti & Wiens, 2014;Kozak & Wiens, 2010), squamates (Knouft, Losos, Glor, & Kolbe, 2006;Pyron & Burbrink, 2009), plants (Evans, Smith, Flynn, & Donoghue, 2008;Meseguer, Lobo, Ree, Beerling, & Sanmartín, 2014), mammals (Dormann, Gruber, Winter, & Herrmann, 2009), and birds (Pearman et al., 2014). These reveal an important role for climate adaptation in the evolution of these organisms. The climate niche reflects the physiological tolerance to climatic conditions (Grinnell, 1917) and may strongly influence the distribution of species over space and time (Soberón, 2007). Unfortunately, climatic niche evolution has been scarcely investigated in phytophagous insects (but see Kellermann et al., 2012;Pitteloud et al., 2017) and we still have a limited understanding of the distribution in climate space of related insect species. This relative neglect may be due to the difficulty of gathering climatic distribution data for an entire clade of phytophagous insects, but it may also stem from the assumption that their distributions simply reflect host plant distributions and are therefore uninteresting in their own right.
Phytophagous insect development, survival, and capacity to exploit plants are known to be sensitive to climate conditions, and insects do exhibit many climate dependent traits (Bale et al., 2002;Gallego, Verdú, Carrascal, & Lobo, 2016;Howe, 1967;Régnière, Powell, Bentz, & Nealis, 2012). For example, Duyck, David, and Quilici (2006) showed that Ceratitis fruit fly pupae exhibit different climatic niches, coexisting on the French island of La Réunion in the southwest Indian Ocean under different humidity conditions. Another telling example is the hemlock woolly adelgid (Adelges tsugae), an invasive insect introduced to eastern North America that is decimating populations of eastern and Carolina hemlock (Tsuga canadensis and T. caroliniana) from Georgia to Connecticut. The success of its invasion is limited by its low cold tolerance: several T. canadensis and T. caroliniana populations distributed in colder regions have not been impacted by the pest (Paradis, Elkinton, Hayhoe, & Buonaccorsi, 2008). It is also well known that several insect species are tracking cooler habitats or hatching earlier than they used to in response to climate warming (Visser & Holleman, 2001). Finally, studies have shown that some phytophagous insects occupy only a part of the range of their host plants (Battisti et al., 2005;Carroll, Taylor, Régnière, & Safranyik, 2003;Godefroid et al., 2016;Hill & Hodkinson, 1992), suggesting that climate variables, such as temperature and precipitation, may play a significant role in defining the insects' habitat and distribution. Altogether, these studies provide evidence that phytophagous insects can display adaptations to specific climate conditions.
Here, we investigated the role of climate adaptation in the longterm evolution of a clade of aphids (Hemiptera; Aphididae) of the genus Cinara, in relation to their host plants. We sought to answer the following general questions: (a) Does this phytophagous insect lineage exhibit long-term climatic specialization? (b) Do aphids occur in the full range of climate conditions occupied by their host plants, or do they only occupy a subset of the host's climatic niche?
Concordance between the insects' and hosts' climatic niches would suggest that observed insect climatic niches are limited by the hosts' presence rather than their own set of physiological constraints.
Aphids are particularly well-suited to address these questions as they are generally host-specific, they spend their whole life cycle on their hosts, and their range of host preference and their geographic distribution are generally well-known due to the pest status of numerous species. Using a georeferenced occurrence dataset and a comprehensive phylogeny, we studied the evolution of climate tolerance of the conifer-feeding aphid genus Cinara Curtis, 1835 (Aphididae: Hemiptera) in relation to their Pinus spp. hosts.

| Study system, distribution and bioclimatic data
The genus Cinara (Aphididae: Lachninae) includes 253 described species, including those of the recently subordinated subgenus Schizolachnus (Chen, Favret, Jiang, Wang, & Qiao, 2016;Meseguer, Coeur d'acier, Genson, & Jousselin, 2015). Cinara species have their native range exclusively within the Northern Hemisphere, and they are generally host-specific, feeding on one or a few plant species of the conifer families Pinaceae and Cupressaceae (Blackman & Eastop, 2000;Favret & Voegtlin, 2004). Associations with host plants are well known for most Cinara species as they have been the subject of many taxonomic investigations (Eastop, 1972;Favret & Voegtlin, 2004;Voegtlin, 1976), and these association have been reviewed and compiled in the reference book by Blackman and Eastop (1994) which is regularly updated (Blackman & Eastop, 2019). We compiled a total of 6,135 Cinara occurrence records for 106 species. These include data from online databases (GBIF with 2,359 records and IBOL with 1,334 records) and aphid taxonomic databases from the entomological collection of the University of Montreal (1,705 records), the Institute of Zoology of the Chinese Academy of Sciences (221 records), and the INRA collection (516 records). All occurrences were reviewed and ambiguous accessions (e.g., occurrences for which the location names did not correspond with the given GPS coordinates or were assigned to county centroids) and duplicate records were removed. Occurrence data for Pinus species hosting Cinara were retrieved from the Conifers of the World database (Farjon & Filer, 2013).
In order to estimate the aphid climatic niches and those of their associated Pinus hosts, we selected a subset of the bioclimatic variables developed by Hijmans, Cameron, Parra, Jones, and Jarvis (2005), at the finest grid resolution (30s), that likely influence aphid biology: the maximum (BIO 5) and minimum (BIO 6) temperatures of the warmest and coldest months of the year; the mean temperature of the warmest and coldest quarters of the year (BIO 10 and BIO 11 respectively). Temperature variation affects insect development and reproduction (Ratte, 1984), including those of Cinara species (Durak & Borowiak-Sobkowiak, 2013). The aphid life cycle is also likely affected if temperatures exceed threshold values, corresponding to heat stress or extreme cold (Campbell, Frazer, Gilbert, Gutierrez, & Mackauer, 1974). Because drought episodes can also have a significant effect on aphid colony development (Rouault et al., 2006), we also selected the precipitation of the driest month (BIO 17).

| Taxonomic validation and association of occurrences with phylogenetic species
The identification of aphids based on morphological criteria can be ambiguous because aphids show considerable overlap in their morphological characteristics (Coeur d'acier et al., 2014). Specimen misidentification or other taxonomic ambiguity can lead to an erroneous geographic occurrence record of a species. We therefore sought to cure the occurrence dataset of misidentifications: We downloaded 743 Cinara barcode sequences (i.e., cytochrome oxidase I DNA sequences available in IBOL [http://www.bolds ystems. org/] and GenBank [https ://www.ncbi.nlm.nih.gov/genba nk/]) representing 105 species and corresponding to 638 records in our occurrence dataset. Cross-referencing the datasets allowed us to assess whether some of the species included in this study were not monophyletic and could thus be subject to taxonomic ambiguities. To this end, we reconstructed a maximum-likelihood tree using a Kimura-2-parameter model of base substitution (Kimura, 1980) and 1,000 bootstrap replications as implemented in RaxML 8.2.10 (Stamatakis, 2014). The COI tree revealed some polyphyletic and paraphyletic species. In most cases, we defined monophyletic species groups of closely related species that seem to be regularly mixed. We then grouped occurrence data accordingly and assigned them to their respective species group. In cases where sequences assigned to a species were found in distant phylogenetic lineages, the species was either entirely removed from our occurrence dataset and the final phylogenetic tree, or just the evidently problematic occurrences were deleted. Details on the curation of the occurrence dataset and clustering of specimens into species groups are given in Appendices S1 and S2.

| Phylogenetic methods and molecular dating
In order to obtain the most comprehensive Cinara phylogeny, we complemented the nucleotide dataset published by Meseguer et al. (2015) with sequences available on GenBank. We used one marker from the aphid nuclear genome (elongation factor 1-α, EF) and two from its mitochondrion (COI; cytochrome b, CytB). Then, we considered two markers from the aphid's primary symbiont, Buchnera aphidicola: a chaperonin molecule involved in protein folding, GroEL, and "His" including the ATP phosphoribosultransferase, HisG, and the histidinol dehydrogenase, HisD. The final dataset included 82 Cinara species, 15 more than the Meseguer et al. (2015) study.
Phylogenetic analyses were performed with the single concatenated DNA matrix with MrBayes 3.2.3 (Ronquist et al., 2012).
We used BEAST V.1.8.2 (Drummond & Rambaut, 2007) to estimate divergence times in the phylogeny using the topology inferred by MrBayes as a fixed tree. Following Meseguer et al. (2015), we used a Yule tree prior and a relaxed lognormal molecular clock. We first applied the GTR+G+I substitution model to all gene partitions but, detecting over-parametrization for the GroEL and EF genes in Tracer (ESS < 200), we used the HKY+G+I model for these two markers. We ran two MCMC chains of 80 million generations each, sampled every 1,000 generations. We discarded 25% as burn-in and combined the two independent runs with LogCombiner 1.8.3. Three calibration points were used to calculate absolute ages of divergence in Cinara, following Meseguer et al. (2015) (Table 1). Additionally, we performed the dating analysis with a random starting tree in order to obtain a posterior sample of trees while accounting for topological uncertainties.

| Climatic niche evolution in Cinara
For the following analyses, we pruned from the tree the species for which we had gathered fewer than five georeferenced occurrence records. The phylogenetic signal of each selected climate variable was estimated with the λ of Pagel (1999) on the time-calibrated tree. λ varies between 0 (no phylogenetic signal, i.e., the character evolved independently of the phylogeny) and 1 (there is phylogenetic signal, i.e., the phylogeny fully represents the trait covariance among species). From the occurrence data, we calculated the median, minimum, and maximum values of each bioclimatic variable for each species or species group included in the tree and optimized the value of λ using the phylosig function of the phytools package in R (R Core Team, 2013;Revell, 2012). To test whether λ was significantly different from zero, we compared the observed likelihood value with the obtained likelihood value with a fixed λ of zero using a likelihood ratio test.
In order to describe the evolution of climate tolerance throughout the Cinara species' evolutionary history, we used the disparity through time (DTT) method (Harmon, Schulte, Larson, & Losos, 2003) that depicts how trait variation is distributed and when trait values diverged along the phylogeny. This method starts at the root of the tree and calculates the disparity among all existent clades up to the most recent node. The disparity at each node is the mean of the Euclidian distances for n dimensions (n being the number of variables) calculated for each pair of species of the clade. These disparity values are plotted against node ages. To compute disparities values, we extracted the median values for each bioclimatic variable for each species and used the dtt function in the geiger package (Harmon, Weir, Brock, Glor, & Challenger, 2007). Each variable was standardized so that each would have the same weight in the disparity calculation. To consider phylogenetic uncertainty, DTT values were estimated for 5,000 sampled trees from the BEAST posterior distribution inferred with a random starting tree. Low values of disparity suggest that species in the considered clades exhibit similar climatic niches, whereas high values mean that the species are adapted to different climates. Observed disparity values were compared with the values obtained under a pure Brownian process by simulating the evolution of climate variables 1,000 times across our tree using the covariance matrix with the dtt function. Loess regression curves were used to represent the observed and simulated disparities.

| Cinara and Pinus niche equivalency tests
The clades B 2 , B 3 , and B 5 of the Cinara tree (see below) encompass species feeding exclusively on pines (44 species in our dataset). We restricted our analyses of niche equivalency to these species because Pinus-Cinara associations are particularly well documented (Blackman & Eastop, 1994). In addition, occurrences available in public databases often specify the Pinus species whereas the identification of other Cinara hosts frequently does not reach the species level. We retrieved the association of each Cinara species with species of Pinus from Blackman and Eastop (2019) and from our own dataset. When Cinara species were clustered into species groups, the host niche represented the combination of the occurrences of all the hosts of all aphid species included in the group.
In order to test whether the realized climatic niches of Cinara species and their hosts were similar, we performed the niche  Broennimann et al. (2012). For each Cinara species and for each host species range, we considered the environmental space as a rectangle with four sides defined as the minimum and maximum latitude and longitude where the species was found. To stay conservative in our analyses and try to encompass the whole environmental space, we broadened the rectangle by 10° in all directions. We considered the five bioclimatic vari-

| Occurrences and taxonomic validation
After filtering for duplicates and dubious locations, 1,756 of the 6,135 occurrences were retained. The COI tree revealed a few polyphyletic and paraphyletic species, suggesting these nominal species are regularly misidentified or harbor cryptic taxonomic entities (Appendix S1). Indeed, they represent well-known cases of ambiguously defined species and represent taxonomic challenges in need of work (Appendix S2). This kind of phenomenon represents a possible source of error when using georeferenced data from unverified databases (e.g., GBIF) to plot on the tips of a phylogenetic tree. We defined 13 species groups that each encompassed closely related species that were regularly mixed in the COI tree and grouped occurrence data accordingly ( Table 2). Each of these species groups gathered nominal taxa that have been the subject of taxonomic discussion and are often recognized as a single species when using species delimitation methods (Jousselin et al., 2013). In addition, we deleted all 31 occurrences that were assigned to C. pinea in China.
Specimens identified as C. pinea in China were split into two distant genetic lineages (they either clustered with C. pinea from the rest of the world or appeared as the sister species of C. atrotibialis). This

| Phylogeny and dating analyses
We obtained a 2,904 bp alignment (with 22% missing data). The phylogenies inferred with MrBayes and BEAST were well-supported overall and consistent with the one estimated by Meseguer et al. (2015). Two main clades emerged (Appendix S3). Clade A is composed of species feeding on Picea, Abies, and Cupressaceae, whereas

| Cinara niche characterization and evolution
After pruning the tree of species represented by fewer than five occurrences, the dataset included 67 species, with 5-93 occurrences each (mean = 26.9, SD = 20.2). Mapping the bioclimatic scores on the tips of the phylogenetic tree showed that none of the clades is restricted to specific climatic conditions (Figure 1). There is no phylogenetic clustering of climatic tolerance observed among the tips, with each species occupying a different subset of the climatic range of Cinara genus as a whole. The phylogenetic signal was low for four out of the five selected bioclimatic variables (Table 3)

| Niche equivalency tests
The geographic distribution of each aphid species and its range of host plants are presented as maps in Appendix S6. About two thirds (n = 24,  (Figure 3). On the other hand, there are 10 Cinara species that exhibit a narrower climatic tolerance than their host plants (Table 4). This is the case, for example, of C. brevispinosa that is known to feed on P. contorta, P. radiata, and
Although they represent a substantial portion of terrestrial biodiversity, little is known about the phylogenetic patterns of climatic niche evolution in insects over long evolutionary time scales (but see Hidalgo-Galiana, Sánchez-Fernández, Bilton, Cieslak, & Ribera, 2014;Kellermann et al., 2012) and this is especially true for phytophagous insects. The relative lack of phytophagous insect studies is due in part to the limited knowledge of the insects' distribution and host association, but also because their host specialization confounds accurate niche estimation (see below). In order to understand the relative role of climatic niche and the evolution of host preference evolution on phytophagous insect history, a species-level association matrix is needed. These matrices are particularly hard to compile for insects that do not feed on their hosts as adults. Aphids feed on their hosts during their entire life and are thus good models for studying plant-insect associations.
Using an occurrence dataset and a robust phylogeny, we show that climatic niche tolerance is evolutionarily labile in the aphid genus Cinara, with the most of the analyzed temperature and precipitation climate variables displaying a weak phylogenetic signal (Table 3). These results suggest that none of the Cinara clades appeared specialized for any specific climate conditions. On the contrary, all clades encompassed most of the climatic niche variation described within the genus and the observed disparity values for closely related species are greater than simulated under neutral evolution (Figures 1 and 2). These results may indicate a rapid evolution of climate tolerance in the genus and suggest that climatic niche divergence, and adaptation may have played an important role in its diversification and distribution. Alternatively, the climatic niche plasticity observed in our study might simply be the result of the aphids' associations with hosts occurring in disparate climatic conditions. Note: The tests were performed only for pine-feeding Cinara represented by at least five occurrences. Gray highlighted species are the ones for which the test is significant.
Under this latter hypothesis, the realized climatic niche of the specialized insect may actually be narrower than its fundamental niche (i.e., the set of climatically suitable conditions without the constraints of biotic components) and determined by the availability of suitable host plants rather than by the aphid's own physiological limits (Soberon & Peterson, 2005). Although correlative approaches based on distribution data cannot indicate whether host-specialized insects have a wider niche than their hosts, they can reveal situations where insects occupy narrower niches.
Our niche equivalency tests found that a majority of the Cinara species included in this study exhibited a realized climatic niche equivalent to that of their hosts (  (66) P. cembroides ( stem from a generally lower number of occurrences for specialist species (mean for specialist 15.4, mean for generalist 34.2) which might limit the power of the niche equivalency test. Although the niche equivalency test we used is supposed to be robust to small sample sizes, it is mainly designed to limit type I error (i.e., incorrectly rejecting niche equivalency, Broennimann et al., 2012). This may also be explained by the specialists' trend to have a more limited available space. In our study, all the species with climate tolerance narrower than their host's occupied the colder range of their host's climatic niche and were absent from warmer regions. It has long been noted that aphids are negatively impacted by heat and drought stress (Hazell, Pedersen, Worland, Blackburn, & Bale, 2008). Heat stress has been shown to impact aphid's development and mortality (Ma, Hau, & Poehling, 2004;Montllor, Maxmen, & Purcell, 2002;Nguyen, Michaud, & Cloutier, 2009), and seasonal changes in photoperiod and temperature impact the induction of sexual reproduction and the production of frost-resistant eggs (Trionnaire, Hardie, Jaubert-Possamai, Simon, & Tagu, 2008). Studies have also shown that aphid population dynamics have been affected by recent global changes, possibly due to temperature increases (Hullé, Coeur d'Acier, Bankhead-Dronnet, & Harrington, 2010), with global aphid diversity and abundance decreasing with increasing temperatures (Dixon, 1987). In addition, other indirect effects of temperature could explain this inverse latitudinal diversity gradient, for example, a lower efficiency of host plant location in warm environments with greater diversity but lower specific abundance of plants, or the presence of more competitors under warmer environments (Dixon, 1987). Alternatively, the narrow range of climatic tolerance in some Cinara species may be due to historical barriers limiting their expansion throughout the full geographic range of their hosts (Meseguer et al., 2015). shows that a large geographic range has been covered for the entire group. We emphasize also that when Cinara species were sampled (at least by the authors of this study) the aphids were sampled broadly on many conifer hosts, with no focus on any particular species in any particular locations. Therefore, the absence of a particular species in a well-covered area in our dataset, such as Southern California, is probably not a sampling artifact.
The broad sampling reported here should largely obviate possible underestimation of the range of climate conditions under which these species occur.
As a perspective, we must keep in mind that our analysis focused on broad climatic conditions at macroevolutionary scales and did not take into account microhabitat conditions. Future studies could focus on other parameters that are important for aphids. For example, some species migrate from branches to tree roots during summer heat (Durak, 2014;Struble, Osgood, & Pepper, 1976). Hence, different species that occupy the same broad climatic regimes may also experience different microclimatic conditions permitting the coexistence of species on the same host (Favret & Voegtlin, 2004).
The natural enemies and competitors encountered by aphids might also drive biogeographic patterns at regional scales. A combined characterization of aphid tolerance at the macro-and microclimatic scales could provide a more complete characterization of the interactions of the insects, their host plants, and the climate they experience.

| CON CLUS IONS
The realized climatic niche of a phytophagous insect is necessarily associated with its host plant distribution. It is therefore likely that occurrence data for these insects represent incomplete characterizations of their potential ranges and climatic niches that can bias inferences of niche evolution at geological time scales (Saupe et al., 2017). Here, we show that most of the Cinara species investigated occupy a climatic niche equivalent to that of their hosts, suggesting that the main constraint underlying Cinara species' distributions is the physical presence of their hosts. However, we also detected several species that occupy a narrower climatic niche than their hosts, showing a preference for colder conditions. These results suggest that even in host-specialized herbivorous insects, the niche cannot always be fully described in terms of an insect's host association and that adaptation to abiotic factors should be taken into account.
Our results also have bearing on discussing the future of biodiversity in the context of current climate change. We are witness to a dramatic increase in atmospheric CO 2 concentration due to human activities and a concomitant increase in global mean temperature (Waters et al., 2016). The speed at which climate is changing today will substantially alter the global distribution of species over the course of this century (WWF, 2016). Different insect lineage have been already found to track cooler habitats in response to climate warming (Paradis et al., 2008;Parmesan, 2006). The fact that many aphid species cannot track their hosts throughout the warmer areas of their geographic distribution suggests they risk being affected by climate warming. This will most probably induce population isolation, especially in mountain regions where aphids are diverse (Huang, Lei, & Qiao, 2008;Palmer, 1952). The geographic distributions some of these aphid species will probably move northward, increasing the potential threat they can represent for their host trees in these regions. Alexandre Dehne Garcia and the CBGP HPC computational platform. We would like to thank three anonymous reviewers for helpful comments on earlier versions of this manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest. This article does not contain any studies with human participants or animals performed by any of the authors.

AUTH O R CO NTR I B UTI O N S
PA, EJ, ASM, and ACA designed the study. PA and ASM performed the analyses with contributions from MG. ACA, EJ, CF, and Qiao G.X. provided the occurrence data. PA, EJ, and ASM wrote the paper with contributions from ACA, CF, and MG.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sampling locations, Cinara/Pinus association matrix, accession numbers and dated tree are available on Dryad: https ://doi.org/10.5061/ dryad.1r7j2r7.