Phenotypic plasticity and the colonization of new habitats: a study of a colonial spider in the Chaco region and the Cerrado

In social animals, group prey capture could facilitate colonization of new areas with low resource availability. Parawixia bistriata is a colonial spider inhabiting seasonal dry forests and mesic habitats in South America. Individuals capture prey as a group, which allows individuals to broaden their foraging niche by incorporating large prey that cannot be subdued in solitary captures. P. bistriata exhibits two behavioural ecotypes a “dry” (plastic) ecotype which modifies individual’s tendency to capture prey in a group depending on food resources and a “wet” (fixed) ecotype, whose tendency to group prey capture is only modulated by the size of the prey but not by prey availability. By reconstructing the range expansion of the species using phylogeographic and species distribution modelling techniques, we indirectly examined whether group prey capture could have helped P. bistriata in colonization of low resource habitats. Based on cytochrome c oxidase subunit I gene genealogy, we found older populations in northern Cerrado in Brazil with more recent populations located further south in Dry and Humid Chaco in Argentina, with the latter being the most derived. Species distribution modelling for each ecotype suggests that suitable habitat for each ecotype started to overlap at some point during the Last Glacial Maximum (21 ky BP). These results suggest that P. bistriata expanded from northern Cerrado south to the Gran Chaco, being able to colonize mesic habitats at a later stage when individuals reached southern territories in the Chaco. This evidence is opposite to the idea that GPC facilitated P. bistriata colonization from mesic to harsher environments. However, plasticity in group prey capture could have been important to allow individuals to establish in mesic habitats by reducing the cost of group capture when under high resource levels.


