Future changes in the distribution of two non-indigenous orchids and their acquired enemy in Puerto Rico

Establishment of new populations is contingent on overcoming abiotic and biotic barriers. While this applies to all species, these hurdles are at the forefront of invasion biology where prediction, prevention, eradication, and control strategies depend on an understanding of these processes. Terrestrial Arundina graminifolia and epiphytic Dendrobium crumenatum are two non-indigenous orchids spreading throughout Puerto Rico. The two species have acquired a native herbivore and seed predator, the orchid-specialist weevil, Stethobaris polita. With recently acquired presence records of the three species, land cover data and downscaled climate variables, we modeled their potential distributions under current conditions and also those projected under the least and most extreme climate scenarios for 2050 and 2070. We show that D. crumenatum flourishes in urban environments which also provide refugia from S. polita, whereas there is currently limited refugia for A. graminifolia from S. polita attack, as this orchid has similar climatic niches as the weevil. Projections into all climate scenarios suggest range retractions for all species, with a decreased extent of both orchid populations subject to S. polita attack. Thus, we illustrate for island invasions how climate change will likely alter the distribution of acquired biotic interactions.


Introduction
Over the past century, our natural world has dealt with new, largely human-induced changes at an unprecedented rate and scale. Increasingly obvious are the changes brought about by biological invasions which have become relentlessly pervasive worldwide (Stohlgren et al. 2008;Seebens et al. 2017). Islands are particularly susceptible to invasions (D'Antonio and Dudley 1995;Denslow et al. 2009;Pulwarty et al. 2010), where some impacts vary from positive to negative, and obvious and dramatic to subtle and unexpected (e.g., Vitousek and Walker 1989;O'Dowd et al. 2003;Sin et al. 2008;Abelleira Martínez et al. 2010;Recart et al. 2013). In Puerto Rico, barrages of invasive species are spreading across the island, including greater numbers of non-native orchids (e.g., Rojas-Sandoval and Acevedo-Rodríguez 2015; Falcón and Tremblay 2018;Ackerman 2007Ackerman , 2017. Although orchids are disproportionately underrepresented among invasive species (Daehler 1998), they can have both positive and negative impacts (e.g., apparent competition, Recart et al. 2013). For this reason, it is important to identify recent invasions and factors associated with their establishment, spread and consequences to invaded communities.
Upon introduction, invaders immediately interact with the local biological community, and the interactions they form can have dramatic consequences for invasion dynamics. Whether positive or negative, they can heavily influence the rate of spread and population health of invasive species (de Araújo et al. 2014;Svenning et al. 2014). Upon departure from home ranges, alien species usually leave their enemies behind which may enhance invasion success through processes involving enemy release (Shea and Chesson 2002;Allendorf and Lundquist 2003;Fagan et al. 2005). On the other hand, the acquisition of negative interactions may hamper or even prevent invasion success providing some degree of biotic resistance (Elton 1958;Thompson et al. 2007;Falcón et al. 2017). Yet the lack of shared evolutionary history between invaders and acquired enemies should prevent these interactions from being spatially homogenous across environmental gradients (Fine et al. 2004). In this way, identifying areas of refuge from acquired enemies is an integral part of accurately predicting invasive spread.
Natural habitat heterogeneity is normal, but human activities such as urbanization disturb environments well beyond what native species generally experience, creating opportunities for pre-adapted, non-indigenous species (Daehler 2003). The divergent abilities among native and exotic species to adjust to urbanized worlds can create refuge from biotic resistance. Indeed, urban areas facilitate non-native prominence by weakening any biotic resistance of native species (Holway 1995;Vila and Pujadas 2001;Burton et al. 2005;McKinney 2006;Soifer and Ackerman 2019).
Besides the more obvious local effects of human activities on environmental conditions, regional climate change may also affect dynamics of invasions, especially along the leading and trailing edges of an invasion front ( Van der Putten et al. 2010;Hillerislambers et al. 2013;Lankau et al. 2015). These changes are expected to alter current species distributions, either through latitudinal/altitudinal shifts, range expansions or contractions (Walther et al. 2002) and the responses are likely to vary among species, altering the outcome of species interactions (Menéndez et al. 2008). In this way, climate change might create or eliminate areas of refuge, while also ameliorating or exacerbating rates of spread in new territories.
Herein we focus on populations of the terrestrial Arundina graminifolia (bamboo orchid) and the epiphytic Dendrobium crumenatum (pigeon orchid), two species native to tropical and subtropical Asia and adjacent islands but are now well-established in Puerto Rico. As both invasions are actively spreading across the island from east to west, understanding their fullest possible extent ahead of time can provide useful information for future research and/or mitigation efforts on either species.
While modelling distributions and recording climatic responses, it is important to contextualize both invasions with their known ecological interactions and with the impending changes in the current climate system. Both invasive orchids have acquired a new enemy in the form of the native herbivore and seed predator, Stethobaris polita (Curculionidae), an orchid-specialist weevil known to compromise fruit set and seed production resulting in strong demographic effects on its hosts (Recart et al. 2013;Falcón et al. 2017). The implications of this interaction for both species on the island has not been quantified; however, the introduction of additional exotic orchids may increase risk to native orchids through apparent competition (Recart et al. 2013). In this way, understanding the extent of these interactions with S. polita across the island can also be particularly useful in understanding their greater ecological impacts.
Species distribution models (SDM) are frequently utilized to project species range shifts under climate change, and recently, they have been used to predict distribution differences between interacting species. Most such studies focus on mutualistic relationships (e.g. plant-pollinator) (Kolanowska and Jakubska-Busse, 2020;Malekia and Sadeghi, 2020), while fewer examine how climate change will influence the scale of negative biotic interactions (Menéndez et al. 2008;Aryal et al. 2016).
Herein, we investigate the current spatial extent of both orchid invasions while identifying areas of refuge from weevil attack. Additionally, we explore how the distributions of all three species, along with the extent of their interactions, might shift under future climate change scenarios. We anticipate the lowland D. crumenatum to exhibit an upward altitudinal range shift, increasing the geographical extent of interactions with S. polita, while the highland A. graminifolia will likely exhibit a range contraction under climate change, possibly creating a natural mitigation for the invasion on the island.

Study species
Current knowledge on both orchid species in Puerto Rico, beyond herbarium collections, is rather sparse. The first naturalized population of A. graminifolia was discovered growing terrestrially on a road bank in Río Grande, Puerto Rico, a few kilometers from an orchid nursery that had a cluster of them on their grounds, so it is likely they arrived via the horticultural trade as has happened elsewhere (Ackerman 1995(Ackerman , 2012. The first recorded naturalized population of D. crumenatum was discovered in 2005 growing epiphytically on Annona glabra (Annonaceae) in the Sabana Seca marsh, Municipality of Toa Baja (Ackerman and Collaborators 2014). Informal polls at meetings of three orchid societies in San Juan within the last 10 years revealed D. crumenatum is widely cultivated in home gardens and often sets fruit. Since their original discoveries, both species have spread across the island, as evidenced by local herbarium records (http://herbario.uprrp.edu/, accessed 19 Jun 2019).
Both orchid species exhibit signs of high climate sensitivity, making their responses to climate change particularly important in modelling their distributions. Arundina graminifolia is found naturalized at higher elevations, indicating a sensitivity to higher temperatures. Previous research has established that under any climate change scenario, habitat suitability declines globally for A. graminifolia, although the specific implications for this loss are unknown for Puerto Rico (Kolanowska and Konowalik 2014).
Arundina graminifolia flowers are self-compatible, but seed production is generally pollinator dependent (Huda and Wilcock 2012). A self-pollinating form exists (Forbes 1885), but it has yet to be reported in its invasive range. In Puerto Rico we have seen Centris haemorrhodalis and Apis mellifera (Hymenoptera: Apidae) enter the lightly fragrant flowers barren and leave with a pollinarium attached to the scutellum. In its native range, A. graminifolia is pollinated by bees and wasps from at least five families (Sugiura 2013;van der Pijl & Dodson 1966). While the buds have extrafloral nectaries, open flowers lack any pollinator reward, legitimate visits are rare, and fruit set is low as is typical of orchids employing deceit to attract pollinator (Tremblay et al. 2005).
All orchids must overcome another life history bottleneck: procurement of a mycorrhizal fungus for successful seed germination. The dust-like wind dispersed seeds lack endosperm so it is imperative to gain a fungus to procure nutrients for growth and development (Rasmussen and Rasmussen 2014;Yeh et al. 2019). Meng et al. (2019) discovered a Tulasnella species (Basidiomycota: Tulasnellaceae) in A. graminifolia seedlings collected from a native population in Yunnan, China, and these were capable of inducing seed germination in vitro. Some orchids exploit a broad spectrum of Rhizoctonia fungi, whereas others are more specific (Otero et al. 2002). We do not know the degree of specificity engaged by A. graminifolia, nor how widespread their mycorrhizal fungi are. However, given the ease at which A. graminifolia has become invasive, we would expect that they would be able to exploit a broad spectrum of orchid mycorrhizae, or would specialize on one or few widespread fungi (e.g., Bayman et al. 2016).
The potential distribution of the epiphytic D. crumenatum or its response to climate change has not been modeled. However, a strong-and perhaps unusual-climatic response may be expected for D. crumenatum because epiphytes, which by virtue of being exposed to the elements, are especially susceptible to vagaries of weather (Benzing 1990;Zotz 2016). The gregarious flowering of D. crumenatum is triggered by rain-induced temperature drops (Brooks and Hewitt 1909;Seifriz 1923;Goh et al. 1982). Changes in the frequency of these temperature drops may influence habitat suitability for D. crumenatum. For current distribution models of D. crumenatum to remain applicable in the coming decades, they must account for scenarios with an altered climate. On the other hand, D. crumenatum may show unusual resiliency since it readily establishes on urban street trees in both its native and invasive ranges (Ackerman 2012;Leong and Wee 2013).
Dendrobium crumenatum is a self-incompatible species (Johansen 1990), but the strongly fragrant flowers produce nectar attracting honeybees. In its native range, D. crumenatum is pollinated by Apis cerana and A. dorsata (Brooks and Hewitt 1909;Leong and Wee 2013), but in its invasive range (Hawai'i and Puerto Rico), non-indigenous Apis mellifera is a capable surrogate (Ackerman 2012(Ackerman , 2017. It has also naturalized in Jamaica, Guadeloupe, and the Seychelles (C. Hamilton in Ackerman and Collaborators 2014;Feldmann and Barré 2001;Jolliffe 2010, respectively).
Seed germination of D. crumenatum in Singapore occurs with Basidiomycota genera Ceratobasidium sp. (Ceratobasidiaceae) and Epulorhiza sp. (Tulasnellaceae), the former of which appears to be a locally common fungus (Ma et al. 2003;Izuddin et al. 2019). There are no reports of fungal associates in the invasive range of D. crumenatum.
Stethobaris polita is a small black weevil (2.8-3 mm long) first described from a Vanilla flower in Puerto Rico. It also occurs in Guadeloupe, Saint Vincent and Dominica of the Lesser Antilles (O'Brien and Turnbow 2011; Recart et al. 2013). Adults eat orchid flowers, occasionally apical meristems, and lay eggs in fruits and sometimes in peduncles. They attack both native and non-indigenous orchids. Once relatively rare, the beetles are now common in populations of non-indigenous orchids Arundina graminifolia and Spathoglottis plicata, both of which flower and fruit all year (with seasonal peaks). However, the weevil has thus far failed to inhabit urban landscapes because they apparently lack the ability to disperse and/or establish within such environments, thereby creating refugia for orchids that have the ability to thrive in such environments (Soifer and Ackerman 2019).

Study site
During June and July 2019, we crisscrossed the Caribbean island of Puerto Rico (18.2208°N, 66.5901°W) searching for Arundina graminifolia and Dendrobium crumenatum. Puerto Rico has six Holdridge Life Zones: subtropical dry forest, subtropical moist forest, subtropical wet forest, subtropical rainforest, subtropical lower montane wet forest, and subtropical lower montane rain forest (Ewel and Whitmore 1973). We covered all but the dry forest habitats which are not known to harbor either orchid species or the weevil.

Collection of occurrence records
We obtained coordinates using Garmin eTrexÒ GPS and recorded plant data for populations along trails and roadsides, and on public and private properties. The locations of individual populations were obtained from herbarium collections and personal communications with local field biologists. Some populations were also identified through incidental discovery, either during travels between known populations, or through intentional scouting in previously unexplored municipalities. Ultimately, we collected population data from the following 19 municipalities: Adjuntas, Aibonito, Barranquitas, Caguas, Canóvanas, Carolina, Cayey, Cidra, Jayuya, Lares, Luquillo, Naranjito, Orocovis, Patillas, Río Grande, San Juan, San Sebastian, Trujillo Alto, and Utuado. Unsurprisingly, most populations of both orchid species were observed along the island's extensive road network (https://data. humdata.org/dataset/hotosm_pri_roads; accessed 13 April 2021). Arundina graminifolia grows naturally in disturbed, open habitats (Dressler 1981;Chen et al. 2009) which roadside banks provide, and Dendrobium crumenatum is often cultivated in home gardens and roadside trees.
Populations were only considered if there was evidence of naturalization and/or recruitment as follows: the presence of reproductive structures (flowers or fruits) on at least one plant; the presence of recruits in the population; or sufficient circumstantial evidence that the population was not cultivated, such as the plant being inaccessible to the general public (i.e. on top of a steep roadside hill or high up a tree). Even if a population was originally cultivated, it was still included if successful recruitment was evident, acknowledging the role of the horticultural trade as a vector in non-native orchid dispersal.
To reduce spatial autocorrelation in the model, orchid populations were only sampled if they were over 1 km distant from each other (i.e. the resolution of 30-arcsecond bioclim layers) from previously sampled populations, as calculated using car odometers.
At every population, we recorded the number of S. polita present, along with the presence of other possible enemies. We categorized weevil presence at a waypoint from either weevil sightings or evidence of weevil damaged reproductive structures (buds, flowers, fruits).
Samples included a total of 30 populations of D. crumenatum, 28 populations of A. graminifolia, and 140 populations of S. polita. Since presence points for S. polita would be limited to the presence points of our two orchids, their modelled distribution might create an incomplete picture of their overall distribution. To counteract this issue, we combined our S. polita presence points with those recorded with a more widespread alien orchid, Spathoglottis plicata (Recart et al. 2013;Soifer and Ackerman 2019). For this reason, the number of populations of S. polita is higher than those of either of our focal orchids.

Distribution modelling and mapping
We modelled the distributions of A. graminifolia, D. crumenatum and S. polita using Maxent v 3.4.1 (Phillips et al. 2019), similar to previous studies on invasive species distribution in both the present and future (Recart et al. 2013;Kolanowska and Konowalik 2014;Soifer and Ackerman 2019). Models were considered with all 19 bioclimatic variables from Worldclim Version 2--Global Climate Data (Fick and Hijmans 2017) along with a recent landcover map for the island (PRGAP 2006). Pearson's correlation tests (supplied through SDM Toolbox; Brown et al. 2017) were run for each pair of variables to identify highly correlated variables ([ 0.7), keeping one of them. To determine the appropriate combination of uncorrelated environmental variables for each species, preliminary models were run (replicates of 5-10) with every combination of uncorrelated variables, and the combination with the highest test and training gain was kept for the final model, as identified by jackknife tests on variable importance. As a result, the same five environmental variables were used for modelling each study species: temperature seasonality (Bio 04), maximum temperature of warmest month (Bio 05), temperature annual range (Bio 07), precipitation seasonality (Bio 15), and land cover. Bioclimatic variables had a spatial resolution of 30 arc-seconds (* 1 km) and land cover had a resolution of 15 m 2 (Gould et al. 2008). We prepared the environmental layers for Maxent on ArcGIS with the Extract By Mask tool with the ''Convert Units'' setting to match the extent, cell size, and coordinate system of the PRGAP land cover map. All layers had a final resolution of 15 m 2 and the coordinate system NAD 1983 Stateplane Puerto Rico Virgin Islands FIPS 5200.
Logistic models were run across 100 replicates with bootstrapping. We applied the ''random seed'' option, which randomly partitions data into testing and training data, with 80% of the presence records selected for training and the other 20% for testing. All other settings were kept at their default, unless otherwise noted. For any cell with more than one population, only one was used to contribute to the model. Sampling efforts for both orchids every 1 km resulted in no populations being removed for either A. graminifolia or D. crumenatum. The previous sampling efforts of S. polita with the orchid S. plicata did not follow this same methodology, resulting in only 84 of the 140 S. polita populations contributing to the models.
Using the same variables and settings, along with clamping, we projected the distribution of the three species into 2050 and 2070 using the general circulation models (GCM), CCSM4 (Community Climate System Model 4) and CNRM-CM5 (Centre National de Recherches Météorologiques-Coupled Model 5), across representative concentration pathways (RCP) 2.6 and 8.5 (Fick and Hijmans 2019). GCMs were chosen based on their use in similar studies (Malekian and Sadeghi, 2019; Kolanowska and Jakubska-Busse, 2020; Tsiftsis and Djordjević, 2020). We chose future climate variables from Worldclim Version 1 because this is the only version that is available in 30-arcseconds, the resolution we used for current climate variables and needed considering the size of Puerto Rico (approximately 8000 km 2 ).
Variation among individual GCMs brings into question how accurate any single GCM is for a particular area (Goberville et al. 2015). For example, Khalyani et al. (2016) found that model recognition of Puerto Rico's bimodal precipitation patterns strongly affected climate projections for the island, with bimodal models projecting more severe climate change than unimodal models. While there is no ''best'' model to use for Puerto Rico, previous studies have used multiple GCMs to account for some of the variation in GCM outputs (Goberville et al. 2015;Dyderski et al. 2018). For this reason, we used one GCM (CNRM-CM5) that models bimodal precipitation patterns, and another GCM (CCSM4) that is unimodal for the Caribbean (Hayhoe, 2013). In doing so, we anticipated having one model (CCSM4) provide more conservative projections than the other (CNRM-CM5), accounting for GCM variation. Doing so also tested the consistency of observed patterns against GCM variability.
We applied the threshold rule of equal training sensitivity and specificity similar to previous studies (Soifer and Ackerman 2019) to produce binary presence/absence maps in ArcGIS through the ''Reclassify'' tool in the Spatial Analyst extension, which automatically changes every cell value in a raster of a certain value or range of values to a new, specified value. Threshold values for each species were their average across each run, available through the maxentresults.csv accompanying each model. These threshold values are as specified: Using the Reclassify function in ArcGIS, we created functions where all values below the thresholds were given a new value of 0 for each species (e.g., for A. graminifolia a range of 0-0.2606999 = 0). Then, all values at or above the threshold were reclassified to a unique number for each species (D. crumenatum = 1; S. polita = 2; A. gram = 4). By combining the rasters through the Raster Calculator Tool in ArcGIS, where the resulting raster's cell values are the sums of the binary orchid and weevil rasters, we created a raster with unique values for absence, occurrence of one species, and occurrence of any combination of two species (i.e., cell values of 2 = only S. polita presence, 4 = only A. graminifolia presence, and 6 = both species present at that cell, etc.). This resulted in rasters containing cells of absence, existence and coexistence of each orchid and the weevil in the present day and under each modelled climate scenario.
Attribute tables for these rasters give the number of absence, presence and coexistence cells in each model. Comparisons in their frequencies between climate scenarios shows the corresponding changes in range size and overlap. Specifically, the percent change in the total range size is quantified as the difference in occupied cells of a specific value between a climate scenario and the present day. Additionally, the percent of range overlap is quantified through the following equation: range overlap = (# of occupied cells of both orchid and weevil / total # of occupied cells for the orchid)*100).
Along with modelling species distributions, Maxent provides multiple statistical tests on model significance and performance. We used the area under the receiver operating characteristic (ROC) curve (AUC) to discern model output from a random prediction. We conducted jackknife tests to evaluate variable importance.

Present and potential species distributions and interactions based on current conditions
Field observations confirm that both A. graminifolia and D. crumenatum have acquired a new enemy, S. polita. We found the weevil more often on A. graminifolia than D. crumenatum. Of the 30 populations of D. crumenatum surveyed, only 6 contained evidence of S. polita damage. In comparison, 20 of the 28 A. graminifolia populations surveyed had evidence of S. polita damage to their sepals, lip, and/or petals, including oviposition in A. graminifolia fruits. Although S. polita has been seen feeding on meristematic tissue of D. crumenatum in the absence of flowers (JDA personal observations), we observed weevils during sporadic flowering events, with similar damage to D. crumenatum flowers as observed on A. graminifolia. Ovipositing holes in fruits of D. crumenatum were observed as well.
For A. graminifolia, Maxent predicts suitable habitat within the higher elevations of the Cordillera Central, Sierra de Cayey and Sierra de Luquillo, which matches our field observations of population occurrence (Fig. 1a).
For D. crumenatum, our models show suitable habitat coinciding with urban areas, particularly San Juan and the Caguas-Cayey region, which is where a high density of D. crumenatum populations were observed (Fig. 1b). Dendrobium crumenatum appears to prefer lower elevations and largely disturbed habitat, while being potentially absent from higher elevations of the Sierra de Luquillo (El Yunque National Forest) and the Cordillera Central. The model is largely influenced by land cover, as the training gain is higher with land cover alone than when it is excluded (Fig. S3b). However, areas of high human land use are not the sole variable necessary for survival. Response curves reveal an affinity towards areas with highly seasonal, yet still high temperatures, along with year-round precipitation (Fig. S5). Such conditions are largely found along the lower elevations of the island, except for the dry habitats on the south side of the island, which remain largely uninhabitable for the orchid, even within large urban areas like Ponce (Fig. 1b).
For S. polita, our models show a negative relationship with urban land cover. The intact forests of Sierra de Luquillo (containing El Yunque National Forest) and Sierra de Cayey (containing Carite State Forest) are deemed the most suitable habitat, with hot spots projected in those areas (Fig. 1c). At the same time, urban areas are almost entirely excluded from the extent of S. polita, particularly San Juan and Caguas, which matches field observations. Additionally, S. polita appears to prefer higher elevations in wet conditions, with their range extending along the Cordillera Central, except for portions with the highest elevations, and into the northern karst region (''mogotes'') in the northwest region of the island (Fig. 1c). The scattered patterns across the northern karst region are likely a product of the microclimate conditions from the abrupt elevation changes of the mogotes (limestone haystack hills).
Response curves for both A. graminifolia and S. polita suggest they occupy highly similar climate niches, with analogous responses to each climate variable (Figs. S4 & S6). This results in largely analogous distributions on the island (Fig. 1). However, their distributions are not identical and differences are likely due to small deviations in variable response and importance. In comparison, response curves suggest D. crumenatum occupies different climate niches than either A. graminifolia or S. polita (Fig. S5), which reinforces field observations of limited instances where naturalized populations of D. crumenatum co-occurred with either A. graminifolia or S. polita.
Co-occurrence maps reveal how variations and similarities in response to the selected environmental variables influence the degree of interspecific interactions. The distribution of S. polita excludes areas of high anthropogenic disturbance, providing a refuge for D. crumenatum in urban environments. (Figs. 2f and  3f). Even beyond urban areas, differences in the climate niches of both species likely leads to sparse overlap between the two populations across the island, with S. polita present in less than a quarter of D. crumenatum's range (Table 1). In comparison, the similar response of A. graminifolia and S. polita to environmental variables results in significant overlap in their potential distributions across the island (Figs. 2a and 3a), with roughly three quarters of A. graminifolia's range overlapping with that of S. polita (Table 1). Higher tolerance for lower temperatures might allow A. graminifolia to extend beyond S. polita at higher elevations in the western end of the Cordillera Central. Conversely, tolerance for higher temperatures by S. polita might allow it to extend to lower elevations than A. graminifolia (Figs. 2a and 3a); however, these possibilities should be tested further.
Future species distributions and interactions based on climatic change models As climate changes over the coming decades, so does the potential distribution of these three species. While there is variability among models, all project a decrease in suitable habitat for each species, regardless of GCM, RCP or time period (Table 1). Furthermore, both GCMs agree that severe climate change (rcp8.5) will bring similar to greater range reductions for all three species in the next 50 years than a future of immediate and drastic climate action (rcp2.6) ( Table 1). Both GCMs suggest the ranges of A. graminifolia and S. polita can be half of their current sizes by 2050, even in the most optimistic climate scenarios (Table 1; Fig. 4). In the same climate scenario (rcp2.6), D. crumenatum may see a quarter to a third of its current range no longer suitable. In comparison, under rcp8.5 there can be an even further 10-20% reduction (from present day) by 2050 for all three species (Table 1). By 2070, further reductions in suitable range size might be expected for each species. Under the most extreme climate scenario (rcp8.5), roughly 90% of the current ranges for A. graminifolia and S. polita are projected to no longer be suitable, along with over half of D. crumenatum's current range (Table 1; Fig. 4). There is less consensus in population trends between GCMs and species for rcp2.6 in 2070. For A. graminifolia, the CNRM model suggests certain regions that may become unsuitable in 2050 may become suitable again by 2070, with a larger projected range size in the latter, and with ''re-establishment'' primarily projected in Sierra de Luquillo (Fig. 2b, c); while the CCSM model shows a range size similar to that in 2050 under rcp2.6 (Table 1; Fig. 4). For both D. crumenatum and S. polita, the CNRM model predicts more drastic reduction in range size than the CCSM model under the rcp2.6 climate scenario. For D. crumenatum, the CNRM model suggests its range will be roughly the same size under both climate scenarios in 2070, while the CCSM model suggests, under rcp2.6, it would be reduced slightly more from 2050 to 2070, but the degree of reduction would be roughly half of that under rcp8.5 (Table 1). For S. polita, the CNRM model predicts its range might continue to diminish by 2070 under the rcp2.6 model to levels slightly smaller than under rcp8.5 in 2050. In comparison, the CCSM model for rcp2.6 suggests a smaller reduction in range size by 2070 than is predicted for 2050 under rcp8.5. For both D. crumenatum and S. polita, these differences in GCM predictions bring some uncertainty to how influential current human action on climate change will be on the range sizes for either species in the next 50 years.
Overall, neither GCM produces consistently more conservative projections than the other, often producing highly similar results. In most climate scenarios, the projected changes in range size for each species are less than 12% different (in relation to present day) between the two GCMs ( Table 1). The exception to this is for rcp2.6 in 2070, where there is the largest disagreement between models upwards of 30%.  Indeed, while all models agree that each of the modelled species might have reduced range sizes by rcp2.6 in 2070, the degree at which such reductions occur, and where, are not as clear. Binary maps of these models reveal consistent projections of suitable habitat for the two GCMs. All models project range contractions for all three species, with only small-scale, isolated instances of ranges shifting into previously unsuitable areas (e.g., D. crumenatum establishing in the northwesternmost corner of the island; A. graminifolia ''re-establishing'' in Sierra de Luquillo in 2070 under rcp2.6 (Figs. 2-3); and S. polita expanding slightly outward in the center of the island under some models (Fig. S1)) all of which are miniscule in comparison to projected losses elsewhere. Indeed, A. graminifolia is projected to lose suitable habitat in much of Sierra de Luquillo, except at the highest altitudes under the CNRM model (Fig. 3), almost all of Sierra de Cayey, and large pockets throughout the Cordillera Central (Figs. 2b-e and 3b-e). Similarly, all models project continued D. crumenatum presence in the San Juan-Caguas corridor and along the northern coastline, with varying degrees of range reduction everywhere else (Figs. 2g-j and 3g-j). Finally, all models project S. polita to lose suitable habitat in the Northwest sector of the island, much of Sierra de Cayey, and pockets of the Cordillera Central, with continued presence primarily projected in Sierra de Luquillo (Fig. S1). The main difference between GCMs for all three species is the size of Extracted elevation values at predicted D. crumenatum suitable habitat do not suggest an altitudinal range shift under future climate scenarios, with many climate scenarios predicting analogous to slightly decreased values in average elevation, rather than an increase (Table 2).
Generally, D. crumenatum has lower projected losses in suitable habitat than A. graminifolia and S. polita (Table 1). While projected losses in suitable habitat are largely proportional for A. graminifolia and S. polita, the larger range size of S. polita means its projected losses cover a larger area of the island (Fig. 4). In general, the total spatial coverage of S. polita on the island is projected to decrease more under climate change than either orchid. As the projected range sizes for all three focal species decrease under climate change, so too does the total area where both orchids interact with S. polita (Fig. 4). In comparison to projected reductions in range size, the decreases in overlap between orchid and weevil vary. For D. crumenatum, all models predict the proportion of its range overlapping with S. polita will decrease to some degree under climate change, with as little as a 21% reduction in overlap under CCSM2050 rcp2.6 where only 18.5% of its range potentially overlaps with S. polita, to as much as an 82% reduction in overlap under CCSM2070 rcp8.5, where as little as 4% of D. crumenatum's range is projected to overlap with S. polita (Table 1). Furthermore, all models but one project at least a 33% reduction in the coverage of S. polita in D. crumenatum's range (Table 1). For A. graminifolia, there is less agreement between GCMs in projecting its proportional interaction with S. polita. Projections into 2050 bring marginal decreases in proportional overlap under the CNRM model and marginal increases under the CCSM model (Table 1), suggesting the decreases in A. graminifolia range size likely leads to proportional decreases in its overlap with S. polita. Moreover, under the CNRM model the proportion of A. graminifolia range overlapping with S. polita is reduced by more than 33% by 2070, while the CCSM model projects only 15% and 27% under rcp2.6 and 8.5, respectively (Table 1). In this way, the proportion of A. graminifolia's range overlapping with S. polita may decrease to some degree by 2070 under climate change, but to a lesser degree than D. crumenatum.

Overview
Global invasions are on the rise and are doing so under an increasingly changing climate. To properly deal with these invasions, it is important to understand their potential scale across time and space. Similarly, it is important to identify biotic interactions and environmental variables that influence the establishment and spread of these invasions. If there is biotic resistance to non-native spread, determining the spatial heterogeneity of this interaction can provide insight into the propagule pressure driving the invasion. However, climate change over the coming decades-and beyond-is expected to alter species distributions and demography; potentially altering the spatial scale of invasions and their biotic interactions in the near future. Unfortunately, little to no previous research on how climate change will alter the distribution of invasions and the scale of biotic resistance. By studying the current and potential ranges of two naturalized orchids Dendrobium crumenatum and Arundina graminifolia, along with their acquired florivore and seed predator, Stethobaris polita, we have successfully modelled how climate change may influence the distribution of the two invaders and the spatial extent of their negative interactions with the weevil which may ultimately affect rate of spread or range contraction.

Factors driving present distributions and biotic interactions
Land cover is the factor that is most associated with the distribution of D. crumenatum across Puerto Rico, with an overwhelming prevalence of naturalized D. crumenatum within urban areas. This abundance in urban centers is likely tied to their spread through the formal and informal horticultural trade, similar to many other ornamental plants (Mack and Erneberg 2002;Mack 2003). The more popular a plant is, the greater the propagule pressure and the higher the probability of invasional success (Dehnen-Schutz et al. 2007). The rapid spread of D. crumenatum across the island over the past decade is certainly connected to its increasing popularity among orchidists and other plant enthusiasts. Rapid establishment and spread of D. crumenatum require two critical processes: exploitation of mycorrhizal fungi and pollinator services. The orchid mycorrhizal fungi are parasitized at the seed germination stage as orchid seeds lack endosperm making them dependent on their fungi for nutrition, at least until the plants become fully photosynthetic Bellgard and Williams 2011;Rasmussen and Rasmussen 2014;Meng et al. 2019). We do not know what fungi are involved, but one study in its native range identified Epulorhiza (Tulasnellaceae) as a fungal associate of D. crumenatum (Ma et al. 2003). Regardless, we expect that either the fungus is widespread in the tropics or the orchid has the capacity to exploit a wide range of fungal taxa. On the other hand, pollinator services are clear. Dendrobium crumenatum in its native habitat is pollinated by Apis cerana and A. dorsata (Hymenoptera: Apidae) (Brooks and Hewitt 1909;Leong and Wee 2013), but throughout much of its invaded range (including Puerto Rico), the orchid is pollinated by a highly effective surrogate, Apis mellifera, which has been introduced worldwide. In Puerto Rico, A. mellifera occurs throughout the island, so pollinator services are not likely limiting anywhere on the island (Ackerman 2017(Ackerman , 2021. Since D. crumenatum is a recent invasion, the evidence of recruitment from initial cultivation is profound. Almost all naturalized populations of D. crumenatum appeared to be sourced from nearby cultivated parent plants. Many orchids were planted in front yards along roadsides, and the role of roads as an effective corridor for non-native seed dispersal has already been identified (Parendes and Jones 2000;Mortensen et al. 2009). Rivers are also likely corridors for D. crumenatum spread because windflow vectors the ultralight seeds and the higher humidity promotes seedling establishment of epiphytic orchids (e.g., Scheffknecht et al. 2010;Olaya-Arenas et al. 2011;Crain 2012). Currently, we are aware of only two populations along riparian habitats, but as the orchid continues to spread over the coming decades, we can expect more riverine populations without an identifiable source population, similar to what we currently observe with most A. graminifolia populations.
Land cover provides widespread refugia for D. crumenatum from S. polita attack, adding another complementary aspect in facilitating D. crumenatum establishment within those areas. An increase in abundance of non-native species with a concurrent decrease in native species towards urban centers is a widely observed phenomenon (Almasi 2000;McKinney 2002;Burton et al. 2005). Our findings reaffirm the absence of S. polita in the urban matrix (Soifer and Ackerman 2019). With a majority of D. crumenatum residing almost entirely within urban and suburban areas, there is little co-occurrence with S. polita and limited biotic interaction (Fig. 1).
However, the modelled proliferation of D. crumenatum primarily within urban centers might also be due to spatial sampling bias, with most populations found within the San Juan-Caguas corridor overemphasizing the role that land use has in determining D. crumenatum establishment. Conversely, it might be underpredicting suitability within natural areas where S. polita is abundant. It may be just a matter of time for D. crumenatum to become a common fixture in forested areas. After all, it has only been noticeably naturalized for less than two decades.
The two processes for successful establishment and spread of A. graminifolia are similar as that of D. crumentatum: exploitation of mycorrhizal fungi and pollinator services. The fungus associated with A. graminifolia is a Tulasnella species retrieved from seedlings in a population in Yunnan Province, China (Meng et al. 2019). We do not know whether this is a widespread fungus, or that A. graminifolia exploits other fungi in its extensive native and non-native range. However, pollinator services are clear. Flowers are self-compatible but not self-pollinating (Huda and Wilcock 2012). Unlike the nectar producing D. crumentatum, A. graminifolia attracts pollinators by deceit, offering no reward to foraging pollinators. In its native range, A. graminifolia is pollinated by multiple species, and is a pollinator generalist (Sugiura 2014). In Puerto Rico, A. graminifolia is pollinated by Apis mellifera, but also by Centris haemorrhoidalis, a larger solitary bee. Both occur throughout the current range of A. graminifolia on the island (Ackerman 2017(Ackerman , 2021. While A. graminifolia is a popular ornamental plant in Puerto Rico, it's proliferation along remote roadsides within mountain forests coincides with localities in its native habitat (Dressler 1981;Chen et al. 2009).
Major urban centers that provide orchids refuge from S. polita attack are all found at lower elevations in Puerto Rico, and do not overlap with the distribution of the mid to high-elevation A. graminifolia. As a result, land cover does not provide refuge for A. graminifolia from S. polita attack, except perhaps in the western end of the Cordillera Central (Fig. 2a). This region appears to have climate barriers to S. polita survival (Figs. 2-3), but it accounts for only a quarter of A. graminifolia's current potential range.
The present degree of overlap between both orchids and the weevil suggest interactions with S. polita are not preventing non-native orchid establishment.
Indeed, prolific populations of A. graminifolia exhibited weevil damage. Nonetheless there is likely some degree of biotic resistance as was shown for another invasive orchid in Puerto Rico, Spathoglottis plicata, which had positive growth rates despite heavy infestations of S. polita. Demographically, the weevils extended lag times thereby slowing orchid population growth and just delaying the inevitable (Falcón et al. 2017). Future research is needed to investigate the specific demographic effects of S. polita attack on A. graminifolia and D. crumenatum, and how they might compare to the previously studied effects on invasive and native orchids (Recart et al. 2013).

Distributions and biotic resistance under climate change
While climate change is largely expected to facilitate global invasions, our study shows climate change decreasing the size of suitable habitat for two nonnative orchids in Puerto Rico. Despite each orchid occupying contrasting climate niches, both species are expected to experience similar trends in range contraction due to climate change (Table 1; Fig. 4). Climate responses are regionally-specific, but invasive range-reduction under climate change has been predicted in other studies (Beaumont et al. 2009;Bradley et al. 2009;Kolanowska and Konowalik 2014). While higher CO 2 levels and increased global commerce should continue to promote plant invasions in the future, climate changes may promote some invasions while restricting others (Bradley et al. 2010). Moreover, the degree at which it restricts invasions is species-specific, and this is particularly evident when comparing the projected responses of A. graminifolia and D. crumenatum; where the former is projected to lose between & 50-90% of its current range by 2070, while the latter is projected to lose between & 27-66% of its range in that same time frame (Table 1). Considering D. crumenatum presently has a larger range size than A. graminifolia ( Fig. 1; Fig. 4), it is projected to continue to remain more widely abundant on the island under future climate change than A. graminifolia.
Multiple studies have illustrated how climate change can invoke independent responses among interacting species, creating novel communities and altering existing ones (Schweiger et al. 2008; Van der Putten et al. 2010;Stewart et al. 2015). In the present case, our models generally suggest that by 2070 climate change will decrease the extent by which each non-native orchid will be susceptible to weevil attack. However, these decreases are greater for D. crumenatum than for A. graminifolia (Table 1). This might be due to A. graminifolia and S. polita already occupying similar climate niches, and subsequently interacting at a far greater degree than D. crumenatum and S. polita, while also experiencing analogous levels of projected range contraction under climate change (Table 1, Fig. 4). Moreover, D. crumenatum is projected to principally lose suitable habitat in regions where it presently interacts with S. polita (Fig. 2-3) while simultaneously projected to lose notably less habitat than S. polita in general (Table 1; Fig. 4).
The future distribution of the three focal species can have important implications for the surrounding community. The invasion of non-native orchids in Puerto Rico has been a demographic boom for S. polita by providing abundant resources throughout the year. Unfortunately, this can have spillover effects on native orchids increasing the frequency of attack and reducing their reproductive output (Recart et al. 2013). However, projected reductions in range size for S. polita upwards of 90% by 2070 questions the overall ecological significance of S. polita on both native and non-native orchid populations in the future. While empirical tests would be necessary to quantify the ecological ramifications of reduced S. polita coverage, such drastic reductions might bolster the population dynamics of orchid species, both native and nonnative., Reduced interactions with S. polita for both non-native orchids would likely increase successful seed production and thus propagule pressure while simultaneously reducing negative effects on native orchids in those same areas. In general, climate change can have cascading effects on community structure and most species interactions (Tylianakis et al. 2008;Gilman et al. 2010), which might further influence species establishment within climatically suitable areas.
When interpreting climate change SDMs, it is important to recognize that models work to show what is possible. The degree of projected range shifts can be further influenced by other environmental factors. For example, dispersal capabilities and conditions for germination and establishment will determine if and where species move within climatically suitable areas. In some cases, dispersal limitations can impede invasions more than biotic resistance (Moser et al. 2011). Schweiger et al. (2008) found that the dispersal capabilities of interacting species are highly influential in determining their future distributions, with significant expansions in distribution with unlimited dispersal (goes everywhere that is predicted suitable), but significant decreases when dispersal is highly restricted (future range only includes overlap between present suitability and future suitability). While models only predict small-scale, and sometimes isolated, cases of range expansion or shift for both orchids, which do little to counteract the range contractions modelled elsewhere, dispersal capabilities can remain broadly relevant considering the highly fragmented distributions of both orchid species in certain climate scenarios, particularly A. graminifolia under rcp8.5 (Figs. 2, 3). While true unlimited dispersal is ecologically unrealistic, the continued stratified spread of A. graminifolia and D. crumenatum in the horticultural trade and their production of many thousands of winddispersed seeds within a single fruit should not make either orchid dispersal-limited. All orchids require mycorrhizal fungi for successful germination, but the apparent ease at which the two orchids spread suggests that their distributions may not be limited by the lack of appropriate fungi, at least not beyond the microsite scale. As for S. polita, it does oviposit within orchid fruits, tying its own dispersal to that of its host plants, but the weevil is not dependent on just one host species. Nevertheless, we know little of the current and future fate of the orchid flora. We do anticipate that the weevil likely has niche flexibility since the orchid family is species-rich with representatives in all terrestrial habitats. Stethobaris polita exploits species in all three subfamilies present within its range, and its native orchid hosts often occur in small, hyperdispersed populations (Stewart et al. 2015;Brewster and Ackerman 2013;Tremblay and McCarthy 2014). Since suitable range sizes are projected to decrease for all three species under climate change, any dispersal limitations that prevent establishment within modelled regions of climate suitability would make the range reductions of the three species greater than what our model predicts.
One important consideration to make in interpreting our results is the ever-changing nature of land use. Since we cannot accurately predict landscape or land use changes in the coming decades, land cover was held constant across all climate projections. However, it is expected for land use to change in the future through either human actions (Martinuzzi et al. 2007), climate change (Khalyani et al. 2016), or natural succession (Chinea, 2002). As urban areas move and change, so might the scale of the interaction between S. polita and D. crumenatum.
Contrary to our predictions, our models showed no evidence of an altitudinal range shift for D. crumenatum under projected climate change. In fact, the distribution of projected D. crumenatum populations along the elevation gradient remains consistent from the present day through various projected climate scenarios. With land use being held constant through future climate scenarios, their continued presence in those areas despite an altered climate might be due to tolerance for heat as long as there is sufficient precipitation, or a spatial sampling bias towards urban centers like San Juan, which are at sea level. Even if our sampling was biased, we have no evidence of an altitudinal range shift towards mid to higher elevations (Table 2).
Species distribution models can provide important insight into the future distributions of invasions and of their acquired interactions; however, they provide an incomplete picture of novel community structure and population dynamics within suitable regions. The next step in understanding the future of A. graminifolia, D. crumenatum and other orchid invasions on Puerto Rico is assessing how climate change will influence their interactions with pollinators, ants, herbivores, seed predators, and mycorrhizal fungi; especially since there is growing literature on climate-induced phenological uncoupling between interacting species (Thackeray et al. 2010;Walther 2010). Finally, future research should monitor the effects that A. graminifolia and D. crumenatum have on native orchids through apparent competition (Recart et al. 2013).
Considering the positive attitudes towards invasive orchids, their mitigation and eradication is likely unpopular and possibly not economically worthwhile (Ackerman 2007). Although the overall ecological effects of climate change are expected to be overwhelmingly negative, it can directly and indirectly mitigate invasions where human intervention will likely not occur. While climate change will not eliminate either orchid from Puerto Rico in the next 50 years, our analyses indicate that it might do so if climate trends continue in the next century.
Climate change can have a profound effect on species distributions and interactions, and these effects have important implications within the context of invasions. Non-analogous responses to climate change between invaders and their acquired enemies can create or alter spatially-heterogenous outcomes to such interactions. To the best of our knowledge, this study is one of the first to model the changing spatial relationships between plant invaders and an acquired enemy that may offer some measure of biotic resistance now and in the future.