Climate change is likely to affect the distribution but not parapatry of the Brazilian marmoset monkeys (Callithrix spp.)

Parapatric distributional patterns can arise from abiotic or biotic factors, or from dispersal barriers. Climate change can potentially affect parapatry by changing species’ potential geographic distribution, and thereby widening or shrinking contact zones. Here, we study the effects of climate change on all six species in the genus Callithrix, a group of small‐sized Neotropical primates that is distributed parapatrically in eastern Brazil, allegedly due to biotic interactions.

Biotic interactions and climate conditions are the main causes of parapatry between closely related species (Bull, 1991;Bull & Possingham, 1995). In some cases, dispersal barriers can also contribute to this pattern. Amazonian rivers, for example, shape the parapatric distributions of sister species of birds and primates (Ayres & Clutton-Brock, 1992;Mercês, Lynch Alfaro, Ferreira, Harada, & Silva Júnior, 2015;Ribas, Aleixo, Nogueira, Miyaki, & Cracraft, 2011;Vale, Cohn-Haft, Bergen, & Pimm, 2008). Species are parapatrically distributed when their ranges are adjacent, mostly with nonoverlapping areas, with the exception of narrow sympatric zones where they have opportunities to interact (Bull, 1991). Two species that are ecologically similar can have their ranges limited by each other, resulting in a parapatric distribution in which species' persistence is given by competition strength, environmental heterogeneity and dispersal ability (Case & Taper, 2000;Case, Holt, McPeek, & Keitt, 2005). A similar scenario may arise from other interspecific interactions, such as parasitism or hybridization with low fitness offspring (Bergstrom, 1992;Gröning & Hochkirch, 2008;Ribeiro & Spielman, 1986). When parapatry is maintained by biotic factors, the presence of one species may play a major role in shaping the range of another, by preventing it from occupying regions where abiotic suitability coincides. Abiotic factors, however, may also drive parapatric distribution patterns. When distributional limits that establish parapatry are maintained by abiotic factors, the two species occupy virtually all their abiotically suitable regions, which happen to be both adjacent to one other and limited by abiotic factors (Martínez-Freiría, Sillero, Lizana, & Brito, 2008).
Climate change has the potential to converge distributions of closely related species to coincident geographic regions, leading to larger overlap of ranges and more intense interspecific interactions (Chunco, 2014;Engler, Rödder, Elle, Hochkirch, & Secondi, 2013;Jankowski, Robinson, & Levey, 2010;Krosby et al., 2015). If climate change displaces the distribution of two competitors to the same geographic area, for example, it is likely that only one species will be able to persist. Therefore, knowing the underlying causes of parapatry is essential to predict distributional responses to climate change on such species.
Here, we studied the effects of climate change on the parapatric distribution of Callithrix Erxleben, 1,777 species, a genus of small-sized Neotropical primates. Studies on the vulnerability of South America's biodiversity to climate change are particularly scarce in comparison with other continents (Pacifici et al., 2015;Vale, Alves, & Lorini, 2009), and primates, a very diverse and highly endangered order in South America (IUCN, 2017), are no exception.
Indeed, primates in Neotropical regions, such as Central America, Amazon and south-eastern Brazil, are among the most vulnerable to climate change (Graham, Matthews, & Turner, 2016). There are two studies focused on predicting likely future impacts of climate change on Neotropical primates, with results predicting future range contraction (Gouveia et al., 2016;Meyer, Pie, & Passos, 2014). This range contraction has also been predicted for other taxa in several biomes in the region, such as trees (e.g., Siqueira & Peterson, 2003), other vertebrates (Loyola, Lemes, Faleiro, Trindade-Filho, & Machado, 2012;Marini, Barbet-Massin, Lopes, & Jiguet, 2009;Souza, Lorini, Alves, Cordeiro, & Vale, 2011;Vale, Lorini, & Cerqueira, 2015) and invertebrates (Ferro et al., 2015;Giannini et al., 2012). If this trend of species' distribution contraction also holds for the genus Callithrix, depending on its location, it may lead to either (a) a reduction in the overlap between distributions, changing the genus's pattern of distribution from parapatric to allopatric; or (b) a convergence of species' distributions to some more climatically stable areas, moving the genus's distribution to a sympatric pattern in future.
Callithrix, family Callitrichidae ), includes six species, all endemic to Brazil, that exhibit a typical parapatric distribution pattern ( Figure 1). These tufted marmosets, F I G U R E 1 Original geographic range of Callithrix species within Brazilian biomes (Caatinga xeric shrubland, Cerrado savanna and Atlantic rain forest), according to the International Union for Conservation of Nature (IUCN, 2017) as popularly known, inhabit rain forests in the Brazilian Atlantic Forest, seasonal forests surrounded by savannas in Cerrado and xeric shrublands in Caatinga biomes (Figure 1). In the Atlantic Forest, Callithrix aurita (É. Geoffroy in Humboldt, 1812) and C. flaviceps (Thomas, 1903) occur in montane regions of the south-east, while C. geoffroyi (É. Geoffroy in Humboldt, 1812) and C. kuhlii (Coimbra-Filho, 1985) occur in the lowlands further north (Rylands, Coimbra-Filho, & Mittermeier, 2009; Figure 1). Callithrix jacchus (Linnaeus, 1758) typically occurs in Caatinga but also in some adjacent Atlantic Forest, while C. penicillata (É. Geoffroy, 1812) occurs in a broad area within the Cerrado and also in some adjacent Atlantic Forest areas Vivo, 1991). Recently, however, Callithrix penicillata and, particularly, C. jacchus have become invasive in south and south-eastern Brazil thanks to intentional or accidental human introduction (Begotti & Landesmann, 2008;Cerqueira, Marroig, & Finder, 1998;Modesto & Bergallo, 2008). Although rivers may not constitute effective barriers to the dispersal of Callithrix species (Rylands, Fonseca, Leite, & Mittermeier, 1996), some species have part of their distribution limited by rivers, such as the São Francisco river for Callithrix jacchus and C. penicillata (Malukiewicz, Hepp, Guschanski, & Stone, 2017), the Jequitinhonha river for C. geoffroyi and C. kuhlii, and the Contas river for C. kuhlii  Figure 1).
Despite fertile offspring, captive animal studies reveal a high mortality in the first years of hybrid infants (Coimbra-Filho & Mittermeier, 1973), suggesting lower fitness of hybrids. The recent diversification within the Callithrix genus, with an estimated appearance 2.5 million years ago (Perelman et al., 2011;Schneider et al., 2012) and the most recent split between Callithrix jacchus and C. penicillata dated to approximately 0.70 MYA (Malukiewicz et al., 2017), implies incomplete reproductive isolation, low phenotypic differentiation and consequent prezygotic isolation between species. Therefore, biotic factors such as hybridization producing low fitness offspring and interspecific competition due to high ecological similarity are likely to play a role in Callithrix's distributional patterns. At the same time, some contact zones are located at environmental transitions (Grelle & Cerqueira, 2006;Rylands et al., 1996), raising the possibility of abiotically driven parapatry, with range limits coinciding with areas of low environmental suitability for species' persistence. It is possible, therefore, that the parapatric distributional pattern of Callithrix species results from a combination of abiotic conditions (climate, topography and vegetation structure) and interspecific interactions. A thorough investigation of the issue, however, is still lacking.
In this study, we tested two working hypotheses for maintenance of parapatry, based on species' current distribution patterns: (a) under the abiotic-factor hypothesis, climate is the main cause of parapatry and species are able to occur in sympatry if their ranges overlap; and (b) under the biotic-factor hypothesis, interspecific interactions maintain range limits and one species excludes the other in overlapping areas. We then predicted the effects of future climate change on Callithrix species' distributions assuming that the causes of their parapatry will be maintained.