Introduction
There are certain behavioural traits that can facilitate colonization of a new environment by a species (Sol et al. 2002;Duckworth and Badyaev 2007). Given that the success of a colonizer depends on its ability to survive and reproduce under the new ecological conditions, individuals should be preadapted to the new environment or be flexible enough to quickly respond in an adaptive manner to the new environmental conditions (e.g., behavioural flexibility, Wright et al. 2010;Webb et al. 2014;Beever et al. 2017).
In group living species, there could be specific traits that facilitate moving and establishing into a new habitat (Spinks et al. 2000;Lövy et al. 2012;Cornwallis et al. 2017). In these species, it is assumed that evolution of sociality was favoured under the environmental conditions of the habitat where the species originated. However, when social species inhabit multiple habitats with different environmental conditions (e.g., resource), it is interesting to know if specific traits could have helped in the colonization of new environments when the source and the new habitats differ in resource availability.
Spiders are a particularly interesting group to develop studies on social behaviour because it mainly comprises solitary and aggressive species. In social species, group behaviours can be important in terms of colonizing ability (Spinks et al. 2000;Lövy et al. 2012;Cornwallis et al. 2017). As opposed to social insects which can be considered reproductive societies (i.e., group characteristics are mainly driven by the reproductive function) social species of spiders are considered foraging societies (Whitehouse and Lubin 2005), although other benefits in reproductive and survival aspects are gained by grouping (Avilés and Guevara 2017;. Thus, foraging related traits could potentially help in the colonization of habitats in social spiders. It is generally assumed that high prey availability levels and large prey size are necessary for groups to occur (Powers and Avilés, 2007;Yip, et al. 2008;Avilés and Guevara 2017;Grinsted et al. 2020) and that spider groups behave as "foraging flocks" that obtain the benefits of a large capture when they group their capture webs (Rypstra 1989). This can explain why most social species of spiders are found in tropical humid areas where insect prey is very abundant (Riechert 1985;Majer et al. 2015). However, there are social spiders inhabiting habitats with different resource levels. In these cases, we could ask whether social related traits could have facilitated the colonization of the new environments.
In some social species of spiders found under low prey availability conditions such as semi-arid habitats, groups of individuals remain in a colony benefiting by an increase in prey availability as they can exploit microhabitats such as open spaces which serve as flying corridors for insect prey (Lubin 1974). In addition, by grouping their webs spiders also profit from an increase in prey capture (ricochet effect;Rypstra 1989;Uetz 1979;Rao 2009). Adding to the benefits of group living, group prey capture (GPC) is another trait expressed by some spider species. Through GPC individuals can capture prey of large size that cannot be subdued by a solitary individual (Avilés and Guevara 2017 and references therein). This results in an increase in the amplitude of their foraging niche and decreases the variability in food intake compared to solitary species Majer et al. 2018) allowing individuals to obtain more food despite the low prey levels found in those environments.
Parawixia bistriata (Rengger 1836) is a colonial orb weaving spider (Araneidae) from South America. It mainly inhabits semi-arid habitats such as the Cerrado in Brazil and the Dry Chaco within the Gran Chaco in Argentina and Paraguay. In addition, it is also found in mesic habitats such as the Humid Chaco. These habitats differ in prey availability levels and the species shows behavioural flexibility in the expression of GPC depending on resource levels. Individuals capture prey in a group when prey is larger than the spider (Fowler and Gobbi 1988;Fernández Campón 2007). GPC occurs when a prey landing in a resident's web attracts neighbours due to the vibrations produced. The larger the prey, the further the vibration travels and the more spiders are attracted to the resident's web. First, the resident spider signals the neighbours not to enter its web by shaking it with its legs. When group capture occurs, this signal is ignored, and the resident and the first neighbours that arrive at the web participate in the capture of the prey. Once the prey is subdued, other individuals join and feed communally. The number of individuals engaging in prey capture and feeding depends on the size of the prey and the habitat type (capture: 2-9 individuals; feeding: 2-22 individuals). Interactions during GPC can sometimes be aggressive (probably modulated by prey levels) (Fernández Campón 2007). In P. bistriata, there are different behavioural ecotypes. A "wet" ecotype typical of mesic environments which shows a fixed tendency to GPC that does not depend on local prey availability conditions but only on prey size, and a "dry" (flexible) ecotype which under low prey availability conditions shows a higher tendency towards GPC compared to high prey availability conditions (Fernández Campón 2008). Studies carried out on this species suggest that GPC is important for their survival in dry habitats, with low prey availability. The absence of flexibility dependent on prey abundance of the wet ecotype impairs their growth and survival as seen when individuals were translocated to dry habitats under low prey conditions (Fernández Campón 2005, 2008. Under these circumstances, wet ecotype individuals show a lower growth rate and fecundity (estimated on the number of eggs per sac) than individuals that can adjust their tendency to capture prey as a group ("dry" ecotype) (Fernández Campón 2010). Interactions with conspecifics during GPC can sometimes be aggressive, implying costs to individuals participating in GPC . Thus, plasticity in this behaviour may be an advantage in habitats with low and more variable prey availability, whereas lack of plasticity in GPC expression is not detrimental in habitats with high and more stable prey levels.
It is possible to gain insight into a possible role of GPC in the colonization of new habitats by the use of phylogeographic studies and paleoclimatic information to determine the trajectory of colonization of the species over evolutionary time. In this study we use this approach to indirectly examine possible pathways in which GPC could have contributed to the colonization history of P. bistriata. We examine two possible scenarios. One in which this species originated and evolved in habitats with high prey conditions and the expression of GPC favoured the colonization of habitats with harsher conditions; and the other, in which P. bistriata originated in low prey condition environments with the flexible expression of GPC and later colonized high prey habitats with a later loss of behavioural flexibility in GPC that might have been favoured by the high prey levels and higher costs of aggressive interactions during GPC or by stochastic events (Lubin 1974;Fernández Campón 2007;Yip et al. 2017;Quero et al. 2020). Following the "high prey level" hypothesis we would expect oldest populations of P. bistriata occurring in mesic habitats and more derived ones in semi-arid environments. Alternatively, under the second scenario, we would expect the populations from mesic habitats to be more derived than those from arid environments.

Sample collection, DNA extraction, amplification and sequencing
We collected adults of P. bistriata in most of the species' range, comprising 64 individuals (one per colony) from 16 locations in the Gran Chaco in Argentina and Paraguay and the Brazilian Cerrado (Table 1; Fig. 1). Opisthosomas and palps were stored in absolute ethanol for species confirmation and deposited in the Arachnology Collection of IADIZA (CAI) (Mendoza, Argentina). Genomic DNA was extracted from one or two legs using an adapted "salting out" protocol (Sunnucks and Hales 1996). We amplified the mitochondrial region of the cytochrome c oxidase subunit I (COI) with the primers C1-J-1718 y C1-N-2776 (Framenau et al. 2010). The PCR conditions followed were: a denaturation step of 3′ at 95 °C, 30 cycles of 94 °C for 1′, 50 °C for 1′ and 72 °C for 2′, and a final extension step at 72 °C for 2′. PCR reactions consisted of 1 μL of DNA, 0.2 μL of 5U Taq DNA Polymerase (Genbiotech), 2-2.8 μL of 25 mM MgCl2 (Genbiotech), 2.5 μL of 10× KCl Buffer (Genbiotech), 0.5-1 μL of 10 mM dNTP mix (GE Healthcare), 0.5 μL of each primer 10 μM, and ultrapure H2O until completing a volume of 25 μL. The PCR products were purified and bi-directionally sequenced with the Sanger method, by means of the Sequencing Service of "Unidad de Genómica de INTA-Castelar" (Buenos Aires, Argentina).

Alignment, genetic diversity and haplotype reconstruction
The sequences were aligned in Geneious® 9.1.8 and were visually inspected; we found the best fitting substitution model according to the AIC on jMODELTEST 2.0 (Darriba et al. 2012). Diversity indices were calculated on Arlequin 3.5 (Excoffier and Lischer 2010). Relationships among haplotypes were visualized on median-joining networks (Bandelt et al. 1999) built on POPART (http://popar t.otago .ac.nz).

Population structure
We estimated the genetic distances among populations in Arlequin 3.5 (Excoffier and Lischer 2010). To verify the presence of isolation-by-distance among our samples, we analysed the correlation between the populations' genetic and geographic distances through Mantel tests. We assessed the population structure between biomes Cerrado and Chaco with F ST values calculated on Arlequin 3.5. We also performed an analysis of molecular variance (AMOVA) using (1) Cerrado and Chaco (dry and humid combined), and (2) Cerrado, Dry Chaco and Humid Chaco as separate groups, to test the structure both between and within these biomes. Population structure was also assessed on BAPS 6.0 (Corander et al. 2008), which determines the most likely number of clusters (k) within a given group of sequences. We allowed k to vary between 1 and 20.

Phylogenetic inferences and divergence times
We conducted a one locus *Beast analysis (Heled and Drummond 2010) in BEAST 1.8.0 to estimate the divergence among major mitochondrial lineages (as determined by BAPS) taking into account incomplete lineage sorting. We used a strict clock with a substitution rate of 0.0112 (SD = 0.001) substitutions/site/million years. Bidegaray-Batista and Arnedo (2011) estimated this rate for the Dysderidae family based on the well-resolved geochronology of the Mediterranean basin. Recently, Kuntner et al. (2013) found that the substitution rates estimated for the orbicularian families (which include Araneidae) overlapped with those from the Dysderidae, allowing these rates to be implemented to estimate the divergence times for orbicularian taxa. Two runs were conducted with 150 million generations, sampling every 15,000 generations. The resulting trees and log files from each of the two runs were combined using the computer program LogCombiner v1.

Paleoclimatic models
We employed Species Distribution Modelling (SDM) to reconstruct the spatial distribution of P. bistriata's dry and wet ecotypes in the past. SDM is a statistical tool that uses the association between different environmental variables and known locations of species presence to define the abiotic conditions within which populations can be maintained, assuming actual conditions were similar in the past (niche conservationism) (Guisan and Thuillier 2005). Because behavioural ecotypes are associated with different habitats (environmental conditions) and we wanted to use SDM as an independent approach for testing directional expansion of P. bistriata between habitat types, we modelled the distribution of each ecotype separately, based on the actual occurrence of populations in dry (Cerrado and Dry Chaco) and wet (Humid Chaco) habitats. We used the presence data collected in field trips, Levi (1992), museum holding databases (Museo Argentino de Ciencias Naturales-MACN and Museo Nacional de Historia Natural de Paraguay-MNHNP), and checked presence data from the web databases iNaturalist (www.inatu ralis t.org) and Sistema Nacional de Datos Biológicos (SNDB; https ://datos .sndb.mincy t.gob.ar) by F. Fernández Campón.
To avoid spatial autocorrelation in the output of the model we only used those records that were 1.5 degrees (150 km) from each other. This was a prerequisite for the model to be robust and reliable when projected to past times. That left us with 19 data points for the dry and six occurrence points for the wet ecotype.
As a first step, we used actual environmental information and species spatial occurrences to fit a model that explain the potential spatial distribution of P. bistriata ecotypes in the present. As second step, the model previously fitted was re-projected to past by using environmental information of Last Glacial Maximum (LGM) and Last Interglacial period (LIG) to reconstruct species potential distribution for those periods of time. The environmental information considered was only climatic because climate has a strong influence on invertebrate physiology and distribution ranges (Brandt et al. 2020;Kellermann et al. 2012;Addo-Bediako et al. 2000) and there exist available climatic data on actual and past conditions. Climate layers at a resolution of 2.5 arc-min were obtained from the World-Clim 1.4 free climate database http://www.world clim.org/), which comprises values for 19 bioclimatic (bc) variables, averaging the 1950-2000 period (Hijmans et al. 2005). Prior to the modelling process, we performed Pearson's correlation analysis on variables in order to avoid spatial correlation (Pearson > 0.90) and reduce over-parametrization. As a result, eight bioclimatic predictors were used: (bc2) mean monthly temperature range; (bc4) temperature seasonality; (bc5) maximum temperature of warmest month; (bc9) mean temperature of driest quarter; (bc13) precipitation of wettest month; (bc14) precipitation of driest month; (bc15) precipitation seasonality; (bc18) precipitation of warmest quarter. Paleoclimate layers representing the LGM (21 ky BP) and LIG (130 ky BP) were those used to reproject to past the species spatial distribution. Paleoclimatic layers were downloaded from the WorldClim website (http://www.world clim.org/), which include downscaled climate data from different Global Climate Models (GCMs), based on original data made available by CMIP5 (Coupled Model Intercomparison Project Phase 5; http://cmip-pcmdi .llnl. gov/cmip5 /); data were calibrated using WorldClim 1.4 as baseline'current' climate.   Table 1 Models were built using the maximum-entropy algorithm MaxEnt version 3.3.3k of the software (http://www.cs. princeton.edu/ ~ schapire/maxent). We selected MaxEnt because it works well with relatively small sample sizes (Hernandez et al. 2006(Hernandez et al. , 2008Wisz et al. 2008;Tognelli et al. 2009). As MaxEnt selects variable features and protects against overfitting internally through regularization (Quinn et al. 2018), we identified optimal levels of model complexity in the model using the R package ENMeval (Muscarella et al. 2014). ENMeval performs automated runs and evaluations of ecological niche models by selecting among different combination of scores of regularization and features to obtain models that maximize predictive ability and avoid overfitting (Muscarella et al. 2014). Specifically, we used a range of values for the regularization multiplier (β) (0.5-4 in increments of 0.5), and considered linear, quadratic, hinge, product threshold and categorical feature classes. Feature classes are functions created by MaxEnt for each environmental variable (Phillips and Dudík 2008) and the regularization multiplier controls the smoothness of the distribution curve. As occurrence data on P. bistriata were relatively low, we ran modelling with leave-one-out data partition as recommended by Pearson et al. (2007) and used 10.000 random background points. Finally, we selected the best fitting model to that with the lowest AIC value.

Genetic diversity
TN93 + G was the best fitting model for nucleotide substitution (Tamura and Nei 1993). The dataset presented 29 haplotypes (Table A.1). Most haplotypes were restricted to one locality, except some which were shared among locations within Chaco (Fig. 1a). Diversity levels were very low (Table 2). Six out of 16 populations were monomorphic, two in Cerrado and four in Chaco (Table A.1).

Genetic structure
Haplotypes were structured into four groups according to BAPS (Fig. 1b; Table A.2): three groups within the Cerrado and one group comprising all Chaco locations plus a locality in the Cerrado that borders the Pantanal (Bo). BAPS groups correspond to clades defined in the phylogenetic tree, with ancestral clades including Cerrado populations (Fig. 2, see  below). The haplotype network (Fig. 3) can be divided into the four groups inferred by BAPs. In the network, group four (Chaco + Bonito) shows a star-like topology, with an haplotype (H_19) from LY and R34 locations in Dry Chaco occupying a central position, and many low frequency haplotypes emerging from it (i.e. younger ones). Haplotypes from MB, RG and Suma locations in the Humid Chaco occupy derived positions within the star, indicating their younger age relative to haplotypes from Dry Chaco. Fst value between biomes showed important structuring (Table 2). AMOVA analyses using biomes Cerrado and Chaco as geographical groups supported the spatial structure detected by BAPs and the Fst: while a third of the variability is explained by differences between biomes indicating strong structure at that level, similar percentages of variability were shown among population within biomes and within populations (Table 3a). In addition, when considering Dry and Humid Chaco as well as the Cerrado (Table 3b), the AMOVA showed that variability among populations within biomes was important (39.85%), however, a significant percentage of variability was found among the three biomes (28.25%). The Mantel tests did not indicate significant correlation between genetic and geographical distances (r = 0.095, p = 0.241 from 1000 randomizations) implying an absence of isolation by distance.  ( Bayesian gene tree inferred for P. bistriata using COI dataset. Colours of branches indicate posterior probabilities according to the reference shown at the left of the tree. Node age is shown at each node with the horizontal bar indicating 95% HPD. First coloured lines at the right of the populations sampled show the biome in which they were found. Red: Cerrado; yellow: Dry Chaco; green: Humid Chaco. Second coloured line and numbers show individuals within each BAPS group, colour of groups as in Fig. 1. Numbers at the nodes indicate both posterior probability/node age

Phylogenetic inferences and divergence times
The gene tree obtained from the *Beast analysis (Fig. 2) recovered groups similar to those inferred by BAPS (Table A.2): groups 1 and 2 from Cerrado both resolve as monophyletic, with highest posterior probabilities; group 3 -the most southern group from Cerradoresolves as paraphyletic, including group 4 from Chaco as a derived and monophyletic lineage. Haplotypes from group 1, located at northern Cerrado can be considered the most ancient, branching off first from the last common ancestor at approximately 5,47 million years before present (My BP) (95% highest posterior density, HPD 6.30-3.41). Conversely, haplotypes belonging to the Humid Chaco emerges as the most derived cluster within  Table A.1. Circle size is proportional to frequency, colours correspond to populations and red dots represent missing intermediate haplotypes group 4, supported by a high posterior probability. The clade including all haplotypes from Chaco populations and Bonito last shared a common ancestor with a clade from Cerrado that would have existed at about 1.27 My BP (95% HPD 1.5-0.77). This could be considered the time of divergence between spiders from both biomes, Chaco and Cerrado. Unlike clades located in Cerrado, some clades within Chaco (which in BAPS form one same cluster) show weaker support values, one exception is the lineage from the Humid Chaco. According to this analysis this latter clade is relatively young, with a time of divergence estimated in 31 ky BP (95% HPD 26-3). An animation of the inferred geographic distribution changes of P. bistriata through time can be viewed in http://mland is.githu b.io/phylo wood/ by opening the txt matrix in Supplementary Materials.

Paleoclimatic models (species distribution model)
The best fitted model was lqhpt 1.5 and had an AIC = 727.87 with a delta AIC of 2.48 compared to the second-best model. It included seven parameters, excluding (bc5) maximum temperature of warmest month from the overall set of environmental variables considered. The model for the LIG (130 ky BP) scenario shows disjointed suitable areas for the dry ecotype occurrence (Fig. 4a). They are in north-western Argentina and Paraguay and eastern Brazil. There are also smaller areas in the north west of the continent in Venezuela and north-eastern Brazil. During the LGM (21 ky BP), results from both MIROC and CCSM models mostly agree (CCSM simulation is shown Supplementary Material A.4). They show an expansion of suitable habitat with the previous disjointed areas now in contact and forming a large area covering much of the Cerrado, the Chaco and further to the north west ( Fig. 4b). At present times, the location of suitable habitat is similar to that of the LIG but more expanded (Fig. 4c). Suitable areas for occurrence of the wet ecotype during LIG are in the Amazon region and north-western Cerrado and do not overlap with those areas of the dry ecotype, especially in Brazil (Fig. 4a). There is another smaller area in Venezuela showing some overlap with suitable areas of occurrence of the dry ecotype and another area with low probability of occurrence (< 0.5 in south east of Brazil and Paraguay). During the LGM the two models show an important reduction of suitable areas for the wet ecotype restricted to what it is now part of the Humid Chaco (in north-eastern Argentina and Paraguay) with some overlap with the dry ecotype in very reduced parts. There is also another small area in northern Brazil in what it is the Atlantic Forest today (Fig. 4b). The current scenario shows one area located in the Humid Chaco similar to that in the LIG model but with a larger extension. Suitable areas for the occurrence of dry and wet ecotypes show more overlap at present times (Fig. 4C).

Discussion
Using indirect approaches by reconstructing the range expansion of P. bistriata and the paleo and current distribution of ecological niches of dry and wet ecotypes, we provided information that serves to indirectly examine the relevance of expression of GPC on colonization of new habitats. We evaluated two contrasting hypotheses. One hypothesis was based on the premise that high prey levels are necessary to allow grouping of individuals (by increasing tolerance among conspecifics). Under these conditions GPC occurred, and colonization of harsher environments was facilitated through GPC as P. bistriata can broaden their foraging niche and survive despite lower prey conditions thanks to the expression of this behaviour.
Contrary to what is expected following the "high prey level" hypothesis, commonly suggested for the evolution of sociality in spiders, we found that P. bistriata's oldest populations are located in north central Brazil within a "dry" biome ( Fig. 2). They later expanded to the Dry Chaco finally colonizing the Humid Chaco. Furthermore, distribution models for both ecotypes support this finding. When modelling environmentally suitable areas for the occurrence of the dry and wet ecotypes, we found that environmental conditions for both ecotypes start to show some overlap during the Last Glacial Maximum Scales at the right of each map indicate the probability of occurrence of the species in a specific area (LGM) at approximately 21 ky BP. At some point between the Last Interglacial period (LIG; ~ 129-116 ky BP) and the LGM, environmental conditions for both ecotypes started to overlap making it possible for the dry ecotype to colonize wet biomes (Fig 4). Divergence times estimated in the COI genealogy show that the clade containing populations of the Humid Chaco diverged from one of the branches corresponding to Dry Chaco populations around 310 ky BP and started to diverge within the Humid Chaco approximately 180 ky ago. These data are not in agreement with the hypothesis that sociality and GPC in P. bistriata originated in mesic habitats and colonization of harsher environments such as the Dry Chaco was facilitated by the expression of group prey capture and feeding.
What is supported by the molecular data is the hypothesis that P. bistriata originated in dry environments and colonization of mesic habitats (Humid Chaco) occurred in the last stages of their range expansion. Divergence times based solely on gene genealogies usually overestimate populations divergence times (Edwards and Beerli 2000;McCormack et al. 2011) and these estimates based on P. bistriata's COI genealogy have large confidence intervals. Taking this into account and despite estimates are not identical, data inferred using both approaches give support to the idea that P. bistriata originated in dry biomes with colonization of wet biomes occurring later in the expansion of the range of the species. We can infer that P. bistriata's ancestor was found in northern South America and ancestral populations of this species expanded from the northern Cerrado towards the south until entering the Chaco region. This suggests that behavioural flexibility or plasticity in GPC is an ancestral trait and this plasticity was later lost when the species colonized and established in mesic habitats with higher levels and less variability in prey abundance. Plasticity in a trait that is relevant for a species survival is expected to be favoured under variable conditions such as the Cerrado and Dry Chaco, both with marked dry and wet seasons. P. bistriata is the only social species so far described within the genus. This is a Neotropical genus with most species found in the Amazone area, with some species present in central America and up to the humid tropical regions of southern Mexico (Levi 1992). In semi-arid and seasonal habitats such as the Cerrado and the Chaco there are patches of open areas that serve as corridors for flying insects. P. bistriata is found in these microenvironments and at the border of forest patches not within the forest canopy. It seems clear that one of the advantages of grouping in this species is the possibility of exploiting this "prey-rich" microenvironments unavailable for solitary species. In addition, although GPC may have enabled P. bistriata to expand its foraging niche and reduce variability in food intake, individuals in the colony still seem to suffer from low prey level as seen by the variability in body condition which is stronger than in wet habitats populations (Fernández Campón 2005).
Shrinking and expansion of suitable habitat during glacial periods, such as expansion of the rain forests towards Cerrado and Chaco (Sobral-Souza et al. 2015;Bartoleti et al. 2017;Trujillo-Arias et al. 2017) may have contributed to the structuring observed in the mitochondrial dataset and may have selected for individuals expressing plasticity in GPC. The results agree with the "plasticity first" hypothesis (Levis and Pfennig 2016) which postulates that while ancestral populations are expected to be plastic, these plastic responses might become genetically accommodated or assimilated in derived populations (Stein and Bell 2019). That is, there is a change in the extent of plasticity due to a shift in the environment. This may result in the absence of plasticity in the trait when selective pressures in the new environment do not favour it. This has been shown in some studies (e.g., diet-induced plasticity in feeding regimes in spadefoot toad tadpoles becoming genetically assimilated as tadpoles colonized new niches; Levis et al. 2018)." Translocations studies suggest that differences in levels of plasticity in GPC in P. bistriata have a genetic basis (Fernández Campón 2008), and thus could be under selection. There is evidence from other studies in which plasticity has been adaptive in extreme novel environments (Wang and Althoff 2019). Based on COI genealogy, populations from Dry habitats in Cerrado and Chaco show higher divergence among them than among populations from Dry and Humid Chaco although the latter differ in GPC plasticity levels. Because COI is assumed to be neutral to selection, we expect the gene genealogy to show divergence based on the time populations have been isolated due to distance and not due to differences in selective forces. Divergence times between populations of dry and humid environments do not seem to be long enough to show in neutral markers such as COI. Finding genes associated to fitness traits will allow us to thoroughly test the adaptive value of plasticity in GPC under different environmental conditions (Tong et al. 2020). So far, the indirect data available seem to go along those lines. P. bistriata show plasticity in other foraging related traits that seem to confer higher fitness to individuals expressing it. Populations studied at the Cerrado, present plasticity in web architecture depending on different type of prey (Sandoval 1994). Thus, it will be interesting to carry out studies to examine whether differences in plasticity also occur in traits related to other aspects of its biology (e.g., physiology).
Our study provides insights into the mechanisms underlying adaptation and diversification in new environments. Phenotypic plasticity in the ancestral populations may have played a role in colonization of new habitats, triggering phenotypic changes that allow the colonizers to deal with the novel conditions (Levis and Pfennig 2016;Stein and Bell 2019). Results from this study show that plasticity in GPC was present in older populations before colonization of humid habitats. This indirect evidence suggests that GPC may have allowed P. bistriata to establish in habitat with different environmental conditions, however, more direct evidence is needed to support this hypothesis.