\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,english]{babel}



\usepackage{float}








\begin{document}

\title{Climate change has different predicted effects on the range shifts of
two hybridising ambush bug (\emph{Phymata}) species}



\author[1]{Vicki Mengyuan Zhang}%
\author[2]{David Punzalan}%
\author[2]{Locke Rowe}%
\affil[1]{University of Toronto - Mississauga}%
\affil[2]{University of Toronto}%


\vspace{-1em}



  
  \date{June 18, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
A universal attribute of species is that their distributions are limited
by numerous factors that may be difficult to quantify. Furthermore,
climate change-induced range shifts have been reported in many taxa, and
understanding the implications of these shifts remains a priority and a
challenge. One approach is to employ species distribution models which
correlates species presence data with a set of predictor variables.
Here, we use MAXENT to predict current suitable habitat and to project
future distributions of two closely related \emph{Phymata} species in
light of anthropogenic climate change. Using species occurrence data
from museum databases and environmental data from WorldClim, we
identified environmental variables maintaining the distribution of
\emph{Phymata americana} and \emph{Phymata pennsylvanica}, and created
binary suitability maps of current distributions for both species on
ArcMap. We then predicted future distributions using the same
environmental variables under different Representative Concentration
Pathways (RCP), created binary suitability maps for future
distributions, and calculated the degree of overlap between the two
species. We found that the strongest predictor to \emph{P. americana}
ranges was precipitation seasonality, while precipitation of the driest
quarter and mean temperature of the coldest quarter were strong
predictors of \emph{P. pennsylvanica} ranges. Future ranges for \emph{P.
americana} are predicted to increase northwestward and southward at
higher CO\textsubscript{2} concentrations. Suitable ranges for \emph{P.
pennsylvanica} are initially predicted to increase, but eventually
decrease with slight fluctuations around range edges. There is an
increase in overlapping ranges in all future predictions. These
differences in optima provide evidence for different environmental
requirements for \emph{P. americana} and \emph{P. pennsylvanica},
accounting for their distinct ranges. Because these species are
ecologically similar and can hybridize, climate change has potentially
important eco-evolutionary ramifications. Overall our results are
consistent with effects of climate change that is highly variable across
species, geographic regions and over time.%
\end{abstract}%



\sloppy


\par\null

\subsection*{INTRODUCTION}

{\label{introduction}}

A universal attribute of all species is that their geographic
distributions are limited. The abiotic and biotic factors that jointly
determine this distribution are expected to be numerous, posing a
serious empirical challenge to their identification and quantification.
One approach is to employ species distribution models (SDM), which
attempt to explain species presence data with a large set of predictor
variables (Elith \& Leathwick, 2009). Although the approach is generally
limited to the use of environmental variables (or their proxies) as
predictors of suitable habitat, SDM have provided new insights into
species requirements that are akin to the ``fundamental niche''
(Hutchinson, 1957) by incorporating constraints set by biotic
interactions (Leach, Montgomery, \& Reid, 2016).

Geographic distributions are, of course, dynamic and dependent upon
changing environmental conditions. Species also vary in their
sensitivity to shifting environmental conditions, and will respond
differently to the same changes (Hickling, Roy, Hill, Fox, \& Thomas,
2006; Malcolm, Markham, Neilson, \& Garaci, 2002), including the
possibility of failure to track new conditions altogether (Loarie et
al., 2009). While there is growing evidence of climate change-induced
range shifts in many taxa, predicting its ecological and evolutionary
implications remains a central challenge (Parmesan, 2006) For example,
climatic variation is undoubtedly linked to natural changes in community
composition over geological timescales, and there is growing evidence of
rapid changes in climate being linked to the invasion and expansion of
alien species (Guo, Sax, Qian, \& Early, 2012). Changing climatic
conditions has also led to more frequent contact between historically
separate species, and this can result in hybridisation (Vallejo-Mar\selectlanguage{ngerman}ín \&
Hiscock, 2016) and, in some cases, species collapse (Njiru, Mkumbo, \&
van der Knaap, 2010). For these reasons, the responses of hybridizing
species to environmental change has been touted as a particularly
important ``window'' on climate change (Taylor, Larson, \& Harrison,
2015).

In the present paper, we took the opportunity to evaluate potential
range shifts in a parapatric pair of insect species that appear to
hybridize in overlapping regions of their respective
ranges.\emph{Phymata} ~\emph{americana} Melin and \emph{P.
pennsylvanica}Handlirsch are two of the most common North American
species in the genus (Family: Reduviidae), with the former more
northerly in distribution, extending west across the American Midwest
and Canadian prairies, and the latter mostly concentrated in the
northeastern United States. Hybridisation in wild populations has been
suspected or inferred (Punzalan \& Rowe, 2017; Swanson, 2013);
consistent with this, current molecular phylogenetic data fail to
distinguish between the two (Masonick, Michael, Frankenberg, Rabitsch,
\& Weirauch, 2017; Masonick \& Weirauch, 2020), despite substantial
morphological divergence (Punzalan \& Rowe, 2017). Both species are
generalist predators occurring in temperate habitats, where they utilize
a wide range of plant species as hunting sites (Balduf, 1939, 1941;
Yong, 2005), suggesting considerable niche overlap. In at least one of
the species, climatic variables (e.g., environmental temperature) play a
key role in thermoregulation, which is linked to mating activity
(Punzalan, Rodd, \& Rowe, 2008a).Thermoregulatory abilities have been
linked to melanic traits, and within- and between-species latitudinal
variation in these traits (Punzalan \& Rowe, 2015) indirectly supports
the importance of climatic variables in the ecology of \emph{Phymata} .
Although there is evidence that their ecological requirements are
consequential to their life histories, there is a shortage of knowledge
regarding their ecology. Thus, there is value in understanding their
habitat requirements and predicting their current and future ranges. We
used biogeographic climate data and SDM to characterize the current and
recent historical range of these two species, and forecast future
distributions under several scenarios of anthropogenic climate change.
We hypothesized that different sets of candidate abiotic factors are
limiting the respective ranges of the two species, resulting in their
current distributions; given anthropogenic climate warming, we
hypothesized that their future overlapping ranges would increase.

\subsection*{MATERIALS AND METHODS}

{\label{materials-and-methods}}

Ambush bug distribution data was compiled from specimens examined by one
of the authors (DP) at the American Museum of Natural History (AMNH),
Carnegie Museums of Pittsburgh, Canadian National Collection of Insects,
Arachnids, and Nematodes, Royal Ontario Museum, Smithsonian National
Museum, University of Guelph Insect Collection, and the University of
Michigan Museum of Zoology. Identifications considered questionable by
DP were excluded from subsequent analyses. These data were supplemented
with information available from museum databases provided by the Spencer
Entomological Collection, at the Beaty Biodiversity Museum
(https://www.zoology.ubc.ca/entomology/) and the Plant Bug Inventory
maintained by the AMNH (http://research.amnh.org/pbi/). We also gathered
data from two citizen science websites, iNaturalist.org and
BugGuide.net. Data collected from iNaturalist.org was included if the
species identity was verified by Paul Masonick (UC Riverside), an
authority on \emph{Phymata} systematics and curator of the iNaturalist
project ``Uncovering the ambush bugs''
(https://www.inaturalist.org/projects/uncovering-the-ambush-bugs). Data
collected from BugGuide.net requires secondary identification
verification before publishing on the \emph{Phymata} webpage and was
assumed to be accurate.

For accessions lacking latitude and longitude data, we supplemented
these data manually using Google Earth. Given that bug sightings
occurred at a specific point, but its accuracy was not reflected on
Google Earth, coordinates were rounded to include degree and hour only
(i.e., minutes and seconds were not used). This is because locality data
were only accurate down to the city level, and, in some cases, were
accurate down to the location of the field station, research station, or
building. For the present purposes, degrees and hours should be
sufficient. The localities of all sightings are mapped in Figure 1.

Although museum and citizen science databases go to lengths to confirm
species identifications, the possibility of misidentifying individuals
is inevitable, as is heterogeneity in sampling methods and reporting. To
mitigate the effects of such errors, the datasets were inspected, and we
removed any data points that were outside of the range of North America
that were a result of errors in the data entry of latitudes or
longitudes, as well as any species identifications that were unreliable
(i.e. indicated as uncertain by the collector/photographer).

Locality data was also divided into two subsets, according to:
``historical'', referring to sightings before 1970, and ``current''
referring to sightings from 1970-present. Preliminary inspection of the
``current'' and ``historical'' distributions of each species suggests
that median latitudes of \emph{P. americana} were decreasing, indicating
a southward shift (Mann-Whitney \emph{U} = 109360 , n = 1075 ,\emph{p}
-value = 0.00019), while \emph{P. pennsylvanica} median latitudes were
increasing, indicating a northward shift (\emph{U} = 127290, n = 970,
\emph{p} -value = 0.00076). This was consistent with our suspicion that
there are opposing changes in the ranges of both species, and motivated
the current study. Additionally, visual comparisons of the historical
and current distributions suggested an increase in the overlap of
\emph{P. americana} and \emph{P. pennsylvanica} distributions (Figure
S1.1 in Appendix 1). For subsequent analyses, we restricted the data
from 1970-2000 (see Table S1.1, Appendix 1 in Supporting Information),
and we refer to this as the current distribution. This was done in order
to prevent temporal mismatch with the bioclimatic data which contained
environmental data from 1970-2000.

We obtained climate data in raster form from the WorldClim database
(Fick \& Hijmans, 2017), which covered all global land areas except for
Antarctica. The grid data was in 2.5 arcminutes (approximately 4.5
square kilometers). The updated (2.0) version of Worldclim's current
environmental data was used to investigate the influence of
environmental variables on ambush bug distributions. This environmental
data ranged from 1970 to 2000 and included 19 bioclimatic variables
(Table S1.2, Appendix 1) derived from monthly temperature and
precipitation measurements. For any pair of variables with Pearson's r
\textgreater{} 0.8, one variable was removed to minimize problems due to
collinearity, based on numerous MAXENT runs (see below). This resulted
in 11 bioclimatic variables used for each species. To map future
distributions, data from future climate variables were collected from
Worldclim's original (1.4) version (as there is currently no updated 2.0
version), using the same 11 bioclimatic variables. Future data were
downscaled from data created using the Community Climate System Model
(CCSM4) global climate model (GCM) from the Fifth Assessment Report by
the Intergovernmental Panel on Climate Change (IPCC). Climate variables
were projected into 2050 and 2070, with four Representative
Concentration Pathways (RCP) trajectories representing four possible
climate change scenarios dependent on atmospheric greenhouse gas
concentration. The best-case scenario for atmospheric greenhouse gases
is represented by the lowest RCP of 2.6. As RCP increases, the
atmospheric greenhouse gas concentration increases as well, up to
RCP8.5, the business-as-usual scenario. Each bioclimatic variable was
scaled to each of the four RCPs; the distributions from these models are
mapped.

MAXENT

We used maximum entropy modeling of species' geographic distributions
(MAXENT), an approach often favoured when restricted to presence data
and considered robust even when data is limited. In the absence of
information about environmental conditions, we assumed the probability
of a species' occurrence within a grid was 0.5 (the default). When a
species was found within a grid for which there is information about
environmental conditions, MAXENT improves the model using the
environmental variables.

We randomly set aside 25\% of the data for the `training' phase, and
assumed a logistic model in all runs, which gives an estimate of the
probability of ambush bug presence within a grid. Our subsequent
analyses were based on the averaged values for 10 runs. A preliminary
run was conducted using all 19 bioclimatic variables to identify
variables that contributed 0\% gain to the model. If these variables
were correlated with other variables that contributed to the model, the
variables with 0\% gain were removed.

We inspected response curves and jackknife plots to evaluate the effect
of different continuous environmental variables, and to determine the
relative importance of these variables. Response curves were generated
for each individual predictor while all other predictors were set at
their average. Geographical projections of the models in the form of
heat maps were used to visualize the predicted probability of ambush bug
occurrence under current and future environmental conditions. Heat maps
were then converted to binary presence/absence maps using thresholds on
ArcMap (ESRI 10.6.1). The threshold used to create binary maps was
``10\textsuperscript{th} percentile training presence logistic
threshold''. This threshold was selected as it assumes that 90\% of the
predicted occurrences will accurately predict the potential range, while
10\% of the predicted occurrences may be erroneous. This results in a
more conservative threshold and is more commonly used with species
distribution data collected over a longer period of time by different
observers (Rebelo \& Jones, 2010).

Raster math calculations were drawn from methods used to calculate the
``suitability status change index'' (SSCI), adopted from (Ceccarelli \&
Rabinovich, 2015). In order to compare the change in suitable habitat,
the future predicted distribution was subtracted from the current
predicted distribution. For both \emph{P. americana} and \emph{P.
pennsylvanica} , current suitable habitat was assigned ``1'' and current
unsuitable habitat was assigned ``0'', while future suitable habitat was
assigned ``2'' and future unsuitable habitat was assigned ``0''. The
difference between current and future predicted distributions resulted
in: ``-1'' = suitable habitat will become unsuitable; ``0'' = unsuitable
habitat remains unsuitable; ``1'' = suitable habitat remains suitable'';
``2'' = unsuitable habitat becomes suitable.

We calculated the percent of overlap between \emph{P. americana}
and\emph{P. pennsylvanica,} projected for 2070\emph{,} to assess changes
in potential contact zones. The predicted ranges of \emph{P. americana}
was subtracted from the predicted ranges of \emph{P. pennsylvanica} and
the same method was used to calculate the change in suitable habitats.
The suitable habitat of \emph{P. pennsylvanica} was reclassified to be
``0'' for unsuitable habitat, and ``2'' for suitable habitat, while the
current suitable habitat and current unsuitable habitat remained ``1''
and ``0'' respectively for \emph{P. americana} . Subtracting rasters
resulted in: ``-1'' = suitable habitat for \emph{P. americana} only;
``0'' = unsuitable habitat for both species; ``1'' = suitable habitat
for both species, indicating potential overlap; ``2'' = suitable habitat
for \emph{P. pennsylvanica} only. The attributes table for each
generated map gives the total number of grids for suitable and
unsuitable habitat. Using these ratios, the percent change in future
predicted distributions under different RCP trajectories and the degree
of overlap between the two species was quantified.

\subsection*{RESULTS}

{\label{results}}

A total of 226 observations were available within the interval for which
climate data was available, 104 for \emph{P. americana} and 122
for\emph{P. pennsylvanica} (Table S1.1, Appendix 1). Using these data,
the models produced by the MAXENT approach were statistically well
supported, as the ratio of true positives (i.e., sensitivity) to false
positives (i.e., 1--specificity) was maximized. Inspection of the
Receiver Operating Curves (ROC), and the Area Under the Curve (AUC) for
both training and test data were greater than 0.90, indicating excellent
model performance (Araújo, Pearson, Thuiller, \& Erhard, 2005).

For both species, precipitation and temperature were identified as the
strongest predictor of occurrence (see Appendix 2). Precipitation
Seasonality (BIO15), the deviation of monthly precipitation from the
annual average, was the environmental variable with the largest relative
percent contribution to the \emph{P. americana} ranges. Response curves
indicate that the highest probability of \emph{P. americana} occurrence
was at a lower precipitation seasonality, or localities with low
variation in monthly precipitation (Figure 2). Isothermality, defined as
a ratio of the diurnal temperature range to the annual temperature
range, was the variable with the greatest permutation importance.
For\emph{P. pennsylvanica} , precipitation was also implicated as an
important factor as indicated in the response curves (Figure 2), though
in this case it was Precipitation of the Driest Quarter (BIO17). The
models suggest an optimal precipitation of about 200 millimeters during
the driest quarter, in which \emph{P. pennsylvanica} has the highest
probability of occurrence. Below precipitation levels of 140 millimeters
and above precipitation levels of 300 millimeters, the probability
of\emph{P. pennsylvanica} occurrence decreases drastically to less than
half of the probability of occurrences at optimal precipitation.
Additionally, there was support for the Mean Temperature of the Coldest
Quarter (BIO11) as it was the variable with the greatest permutation
importance. The probability of \emph{P. pennsylvanica} occurrences was
greatest at a temperature of about -2°C during the coldest quarter.
Above a temperature of 3°C and below a temperature of -6°C, the
probability of \emph{P. pennsylvanica} occurrences decrease to less than
half of the probability of occurrences at optimal temperatures.

Predicted future distributions were mapped as habitat suitability models
for \emph{P. americana} and \emph{P. pennsylvanica} at RCP8.5 (Figure
3), and the predicted changes in suitable habitats are summarized as
percentages (Table 1). Generally, the percentage of suitable habitats is
predicted to increase for \emph{P. americana} . The greatest percent
increase of suitable habitats occurs at RCP6.0, with a 2.7\% increase in
2050 and a 2.0\% increase in 2070, compared to current suitable
habitats. The direction of the range increase is mostly northwestward
and southward. The increase in the percentage of suitable habitats
for\emph{P. pennsylvanica} occurs at a lesser extent when projected for
2050 and 2070. Similar to \emph{P. americana} , the greatest percent
increase of suitable habitats for \emph{P. pennsylvanica} is a 0.8\%
increase which occurs at RCP6.0 in 2050, and the direction of the
predicted range expansion is southward. Notably, there is no change in
the percentage of suitable habitats at an RCP of 4.5 predicted in 2050.
However, predictions for 2070 indicate that the percentage of suitable
habitat either will decrease at an RCP2.6 and RCP4.5 or remain constant
at RCP6.0 and RCP8.5 relative to current suitable habitats. Although
there is no predicted change in the percentage of suitable habitats,
there are some fluctuations around range edges where suitable habitat is
expected to become unsuitable and vice versa.

At all RCP trajectories, there is a slight increase in overlapping
ranges of the two species (Figure \emph{4} ). The largest increase in
overlap is 0.7\%, which occurs at RCP6.0 and RCP8.5. This translates to
a contraction of regions that contain only a single species. For\emph{P.
pennsylvanica} the largest reduction in predicted range area of 2.8\%
occurs at RCP6.0. However, at an RCP6.0, habitats that currently only
contain \emph{P. americana} is predicted to increase. This suggests that
the changes in the amount of suitable habitat will result in the range
shift of \emph{P. americana} towards the range of \emph{P.
pennsylvanica} , resulting in a larger degree of overlap. Additionally,
at RCP2.6, 4.5 and 8.5, the amount of unsuitable habitat decreases, with
the greatest decrease of 0.4\% at RCP8.5 This indicates that a small
amount of habitat that was previously unsuitable is now an area where
either one or both species can inhabit. Conversely, at RCP6.0, there is
a 0.1\% increase of habitat unsuitable for both species.

\subsection*{DISCUSSION}

{\label{discussion}}

Distributions in a changing world

Our models predict different responses of \emph{P. americana}
and\emph{P. pennsylvanica} to anthropogenic climate change, which may
correspond to their respective niche requirements. Our forecasts predict
range expansions of both species is forecasted into 2050, but with a
larger range expansion for \emph{P. americana} (see Appendix 3). This
suggests that the ranges of both species may be able to keep up with
short-term predicted climate change. But by 2070, the range expansion
of\emph{P. americana} is predicted to be more modest, while \emph{P.
pennsylvanica} ranges were predicted to remain the same or contract. It
should be noted that the occurrence and environmental data used in this
study spans 1970 to 2000, which we referred to as the ``current'' range.
Naturally, validating the predictions derived from the models will
require observation and updating as the movement (presumably) proceeds.

A sometimes underappreciated underlying assumptions of the models is
that species are currently in equilibrium with the environment, and that
species ranges are expected to shift as a consequence of changing
environmental conditions (Elith, Kearney, \& Phillips, 2010; Guisan \&
Thuiller, 2005). For example, despite our data suggesting that \emph{P.
pennsylvanica} has a relatively restricted range, it is possible that
both \emph{Phymata} species are still responding to a recent climatic
event, or have already begun responding to recent climate change at the
time that occurrence and climatic data was collected. The predicted
distributions for \emph{P. pennsylvanica} indicates particularly
prominent variation around this species' range edges, possibly
indicating that \emph{P. pennsylvanica} is near its climatic optimum
(Araujo \& Pearson, 2005; Hutchinson, 1957). In comparison, there
appears to be more habitat that satisfies the niche of \emph{P.
americana} given different climate change scenarios. Although it seems
contradictory, at first, that the range of suitable habitat for \emph{P.
pennsylvanica} is not predicted to shift northward in response to
forecasted environmental conditions, we propose several explanations. As
suitable habitats for \emph{P. pennsylvanica} are predicted to be
strongly dependent on winter temperature and precipitation, it is
possible that future winter conditions do not change as drastically as
other environmental variables that have lesser effects, as predicted by
the model. Although it remains to be seen whether the forecasted changes
in climatic conditions are realized, the predicted range expansions (and
potential range overlap) suggest increased hybridisation opportunities
and a larger arena for competition. The potential consequences are
difficult to predict, and depend on a number of factors including rates
of dispersal, the fitness of hybrids, and the possibility of character
displacement (Goldberg \& Lande, 2007; Pfennig, Kelly, \& Pierce, 2016).
Such biotic interactions are not accounted for in our models, but will
almost certainly have important influence on realized distributions
(Hof, Jansson, \& Nilsson, 2012). This also highlights a critical
limitation of SDM in that they typically omit biotic interactions
(Bulgarella, Trewick, Minards, Jacobson, \& Morgan-Richards, 2014), and
the integration of biological interactions with abiotic information
remains one of the frontiers in modeling species distributions
(Anderson, 2017; Elith \& Leathwick, 2009). Furthermore, the models only
make projections of potentially suitable habitats, but do not exclude
the possibility that some populations may successfully persist at or
beyond the predicted range margins of the `preferred' habitat (e.g. due
to local adaptation and/or metapopulation dynamics). There is also no
guarantee that populations will always successfully track spatial shifts
in environmental regimes, in which case the models may underestimate the
possibility and rate of local extirpation. Nevertheless, our models
provide a starting point for generating hypotheses, and adds to a
growing recognition that the current trajectory of climate warming can
have important eco-evolutionary ramifications.

Overall, our results are consistent with effects of climate change that
is highly variable across species, geographic regions and over time
(Menzel et al., 2002). In other taxa, a diverse spectrum of range shifts
have been well documented (Chen, Hill, Ohlemüller, Roy, \& Thomas,
2011). Variability in responses to different climate change scenarios at
different timepoints in the future is seen in studies that have
investigated both individual species (Dowling, 2015; Ning, Wei, \& Feng,
2017) and groups of species (Rebelo, Tarroso, \& Jones, 2010; Urbani,
D'Alessandro, \& Biondi, 2017). Different emissions scenarios (i.e.,
different RCPs) may have opposite effects on distributions, where a
lower RCP induces range expansions and higher RCP projections lead to
range contractions (Wang et al., 2018). Temporal variation has also been
reported, where species were predicted to face extinction due to climate
change at the end of the century, even though current distributions were
predicted to expand (Rebelo et al., 2010). Additionally, predicted
trends of range shifts may also be dependent on the amount of
uncertainty incorporated in climate data sets (Parra \& Monahan, 2008).
For instance, our present study used four climate change projections in
order to capture several potential future distributions, but there are
several other projected concentration pathways that encompass a wider
range of possible future greenhouse gas emissions. Due to the
variability present in these predictions, modelled scenarios should be
used as guides that are ultimately supplemented by additional sampling
or modelling; any long-term trends may be obscured by short-term range
expansions or contractions. The use of SDM such as MAXENT are critical
tools for predicting range shifts but these distributions are contingent
upon the emission scenarios used.

Abiotic determinants of a niche

Species that inhabit the same geographic range may exhibit high
ecological similarity, but imperfect niche overlap will permit
coexistence (Darwell \& Althoff, 2017). The distinct yet overlapping
distributions of \emph{P. americana} and \emph{P. pennsylvanica}suggests
that different bioclimatic variables act to limit ranges. Here, we
identify the variables that are candidates for determining the ranges of
\emph{P. americana} and \emph{P. pennsylvanica} .

Our analyses consistently implicate precipitation as an important
determinant of the abiotic limits of both species, whereby \emph{P.
americana} and \emph{P. pennsylvanica} have different optima. Our
results highlight the possibility that natural selection mediated by
abiotic factors may be specific to life stage (Arnold \& Wade, 1984).
For example, the driest quarter corresponds to the period when eggs of
both species are dormant in winter diapause. During this period, eggs,
which are laid on plant material at the end of summer, are likely to be
near the ground, and the amount of precipitation might translate to
potentially vulnerability to flooding or inundation in the following
spring. Previous work in one species (Mason, 1976) has also invoked
winter conditions during the egg stage as an important period for
triggering phenotypically plastic responses later in life. For many
organisms, strong fluctuations in the environment can be the source of
severe selection, possibly explaining why our analyses also recovered
variability in precipitation (during the summer, when juvenile and adult
bugs are present) as important in predicting occurrences (i.e.,
exclusion).

Mirroring the results for precipitation, our analyses also identified
mean and variability of temperature as potentially important
determinants of geographic distribution. The distribution of \emph{P.
pennsylvanica} was particularly dependent upon the mean temperature of
the coldest quarter, possibly pointing to challenges in the
overwintering of eggs as a mechanism restricting \emph{P.
pennsylvanica}to southern latitudes. Although identification of specific
mechanisms is beyond the scope of the present work, the importance of
temperature and fluctuating environmental conditions is consistent with
an extensive body of literature on thermal ecology in insects, including
a series of studies demonstrating fluctuating and temperature-dependent
selection in ambush bugs (Punzalan et al., 2008a; Punzalan, Rodd, \&
Rowe, 2008b, 2010).

Spatial errors in species distribution models

Errors in species occurrence data is virtually inevitable, resulting
from inaccuracies in georeferencing, imprecision in latitude and
longitude co-ordinates, or uncertainty in locality descriptions.
However, relative to other species distribution modelling methods based
on occurrence data, MAXENT has been found to maintain predictive
accuracy even with locational errors. MAXENT is also less sensitive to a
limited sample compared to other SDM (Wisz et al., 2008), and performs
well, so long as the data is comprised of widely distributed localities.
If, however, the subset of the data used is spatially biased, it may
exacerbate bias in estimates. Citizen science data may be particularly
prone to opportunistic collection, and hence, biased occurrence data.
However, in the present study, the contributions of citizen science data
limited to only three data points retained in any of the model; the
subset of the observations that temporally overlapped with the available
environmental data happened to be comprised mostly of museum data. The
data used in the models originated from multiple sources and databases
and consist of samples across much of the previously assumed range
of\emph{Phymata} .

A potentially more pressing concern is the potential errors arising from
species misidentification (i.e., misidentifying \emph{P.
pennsylvanica}individuals as \emph{P. americana} or vice versa), as it
may result in seemingly robust but inaccurate models (Lozier, Aniello,
\& Hickerson, 2009). MAXENT has been found to be robust against these
systematic biases relative to other SDM, but nevertheless, it is for
this reason that we attempted to remove data corresponding to
questionable identifications, which included a large portion of the
citizen science data. Future studies involving species distribution
modelling could surely benefit from the addition of citizen science data
(Tiago, Pereira, \& Capinha, 2017), as these databases improve and
provided that species occurrence data are sufficiently widely
distributed.

\subsection*{FIGURE LIST}

{\label{figure-list}}

Figure 1. Locations of \emph{P. americana} sightings in red and \emph{P.
pennsylvanica} sightings in blue; this includes both museum and citizen
science data and spans a time frame from 1864-2018.

Figure 2. Response curves of P. americana (top) and P. pennsylvanica
(bottom) to their strongest respective predictors. Red indicates the
mean response averaged over the 10 replicate MAXENT runs, while blue
indicates one standard deviation. For P. americana, BIO15 (Precipitation
Seasonality) had the largest percent contribution, while BIO3
(Isothermality) had the largest permutation importance. For P.
pennsylvanica, BIO17 (Precipitation of the Driest Quarter) had the
largest percent contribution, and BIO11 (Mean Temperature of the Coldest
Quarter) had the largest permutation importance.

Figure 3. Distributions of P. americana in red and P. pennsylvanica in
blue. The top row depicts predicted current (1970-2000) distributions
with points. The second and third row shows projected for 2050 and 2070
respectively at RCP8.5, with the lighter shades of blue and red indicate
previously suitable habitats becoming unsuitable, and the darker shade
of blue and red indicate previously unsuitable habitats becoming
suitable.

Table 1. The percentage of change of suitable habitat under different
RCP trajectories in comparison to current predicted distributions.

Figure 4. Projected overlap of distributions of P. americana and P.
pennsylvanica in 2050 at RCP2.6 (top left), RCP4.5 (top right), RCP6.0
(bottom left), and RCP8.5 (bottom right). Red indicates locations that
are only suitable for P. americana, blue indicates locations that are
only suitable for P. pennsylvanica, and purple indicates locations that
are suitable for both taxa.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

Figure 1. Locations of P. americana sightings in red and P.
pennsylvanica sightings in blue; this includes both museum and citizen
science data and spans a time frame from 1864-2018.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

Figure 2. Response curves of P. americana (top) and P. pennsylvanica
(bottom) to their strongest respective predictors. Red indicates the
mean response averaged over the 10 replicate MAXENT runs, while blue
indicates one standard deviation. For P. americana, BIO15 (Precipitation
Seasonality) had the largest percent contribution, while BIO3
(Isothermality) had the largest permutation importance. For P.
pennsylvanica, BIO17 (Precipitation of the Driest Quarter) had the
largest percent contribution, and BIO11 (Mean Temperature of the Coldest
Quarter) had the largest permutation importance.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

Figure 3. Distributions of P. americana in red and P. pennsylvanica in
blue. The top row depicts predicted current (1970-2000) distributions
with occurrence points. The second and third row shows projected for
2050 and 2070 respectively at RCP8.5, with the lighter shades of blue
and red indicate previously suitable habitats becoming unsuitable, and
the darker shade of blue and red indicate previously unsuitable habitats
becoming suitable.

Table 1. The percentage of change of suitable habitat under different
representative concentration pathway (RCP) trajectories in comparison to
current predicted distributions.\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
& \textbf{Year} & \textbf{RCP2.6} & \textbf{RCP4.5} & \textbf{RCP6.0} &
\textbf{RCP8.5}\tabularnewline
\midrule
\endhead
P. americana & 2050 & + 1.7\% & + 1.2\% & + 2.7\% & +
2.3\%\tabularnewline
& 2070 & + 1.6\% & + 1.1\% & + 2.0\% & + 0.7\%\tabularnewline
P. pennsylvanica & 2050 & + 0.5\% & 0 & + 0.8\% & + 0.4\%\tabularnewline
& 2070 & -- 0.4\% & -- 0.4\% & 0 & 0\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

Figure 4. Projected overlap of distributions of P. americana and P.
pennsylvanica in 2050 at RCP2.6 (top left), RCP4.5 (top right), RCP6.0
(bottom left), and RCP8.5 (bottom right). Red indicates locations that
are only suitable for P. americana, blue indicates locations that are
only suitable for P. pennsylvanica, and purple indicates locations that
are suitable for both taxa.

Data Accessibility:

The data used in this study have been uploaded to the Open Science
Framework, and can be accessed at the following doi:
https://doi.org/10.17605/OSF.IO/B5H3A.

\subsection*{REFERENCES}

{\label{references}}

Anderson, R. P. (2017). When and how should biotic interactions be
considered in models of species niches and distributions? \emph{Journal
of Biogeography} , \emph{44} (1), 8--17. doi: 10.1111/jbi.12825

Araujo, M. B., \& Pearson, R. G. (2005). Equilibrium of species'
distributions with climate. \emph{FORUM} , \emph{28} (5), 693--695.

Ara\selectlanguage{ngerman}újo, M. B., Pearson, R. G., Thuiller, W., \& Erhard, M. (2005).
Validation of species-climate impact models under climate
change.\emph{Global Change Biology} , \emph{11} (9), 1504--1513. doi:
10.1111/j.1365-2486.2005.01000.x

Arnold, S. J., \& Wade, M. J. (1984). On the Measurement of Natural and
Sexual Selection: Theory. \emph{Evolution} , \emph{38} (4), 709--719.
doi: 10.1111/j.1558-5646.1984.tb00344.x

Balduf, W. V. (1939). Food habits of phymata pennsylvanica Americana
melin (hemip.). \emph{The Canadian Entomologist} , \emph{71} (3),
66--74. doi: 10.4039/Ent7166-3

Balduf, W. V. (1941). Life History of Phymata Pennsylvanica Americana
Melin (Phymatidae, Hemiptera). \emph{Annals of the Entomological Society
of America} , \emph{34} (1), 204--214. doi: 10.1093/aesa/34.1.204

Bulgarella, M., Trewick, S. A., Minards, N. A., Jacobson, M. J., \&
Morgan-Richards, M. (2014). Shifting ranges of two tree weta species
(Hemideina spp.): Competitive exclusion and changing
climate.\emph{Journal of Biogeography} , \emph{41} (3), 524--535. doi:
10.1111/jbi.12224

Ceccarelli, S., \& Rabinovich, J. E. (2015). Global Climate Change
Effects on Venezuela's Vulnerability to Chagas Disease is Linked to the
Geographic Distribution of Five Triatomine Species. \emph{Journal of
Medical Entomology} , \emph{52} (6), 1333--1343. doi: 10.1093/jme/tjv119

Chen, I. C., Hill, J. K., Ohlemüller, R., Roy, D. B., \& Thomas, C. D.
(2011). Rapid range shifts of species associated with high levels of
climate warming. \emph{Science} , \emph{333} (6045), 1024--1026. doi:
10.1126/science.1206432

Darwell, C. T., \& Althoff, D. M. (2017). The relative contributions of
competition and abiotic tolerances in determining the geographical
distributions of four closely related Yucca species in
Texas.\emph{Journal of Biogeography} , \emph{44} (6), 1373--1382. doi:
10.1111/jbi.12907

Dowling, C. R. (2015). \emph{Using Maxent modeling to predict habitat of
mountain pine beetle in response to climate change} . doi:
10.1017/CBO9781107415324.004

Elith, J., Kearney, M., \& Phillips, S. (2010). The art of modelling
range-shifting species. \emph{Methods in Ecology and Evolution}
,\emph{1} (4), 330--342. doi: 10.1111/j.2041-210x.2010.00036.x

Elith, J., \& Leathwick, J. R. (2009). Species Distribution Models:
Ecological Explanation and Prediction Across Space and Time.\emph{Annual
Review of Ecology, Evolution, and Systematics} ,\emph{40} (1), 677--697.
doi: 10.1146/annurev.ecolsys.110308.120159

Fick, S. E., \& Hijmans, R. J. (2017). WorldClim 2: New 1-km spatial
resolution climate surfaces for global land areas. \emph{International
Journal of Climatology} , \emph{37} (12), 4302--4315. doi:
10.1002/joc.5086

Goldberg, E. E., \& Lande, R. (2007). Species' Borders and Dispersal
Barriers. \emph{The American Naturalist} , \emph{170} (2), 297--304.

Guisan, A., \& Thuiller, W. (2005). Predicting species distribution:
Offering more than simple habitat models. \emph{Ecology Letters}
,\emph{8} (9), 993--1009. doi: 10.1111/j.1461-0248.2005.00792.x

Guo, Q., Sax, D. F., Qian, H., \& Early, R. (2012). Latitudinal shifts
of introduced species: Possible causes and implications.\emph{Biological
Invasions} , \emph{14} (3), 547--556. doi: 10.1007/s10530-011-0094-8

Hickling, R., Roy, D. B., Hill, J. K., Fox, R., \& Thomas, C. D. (2006).
The distributions of a wide range of taxonomic groups are expanding
polewards. \emph{Global Change Biology} , \emph{12} (3), 450--455. doi:
10.1111/j.1365-2486.2006.01116.x

Hof, A. R., Jansson, R., \& Nilsson, C. (2012). How biotic interactions
may alter future predictions of species distributions: Future threats to
the persistence of the arctic fox in Fennoscandia. \emph{Diversity and
Distributions} , \emph{18} (6), 554--562. doi:
10.1111/j.1472-4642.2011.00876.x

Hutchinson, G. E. (1957). Concluding remarks. \emph{Cold Spring Harbor
Symposia on Quantitative Biology} , 415--427. doi:
http://dx.doi.org/10.1101/SQB.1957.022.01.039

Leach, K., Montgomery, W. I., \& Reid, N. (2016). Modelling the
influence of biotic factors on species distribution
patterns.\emph{Ecological Modelling} , \emph{337} , 96--106. doi:
10.1016/j.ecolmodel.2016.06.008

Loarie, S. R., Duffy, P. B., Hamilton, H., Asner, G. P., Field, C. B.,
\& Ackerly, D. D. (2009). The velocity of climate change. \emph{Nature}
,\emph{462} (7276), 1052--1055. doi: 10.1038/nature08649

Lozier, J. D., Aniello, P., \& Hickerson, M. J. (2009). Predicting the
distribution of Sasquatch in western North America: Anything goes with
ecological niche modelling. \emph{Journal of Biogeography} ,\emph{36}
(9), 1623--1627. doi: 10.1111/j.1365-2699.2009.02152.x

Malcolm, J. R., Markham, A., Neilson, R. P., \& Garaci, M. (2002).
Estimated migration rates under scenarios of global climate
change.\emph{Journal of Biogeography} , \emph{29} (7), 835--849. doi:
10.1046/j.1365-2699.2002.00702.x

Mason, L. G. (1976). Habitat and Phenetic Variation in Phymata Americana
Melin (Heteroptera: Phymatidae). II. Climate and Temporal Variation in
Color Pattern. \emph{Systematic Zoology} , \emph{25} (2), 123. doi:
10.2307/2412738

Masonick, P., Michael, A., Frankenberg, S., Rabitsch, W., \& Weirauch,
C. (2017). Molecular phylogenetics and biogeography of the ambush bugs
(Hemiptera: Reduviidae: Phymatinae). \emph{Molecular Phylogenetics and
Evolution} , \emph{114} , 225--233. doi: 10.1016/j.ympev.2017.06.010

Masonick, P., \& Weirauch, C. (2020). Integrative species delimitation
in Nearctic ambush bugs (Heteroptera: Reduviidae: Phymatinae): insights
from molecules, geometric morphometrics and ecological
associations.\emph{Systematic Entomology} , \emph{45} (1), 205--223.
doi: 10.1111/syen.12388

Menzel, A., Bairlein, F., Post, E., Beebee, T. J. C., Convey, P.,
Fromentin, J.-M., \ldots{} Parmesan, C. (2002). Ecological responses to
recent climate change. \emph{Nature} , \emph{416} (6879), 389--395. doi:
10.1038/416389a

Ning, S., Wei, J., \& Feng, J. (2017). Predicting the current potential
and future world wide distribution of the onion maggot, Delia antiqua
using maximum entropy ecological niche modeling. \emph{PLoS ONE}
,\emph{12} (2), 1--15. doi: 10.1371/journal.pone.0171190

Njiru, M., Mkumbo, O. C., \& van der Knaap, M. (2010). Some possible
factors leading to decline in fish species in Lake
Victoria.\emph{Aquatic Ecosystem Health \& Management} , \emph{13} (1),
3--10. doi: 10.1080/14634980903566253

Parmesan, C. (2006). Ecological and Evolutionary Responses to Recent
Climate Change. \emph{Annual Review of Ecology, Evolution, and
Systematics} , \emph{37} (1), 637--669. doi:
10.1146/annurev.ecolsys.37.091305.110100

Parra, J. L., \& Monahan, W. B. (2008). Variability in 20th century
climate change reconstructions and its consequences for predicting
geographic responses of California mammals. \emph{Global Change Biology}
, \emph{14} (10), 2215--2231. doi: 10.1111/j.1365-2486.2008.01649.x

Pfennig, K. S., Kelly, A. L., \& Pierce, A. A. (2016). Hybridization as
a facilitator of species range expansion. \emph{Proceedings of the Royal
Society B: Biological Sciences} , \emph{283} (1839), 20161329. doi:
10.1098/rspb.2016.1329

Punzalan, D., Rodd, F. H., \& Rowe, L. (2008a). Contemporary sexual
selection on sexually dimorphic traits in the ambush bug Phymata
americana. \emph{Behavioral Ecology} , \emph{19} (4), 860--870. doi:
10.1093/beheco/arn042

Punzalan, D., Rodd, F. H., \& Rowe, L. (2008b). Sexual selection
mediated by the thermoregulatory effects of male colour pattern in the
ambush bug \emph{Phymata americana} . \emph{Proceedings of the Royal
Society B: Biological Sciences} , \emph{275} (1634), 483--492. doi:
10.1098/rspb.2007.1585

Punzalan, D., Rodd, F. H., \& Rowe, L. (2010). Temporally Variable
Multivariate Sexual Selection on Sexually Dimorphic Traits in a Wild
Insect Population. \emph{The American Naturalist} , \emph{175} (4),
401--414. doi: 10.1086/650719

Punzalan, D., \& Rowe, L. (2015). Evolution of sexual dimorphism in
phenotypic covariance structure in \emph{Phymata} . \emph{Evolution}
,\emph{69} (6), 1597--1609. doi: 10.1111/evo.12680

Punzalan, D., \& Rowe, L. (2017). Hybridisation and lack of prezygotic
barriers between Phymata pennsylvanica and americana. \emph{Ecological
Entomology} , \emph{42} (2), 210--220. doi: 10.1111/een.12380

Rebelo, H., \& Jones, G. (2010). Ground validation of presence-only
modelling with rare species: A case study on barbastelles Barbastella
barbastellus (Chiroptera: Vespertilionidae). \emph{Journal of Applied
Ecology} , \emph{47} (2), 410--420. doi:
10.1111/j.1365-2664.2009.01765.x

Rebelo, H., Tarroso, P., \& Jones, G. (2010). Predicted impact of
climate change on European bats in relation to their biogeographic
patterns. \emph{Global Change Biology} , \emph{16} (2), 561--576. doi:
10.1111/j.1365-2486.2009.02021.x

Swanson, D. R. (2013). A Review of the Ambush Bugs (Heteroptera:
Reduviidae: Phymatinae) of Michigan: Identification and Additional
Considerations for Two Common Eastern Species. \emph{The Great Lakes
Entomologist} , \emph{46} (3--4), 154--164.

Taylor, S. A., Larson, E. L., \& Harrison, R. G. (2015). Hybrid Zones:
Windows on Climate Change. \emph{Trends in Ecology \& Evolution}
,\emph{30} (7), 398--406. doi: 10.1016/j.tree.2015.04.010

Tiago, P., Pereira, H. M., \& Capinha, C. (2017). Using citizen science
data to estimate climatic niches and species distributions. \emph{Basic
and Applied Ecology} , \emph{20} , 75--85. doi:
10.1016/j.baae.2017.04.001

Urbani, F., D'Alessandro, P., \& Biondi, M. (2017). Using Maximum
Entropy Modeling (MaxEnt) to predict future trends in the distribution
of high altitude endemic insects in response to climate
change.\emph{Bulletin of Insectology} , \emph{70} (2), 189--200.

Vallejo-Marín, M., \& Hiscock, S. J. (2016, September 23). Hybridization
and hybrid speciation under global change. \emph{The New Phytologist} ,
Vol. 211, pp. 1170--1187. doi: 10.1111/nph.14004

Wang, R., Li, Q., He, S., Liu, Y., Wang, M., \& Jiang, G. (2018).
Modeling and mapping the current and future distribution of Pseudomonas
syringae pv. Actinidiae under climate change in China. \emph{PLoS ONE}
,\emph{13} , 1--21.

Wisz, M. S., Hijmans, R. J., Li, J., Peterson, A. T., Graham, C. H.,
Guisan, A., \ldots{} Zimmermann, N. E. (2008). Effects of sample size on
the performance of species distribution models. \emph{Diversity and
Distributions} , \emph{14} (5), 763--773. doi:
10.1111/j.1472-4642.2008.00482.x

Yong, T. (2005). Prey Capture by a Generalist Predator on Flowering and
Nonflowering Ambush Sites: Are Inflorescences Higher Quality Hunting
Sites? \emph{Environmental Entomology} , \emph{34} (4), 969--976. doi:
10.1603/0046-225x-34.4.969

\subsection*{BIOSKETCH}

{\label{biosketch}}

VMZ and DP conceived of the study and contributed to the study design.
DP collected the \emph{Phymata} species occurrence data. VMZ collected
the environmental data, conducted the modelling work and created the
ArcMap figures. All authors contributed to the writing of the
manuscript.

\selectlanguage{english}
\FloatBarrier
\end{document}