| ME THODS
To assess and anticipate impacts of climate change on Callithrix parapatry, we divided the analysis into three steps ( Figure 2): (a) we tested if climate is the likely cause of the parapatric pattern of the genus Callithrix; (b) generated ecological niche models (ENMs) for each species, in current and future scenarios, to map estimated climatically suitable areas for the occurrence of Callithrix species; and (c) based on the tests from the previous step, assumed the local exclusion in geographic space (in both current and future scenarios) if biotic hypothesis was supported. These methods have been employed before (e.g., Engler et al., 2013;Glor & Warren, 2011;Jovanović et al., 2018;Schiaffini, 2017;Suárez-Mota, Villaseñor, & López-Mata, 2015), but never to improve predictions for parapatric species.

| Species and environmental data
We compiled occurrence from the scientific literature and online databases: speciesLink (https://splink.cria.org.br/) and Global Biodiversity Information Facility (https://www.gbif.org/). As some of these species have been introduced in many areas of south-eastern Brazil (Supporting Information Figure S1.2), we considered only native-range occurrences in modelling procedures. We also checked every record for taxonomic and/or geographic errors, removing inconsistent records. In most cases, information about the date of the records was not available. We infer that the vast majority of records ranges from 1980s to the present, with few very old museum records from early 1900s. We gathered 123 occurrence records for Callithrix aurita, 57 for C. flaviceps, 84 for C. geoffroyi, 79 for C. jacchus, 83 for C. kuhlii and 137 for C. penicillata. To avoid the effects of spatial autocorrelation while retaining an adequate number of occurrence records for model calibration, we used spThin package in R environment, version 3.3.1 to select a subset of the data, with ≥10 km between records (Aiello-Lammens, Boria, Radosavljevic, Vilela, & Anderson, 2015;R Core Team, 2016).
For each species, we delimited a calibration area for ENMs based on the terrestrial ecoregions (Olson et al., 2001) in which species were present. These areas were defined by the limits of occupied ecoregions plus a 200 km buffer (ignoring areas where species had been recently introduced by humans). In total, we used four calibration areas, assuming these to be reasonable accessible areas

| Testing for the causes of parapatric distributional patterns
To elucidate the factors responsible for limiting species' ranges, we tested the degree of overlap between the climatic niches of  (1) Testing the explanation for the parapatric pattern of distribution involves the estimation and comparisons of climatic niches and expected null niches for a species pair, in which niche equivalency and niche similarity tests support the abiotic-factor (i.e., climate explains the parapatric pattern) or the biotic-factor hypothesis (i.e., climate is unable to explain the parapatry, and it is assumed that biotic interactions are likely to explain this pattern). Then, (2) ecological niche modelling is used to estimate the potential geographic distribution of the species. Finally, if the biotic-factor hypothesis is supported in the first step, we propose the use of (3) postmodelling processing in order to avoid the large overlapping distributions between species pairs. Otherwise, if the abiotic-factor hypothesis is supported, we propose the use of nonprocessed ENMs to estimate species' geographic distributions We characterized the environmental space available for the Callithrix genus as the first two orthogonal axes, given by a principal component analysis (PCA), based on the scores of bioclimatic variables from a set of 10,000 randomly located points throughout the three calibration areas together. We then extracted PCA scores of presence records for each species in order to estimate their climatic niches. Additionally, in order to carry out the similarity test we needed to characterize the environmental conditions available within each species' distribution. For each species, we extracted PCA scores from 1,000 randomly located points within the minimum convex polygon (MCP) defined by the species' presence records, plus a 100 km buffer.
We divided the environmental space into a grid of 100 × 100 cells to estimate the climatic niche for each species using a kernel smoother function for densities of occurrence (see Broennimann et al., 2012). We calculated the degree of niche similarity between each species pair using Schoener's D (D OBS ) as a metric of niche overlap (Rödder & Engler, 2011), which ranges from 0 (no niche overlap) to 1 (identical niches).
The niche equivalency test assesses whether the climatic niches of a given pair of species are identical (null hypothesis H 0 ) or not. Because equivalency tests do not adequately consider the availability of climate conditions within the areas accessible to each species (Peterson, 2011;Warren, Glor, & Turelli, 2010), we also used the niche similarity test. This test compares the observed similarity between the niches of two species with what is expected as the similarity between the observed niche of one species and densities of occurrence randomly generated within the climatic conditions available to the other species. For this purpose, we (a) generated 1,000 null densities of occurrence by randomly locating species 1's centre among the climatic conditions available (or background) for the species 2; (b) computed the D-value between these null densities of occurrence for species 2 and the climatic niche of species 1, producing a null distribution (D NULL ); and (c) tested whether the observed similarity between the climatic niches of the two species (D OBS ) is 95% greater than D NULL , using the probability (p) of obtaining D OBS given D NULL . Notice that each species pair had two tests (one in each direction): comparing the climatic niche of species 1 with null densities of occurrence of species 2, and vice versa. To simplify interpretation of the results, we analysed the probability of D OBS given both null densities. In this study, we were particularly interested in identifying species pairs for which climate cannot explain parapatry, so we tested the null hypothesis (H 0 ) that the species' climatic niches are as similar as would be expected by the differences between the available conditions at the geographic regions that each species occupies, using a one-tailed test (Warren, Glor, & Turelli, 2008).
Finally, we assumed that species pairs with equivalent climatic niches (i.e., niche equivalency test's H 0 not rejected) and/or those more similar than expected by chance (niche similarity test's H 0 rejected) are parapatrically distributed due to biotic factors (such as competition and hybridization). Otherwise, we assumed, (i.e., H 0 of equivalency test rejected and H 0 of similarity tests not rejected) that climate conditions are the best explanation for parapatry.
We ran 10 replicates for each algorithm, with 20% of data for model testing, using the cross-partitioning method (Hastie, Tibshirani, & Friedman, 2009) and fixed prevalence at 0.5. The significance of each replicate was tested by the area under the curve (AUC; Fielding & Bell, 1997), discarding those with AUC ≤0.7. Replicates that met this criterion were aggregated in a unique model for each algorithm using weighted-by-AUC mean. A final ensemble model was generated for each species, also by the weighted-by-AUC mean, of the four algorithm models.
ENMs were projected in the current and a future scenario for 2050 given the mean between both GCMs (HADGEM2-ES and MIROC5). Species' distributions were made binary (suitable vs. not suitable) using the least training presence as threshold value (Pearson, Raxworthy, Nakamura, & Peterson, 2007), a presenceonly method, due to our lack of true absence information. Final ensemble models were evaluated by the partial AUC (pAUC) by excluding the area where sensitivity equals zero (Peterson, Papeş, & Soberón, 2008) and using all records of the species as presences and records of the congeners as absences within the species' calibration area, because: (a) species are known to have parapatric distributions; and (b) if there was an inventory that recorded one species, it would have recorded congeners, if they were indeed present.

| Species distributions and postmodelling procedures
Because our main interest was in the shifts in current native distribution of the species, ENMs of each species were analysed only within their own calibration area. The geographic area delimited by the binary ensemble models (hereafter called untrimmed models) was interpreted as an estimate of the geographic distribution of those pairs of species whose parapatry is explained by the climate. However, if a species pair distribution was parapatric due to biotic factors, we assumed local exclusion between the species pairs and generated trimmed models for each one in order to better estimate the species' distribution.
These trimmed models were generated through a pixel-by-pixel analysis to prevent large overlaps between species' distributions.
First, we rescaled the predicted suitability values (from the ENMs analysis) within untrimmed models to range between 0 and 1, in order to make predictions comparable across species (Csergő et al., 2017). Then, within the overlapping area of a species pair, the species with the higher (rescaled) suitability score was considered as present in the pixel at the expense of the other. We recognized, however, that parapatric species generally have narrow sympatric zones in which species co-occur. We defined these zones based on the environmental suitability scores between the species: For low differences (<0.5), co-occurrence was allowed. When trimmed models were produced, their performances were compared with untrimmed models based on their sensitivity scores (i.e., the proportion of observed presences correctly predicted by models), specificity (i.e., the proportion of observed absences correctly predicted) and accuracy (i.e., the overall proportion of observed data correctly predicted).
We restricted four of the six species' geographic distribution assuming rivers as dispersal barriers. Although some southeastern Brazilian rivers likely act as dispersal filters rather than effective barriers, we chose a more conservative and realistic analysis based on species' current known distributions. Therefore, Callithrix jacchus' distribution was restricted to the north of São Francisco river and C. penicillata to the south of this river; C. kuhlii's distribution was limited between the Contas and Jequitinhonha rivers, and C. geoffroyi to the south of Jequitinhonha river.
Finally, to quantify the impacts of climate change on species' range, we calculated the area (km 2 ) of each species' distribution, as well as the minimum and maximum elevation values (metres above sea level) within this area, under current and future scenarios. Comparisons between current and future distributions were conducted based on trimmed models for species whose parapatry is maintained by biotic factors, or based on untrimmed models for those maintained by abiotic factors. We used the International Union for Conservation of Nature (IUCN) criteria B1b(i) to assess extinction risk associated with future climate change, which is based on a continuing projected decline in any of the extent of occurrence, considering the modelled distributions as species' extent of occurrence.
Under this criterion, a species that is projected to have a decline reaching an extent of occurrence of 100 km 2 qualified as Critically Endangered, 5,000 km 2 qualified as Endangered, and 20,000 km 2 qualified as Vulnerable (IUCN, 2012).

| Similarities between climatic niches
Environmental space of Callithrix species was summarized in two orthogonal axes that explained most of the data variation (PC1 = 39.5% and PC2 = 24.9%, Figure 3). PC1 was highly correlated with precipitation of warmest quarter (r = 0.84), temperature seasonality (r = 0.81) and mean temperature of warmest quarter (r = −0.72), while PC2 was correlated with mean diurnal range F I G U R E 3 Environmental space of Callithrix species summarized by the principal component analysis (PC1 = 39.51% and PC2 = 24.92%) using six bioclimatic variables. Colour crosses: species' occurrence records. Colour lines: convex hull including all occurrence records for each species. Grey crosses: climatic conditions for the area accessible to Callithrix, delimited by the three biomes where the genus is known to occur, and sampled through 10,000 randomly distributed points (r = 0.87), precipitation seasonality (r = 0.59) and annual precipitation (r = 0.44). In general, each Callithrix species occupied a specific area within the environmental space, with varying degrees of overlap When we took into account the available climatic conditions for species occurrence, niche similarity tests (Table 1) were unable to reject H 0 in any of the species pairs (i.e., the climatic niches of species pairs are as similar as expected by the climatic conditions available for them). Considering the results of both tests, we assumed that biotic factors play a major role in parapatry between three species pairs (Callithrix aurita and C. flaviceps, C. geoffroyi and C. penicillata, and C. flaviceps and C. penicillata), and overlap between their ENMs was restricted by the pixel-by-pixel analysis, producing trimmed models. By counting pixels within the overlap area between their ENMs, sympatry was assumed for 83% of the pixels between trimmed models of C. aurita (which prevailed in 11% of the pixels) and C. flaviceps (prevailed in 6%). For trimmed models of C. geoffroyi and C. penicillata, the former prevailed in 8% and the latter in 37% of pixels, while 55% were assumed as the sympatric zone. Almost all pixels were assumed within the sympatric zone (93%) for C. flaviceps and C. penicillata, which prevailed in 3% and 4% of pixels, respectively.

| Geographic distribution of Callithrix species
All final ensemble models showed good to excellent significance scores (pAUC > 0.83, Table 2). Two of the four trimmed models yielded more accurate results and higher specificity scores, but lower sensitivity scores, when compared with untrimmed models (Table 2). For future predictions, the two GCMs yielded slightly different predictions (Supporting Information Figure S1.1).
Under the climate change scenario, species tended to maintain most of their core current distribution area, with the exception of Callithrix flaviceps (Figure 4). On the borders of species' range, however, contractions and expansions are predicted, translating into a net reduction in distribution size for three species, and a net increase for four (Figure 4). Predicted percentage area reduction was larger, on average, than predicted area expansion. Callithrix jacchus, C. geoffroyi and C. kuhlii showed increases in area of 34%, 14% and 1%, respectively. Callithrix flaviceps showed 95% area reduction (1,505 km 2 in future), qualifying as Endangered under the IUCN, followed by C. penicillata (1,279,554 km 2 ) and C. aurita (187,625 km 2 ) with 37% and 27% reduction, respectively. In addition to Callithrix flaviceps' very large predicted reduction in geographic distribution size, we also predicted a reduction in its elevational distribution, which currently ranges from 42 to 2,079 m a.s.l. and is predicted to range from

TA B L E 1 Niche similarity test results
showing the probability of occurrence (p) of the observed similarity (D OBS ) between the species pair being significantly greater than expected by the null hypothesis (D NULL ). Only pairs of species with parapatric distributions were tested TA B L E 2 Significance of the ecological niche models given by the partial area under the curve (pAUC) scores for each species and comparisons between untrimmed and trimmed models based on their performance scores (i.e., sensitivity, specificity, and accuracy) (between C. kuhlii and C. penicillata). We predicted no changes in sympatric zone between Callithrix geoffroyi and C. kuhlii because the Jequitinhonha river was assumed as a dispersal barrier.

| D ISCUSS I ON
As discussed by many authors (Rylands, 1996;Rylands et al., , 1996Vivo, 1991), Callithrix species occur at distinct environments (e.g., phytophysiognomies, biomes, elevation and precipitation regimes). The montane species, Callithrix aurita and C. flaviceps, tend to occupy colder and more seasonal environments, while C. jacchus and C. penicillata have the broadest occurrence along environmental space (Figure 3), occurring in the warmer conditions available for the genus. This supports the idea that these latter species have the highest invasive potential within the group and, indeed, both are well known as invasive species in many areas in south-eastern and southern Brazil (Begotti & Landesmann, 2008;Cerqueira et al., 1998;Modesto & Bergallo, 2008;Pereira, Oliveira, & Ruiz-Miranda, 2008;Zago, Miranda, Neto, Santos, & Passos, 2013).
Alternatively, the invasive capacity of C. jacchus and C. penicillata may be explained by their persistence in highly seasonal habitats due to their ability to obtain secondary resources (i.e., exudates) in times of scarcity (Castro, 2003;Pinheiro & Pontes, 2015;Rylands & Faria, 1993), rather than climatic requirements. Additionally, Because most of Callithrix species do not have equivalent or more similar than expected climatic niches, we can assume that they have a certain degree of niche differentiation, leading to partially or nonoverlapping climatic requirements. Given the recent diversification of this taxa (Malukiewicz et al., 2017;Perelman et al., 2011;Schneider et al., 2012) and the current hypothesis that Callithrix diversification arose predominantly from the isolation of ancestral forms (Kinzey, 1982;Marroig, Cropp, & Cheverud, 2004), our results suggest a rapid climatic niche differentiation after species' isolation.
Niche similarities tests, however, can underestimate climatic niches (and consequently the overlap between them) because they are based on the current occupied distribution of organisms, which is potentially reduced by biotic factors and dispersal limitations.  (Traad, Leite, Weckerlin, & Trindade, 2012;Vilela & Del-Claro, 2011;Zago et al., 2013) may buffer the negative effects of climate change. Because variables possibly linked to Callithrix invasiveness were not included in the analyses, interpretations and inferences of our results are readily applicable to the species' native distribution, but not its invaded range.
The contrary is observed for others montane species (Beckage et al., 2008;Kelly & Goulden, 2008;Walther et al., 2002), and under future climate change, Callithrix aurita and C. flaviceps are not predicted to move to higher elevations. In the case of Callithrix flaviceps, the species will reduce its elevational distribution and, simultaneously, suffer severe losses of climatic suitability within its range.
This scenario is worrisome because Callithrix flaviceps is the most threatened species within the genus, already being classified as Endangered under the IUCN Red List (IUCN, 2017) due to its population size smaller than 2,500 mature individuals and ongoing decline (Rylands, Ferrari, & Mendes, 2008). The species is threatened by its  (Aximoff, Soares, Pissinatti, & Bueno, 2016;Carvalho et al., 2013;Nogueira et al., 2011).
Although considering the potential interference of congeners on species' native distributions, our study does not explicitly consider the implication of the two invasive species in the genus, Callithrix jacchus and C. penicillata, on the conservation of native species at the local scale. Additionally, the widespread deforestation of the Atlantic Forest (Ribeiro, Metzger, Martensen, Ponzoni, & Hirota, 2009) may act in synergy with future climate change, reducing the ability of species to adapt to that change (Hof, Levinsky, Araújo, & Rahbek, 2011). To top it all, the frequent and recent outbreaks of yellow fever in the Atlantic Forest (Moreno et al., 2011;Paules & Fauci, 2017) and its severe impacts on populations of primate species (Fialho et al., 2012;Holzmann et al., 2010;Moreno et al., 2015) worsen the situation of future conservation of endangered marmosets.
Contrary to our expectations, except for the Callithrix aurita-C. flaviceps pairing, we did not observe an overall contraction of climatically suitable area or reduction of sympatric zones for the species, probably maintaining the parapatric patterns of the genus in future scenarios. Furthermore, Krosby et al. (2015)  here are limited to pairwise comparisons. When more than two species' ranges overlap, therefore, our method is likely to underestimate more complex interactions due to increased uncertainties regarding the interactions' outputs within a multispecies context (Case et al., 2005).
Finally, our methodological approach provides the possibility of delimiting sympatric zones where species are predicted to cooccur. Here, sympatric zones were defined as the region with low difference in climatic suitability between species (arbitrarily defined as <0.5), but the size of these zones will depend not only on the extent of the area with low suitability difference, but also on the strength of interspecific interactions and dispersal capacity (Bull & Possingham, 1995;Case et al., 2005). For the genus Callithrix, the sympatric zones delimited by our approach are consistent with the known occurrences of native hybrids (see Figure 5). In the cases of species whose hybrids' occurrences are well known, the threshold value for the definition of sympatric zones can be estimated based on hybrids' presence (i.e., threshold value is the highest difference in climatic suitability in which hybrids occur). Although our method does not require a priori ecological information about species interactions and dispersion abilities, this information can and should be included whenever possible. Certainly, more detailed studies about the mechanisms behind species' parapatry (e.g., removal experiments; Cunningham, Rissler, & Apodaca, 2009) and sympatric zones could improve the predictions. In this study, we used a 0.5 threshold value for sympatric zones, but this value can be more or less restrictive according to available information or even be tuned by the occurrence data. For example, a coexistence threshold that optimizes the accuracy score, or any other performance metric, can be estimated by a sensitivity analysis using observed occurrences within the overlap area between models.
Here, we emphasized the potential outcomes of climate change on the parapatric distribution of closely related species when considering their biotic interactions. We presented a novel approach that requires only the data used in conventional ENMs and encourage its use to achieve better predictions of parapatric species' geographic distributions. In addition to its applicability to any species with sufficient data for conventional ENMs, our method is able to predict responses that are undetectable when using procedures that disregard the potential biotic interactions between parapatric species. Compared to the approach that considers only climatic factors, our approach represents a step forward in ENM methods by providing a spatially explicit way in which to address the biotic factors in species' distribution predictions, that is not data hungry and has a wide applicability.
When considering the biotic interactions in future projections, climate change is unlikely to affect the parapatric distribution patterns of Callithrix species or intensify these interactions by increasing the amount of distribution overlap. Callithrix species, however, will respond idiosyncratically to global changes: (a) a severe reduction in potential distribution area would maintain Callithrix flaviceps as Endangered; and (b) area loss for potential distributions for C. aurita and C. penicillata for the next 30 years; (c) slight and no increase in potential distribution area for C. geoffroyi and C. kuhlii, respectively; and (d) a substantial increase in C. jacchus' distribution. We strongly recommend the development of monitoring programs and conservation plans for Callithrix flaviceps, with a special focus on demographic changes and gene introgression with congeners.

ACK N OWLED G EM ENTS
We thank A.

DATA ACCE SS I B I LIT Y
The species' occurrence compiled in this study and the GIS layers of species' environmental suitability generated in this study are available at https://doi.org/10.6084/m9.figshare.c.4278242.
All data sources consulted in this study are properly cited in this article.