Edited by Rodolfo Dirzo, Stanford University, Stanford, CA, and approved December 10, 2015 (received for review July 2, 2015)
The tropical forests of South America are among the most diverse and unique habitats in the world in terms of plant and animal species. One of the most popular explanations for this diversity and endemism is the idea that forests retracted and fragmented during glacial periods, forming ecological refuges, surrounded by dry lands or savannas. These historically stable forest refuges would have been responsible for maintaining the pattern of diversity and endemism observed today. Here, we show that the Atlantic Forest of eastern South America probably expanded, rather than contracted, during the last glacial period. In addition, the emerged Brazilian continental shelf played a major, yet neglected, role on the evolution of this biodiversity hotspot during the last glacial period.
The forest refuge hypothesis (FRH) has long been a paradigm for explaining the extreme biological diversity of tropical forests. According to this hypothesis, forest retraction and fragmentation during glacial periods would have promoted reproductive isolation and consequently speciation in forest patches (ecological refuges) surrounded by open habitats. The recent use of paleoclimatic models of species and habitat distributions revitalized the FRH, not by considering refuges as the main drivers of allopatric speciation, but instead by suggesting that high contemporary diversity is associated with historically stable forest areas. However, the role of the emerged continental shelf on the Atlantic Forest biodiversity hotspot of eastern South America during glacial periods has been ignored in the literature. Here, we combined results of species distribution models with coalescent simulations based on DNA sequences to explore the congruence between scenarios of forest dynamics through time and the genetic structure of mammal species cooccurring in the central region of the Atlantic Forest. Contrary to the FRH predictions, we found more fragmentation of suitable habitats during the last interglacial (LIG) and the present than in the last glacial maximum (LGM), probably due to topography. We also detected expansion of suitable climatic conditions onto the emerged continental shelf during the LGM, which would have allowed forests and forest-adapted species to expand. The interplay of sea level and land distribution must have been crucial in the biogeographic history of the Atlantic Forest, and forest refuges played only a minor role, if any, in this biodiversity hotspot during glacial periods.
The extreme biological diversity of tropical forests has inspired and puzzled naturalists and scientists for centuries, and the forest refuge hypothesis (FRH) has long been one of the major paradigms to explain it. According to the FRH, forest retraction and fragmentation during glacial periods would have promoted isolation and consequently allopatric speciation in forest patches, or ecological refuges, surrounded by open habitats in the Amazon (1). Although originally based on climate fluctuations in the Pleistocene, the FRH was subsequently invoked for climate changes irrespective of the time period (2). The FRH was also applied to South America’s Atlantic Forest (3), one of the top-five biodiversity hotspots on Earth (4). The FRH gained broad acceptance during the 1980s when empirical paleoecological data from neotropical rainforests were still lacking. Nevertheless, heavy criticism came upon the FRH because some paleobotanical data showed that forests had persisted throughout glacial cycles (5). As paleoclimatic models of species and habitats became widely used, recent studies revitalized the FRH, not by considering refuges as the main drivers of allopatric speciation, but instead by suggesting that high contemporary diversity and endemism are associated with historically stable Atlantic Forest areas (6). This hypothesis is based on the assumption that populations were restricted to refugia within the bounds of the current observed distribution of a species, thus disregarding paleogeography. Here, we used coalescent simulations to test alternative demographic models and found potentially larger past distributions. Contrary to the refuge models, here, we show that forest specialist small mammal populations actually expanded during the last glacial period, after the expansion of the Atlantic Forest onto the Brazilian continental shelf, a process that has been neglected in paleoenvironmental reconstructions. This new idea is what we call the Atlantis Forest hypothesis, in homage to Plato’s legendary continent, which was eventually swallowed up by the sea.
We combined the results of distribution models with coalescent simulations based on DNA sequences to explore the congruence between scenarios of forest dynamics through time and the genetic diversity of mammal species cooccurring in the central region of the Atlantic Forest. First, we used climatic models to estimate range shifts from the last interglacial (LIG) [∼120,000 y before present (120 kybp)], through the last glacial maximum (LGM) (∼21 kybp), to the present. For LGM simulations, we extended the baseline climate data onto the exposed continental shelf, in contrast to other studies (6). Then, we used coalescent simulations to test the genetic data for congruence with historic demographic models of population fragmentation, bottlenecks, and expansions during the LGM, the LIG, and the penultimate glacial maximum, which corresponds to marine isotope stage 6 (MIS6) (∼160 kybp) (7). To generate and test alternative biogeographical hypotheses, herein we used data from five forest specialist small mammal species, whose occurrence is restricted to forested habitats. Considering their habitat fidelity (Table S1), these species should provide appropriate proxy records of forest changes through time.
Number and percentage of sites where the five small mammal species occurred within two distinct habitats (forests and open areas) and the calculated index of forest fidelity
The results of species distribution models were in general congruent among species, and three major trends emerged (Fig. 1). First, species showed suitable habitat displacement toward northern and lower elevation (i.e., warmer) areas during the LGM. Second, we found more fragmentation of suitable habitats during the LIG and the present than in the LGM. Third, we detected expansion of suitable climatic conditions onto the emerged continental shelf during the LGM, which would have allowed forests and forest-adapted species to expand as well. The first trend is a well-known biogeographic response to climate change, but the other two were unexpected because they are in the opposite direction to that predicted by the FRH, prompting us to hypothesize alternative biogeographic scenarios that incorporate forest expansion during glacial periods.
Habitat suitability of Atlantic Forest small mammals through time. Modeled ranges of small mammal species sorted according to their habitat fidelity (Table S1) were projected to the present, last glacial maximum (LGM) including the exposed continental shelf, and last interglacial (LIG). Habitat fidelity reflects the proportion of forest sites where each species occurs compared with open areas, ranging from 0 (lack of forest fidelity) to 1 (complete forest fidelity).
We used intraspecific Bayesian coalescent simulations to evaluate the fit of the empirical genetic data to distinct historical scenarios (Fig. S1 and Tables S2–S6). We first tested a single panmictic population against subdivision into two to five populations, either in a stable island system or including bottleneck or expansion during the LGM. Because all scenarios of population subdivision were rejected, we explored demographic changes in a single panmictic population. The stability scenario was rejected for all species, as well as bottlenecks [1–20% of the effective population size (Ne)] in any of the three events (LGM, LIG, or MIS6) (Fig. 2). Population expansions (1–20% of the Ne) during either LGM (which indicate population expansion onto the emerged continental shelf) or LIG, or both, were not rejected for all forest specialist species, usually confirming population expansion during both LIG and LGM (Fig. 2).
Species responses to demographic scenarios of forest refuges (red) and forest expansion (green) during climatic oscillations from 200 kybp. (A) Benthic d18O records reflecting changes in seawater temperature and ice volume (7), and corresponding demographic scenarios of population bottleneck (BN) and/or expansion (EX) in different periods: last glacial maximum (LGM), last interglacial (LIG), and penultimate glacial maximum (MIS6). (B) Results of the coalescent simulations for each species, confirming population expansions during both LGM and LIG. The horizontal line represents the threshold P value (0.05) for each species. Scenarios resulting in P values at or below this level were rejected (gray), and those above this level were not rejected (black). See Fig. S1 for additional scenarios tested and rejected.
Historical demographic models evaluated by coalescent simulations and corresponding parameters. For all models, we simulated populations that suffered three types of historical events: bottleneck, expansion, or population fragmentation. The prior distribution of time of events (t), migration (m), and of the demographic growth rate (r) are specified. Time is expressed in thousand of years (ky) before the present. LGM, last glacial maximum; LIG, last interglacial; MIS6, marine isotope stage 6 (penultimate glacial); Ne, effective population size.
Number of individual sequences (n), localities (loc), base pairs (bp), haplotypes (h), polymorphic sites (S), haplotype diversity (Hd) and nucleotide diversity (π), raggedness statistics (r), and deviation from neutrality tests (D, Fs, R2) for each small mammal species
GenBank accession numbers of the sequences used in estimating substitution rates of marsupials (Didelphimorphia) and rodents (Rodentia)
Fossil calibration points used in the substitution rate estimation of didelphid marsupials and sigmodontine rodents
Models of nucleotide substitution for each species as selected using MrModeltest
Results of the coalescent simulations in BayeSSC of the five small mammal species
The FRH posits that forest specialist species went through cyclical population subdivision and bottlenecks during past glacial periods and population stability or expansion during interglacials. Our results, however, support the alternative hypothesis that populations of forest specialist small mammals actually expanded during both LIG and LGM. We therefore suggest that the Atlantic Forest also expanded during both periods, especially onto the exposed continental shelf during the LGM. Global sea levels dropped to a maximum of −150 m during Quaternary Glacial Maxima, and the average Quaternary sea level was around −62 m (8). The southeastern South American coastline thus would have been displaced hundreds of kilometers to the east during glacial periods, exposing one of the widest shelves in the world, the Brazilian continental shelf (9) (Fig. 3). This expanded available surface area has not been taken into account in any regional paleoecological reconstruction to date. Rather, the area covered by Atlantic rainforest during the LGM has been estimated at only 150,000 km2 (excluding the continental shelf), which is nearly half of its current size (320,000 km2) (10). But the shelf represents an additional 270,000 km2 (9) of land available for potential expansion. A key question, therefore, is what kind of vegetation covered this shelf during the LGM. To address this question, we projected modeled climatic conditions for the Atlantic Forest and confirmed that extensive suitable areas of broadleaf evergreen rainforest were available on the shelf during this period (Fig. 3A), corroborating simulated LGM biome distributions based on global vegetation models (11). Analyses of a marine sediment core from the southeastern Brazilian continental margin revealed a wide spectrum of pollen types, including those representative of dense and semideciduous forests, coastal scrublands, and a low frequency of grasses (12). In contrast, another marine core from the same region suggested maximal grassland expansion and rainforest reduction during the LGM (13), with the two datasets thus reflecting regional variation in forest cover on the shelf. Other proxy data, such as changes in the oxygen isotopic composition of planktonic foraminifera from the same continental margin imply that maximum precipitation actually occurred around the LGM (14), when the shelf was exposed to its maximum extent. Paleochannels are also evidence of active continental sedimentary processes along the exposed shelf during low sea level periods (15), which point to the presence of riparian forests on the shelf, further supporting our results.
Connections between the Brazilian continental shelf and the Atlantic Forest. (A) Projected extent of suitable broadleaf evergreen rainforest areas (green) during the last glacial maximum (LGM) and the present. The LGM map shows the overlap of three models (dark green), two models (medium green), and one model (light green). Suitable areas for rainforest on the continental shelf during the LGM (green) are submersed in the present (blue). (B) Topographic map of the eastern Brazilian coast emphasizing the continental shelf (light blue, bounded by the −150 m isobath) and four key features mentioned in the text: the Abrolhos Bank (site 1), Tubarão Bight (site 2), Doce River (site 3), and Espírito Santo Mountains (site 4). The lowlands (<500 m) are represented in green, the highlands (>500 m) in orange, and ocean depth increases progressively toward darker blue areas.
Interglacial forest expansion on the continent is expected and confirmed by a long-term (140 kybp) pollen record from the Atlantic Forest (16). Forest extent on the continent during the LGM is, however, a controversial topic. Pollen records suggest that grasslands dominated south and southeastern Brazil during this period, indicating a markedly drier and cooler climate (17), and giving support to the FRH. In contrast, other pollen and diatom records combined with carbon and nitrogen isotopes suggest instead that cool, but humid, forests persisted, with no evidence of forest retraction throughout the LGM (18). Oxygen isotopes from speleothems, which are secondary carbonate cave formations, also indicate wetter climate from 70 to 17 kybp, especially during the LGM (19), thus supporting the Atlantis Forest hypothesis and refuting the FRH. We argue that, even if large forested areas were replaced by grassland and the tree line shifted 800–1,000 m downward on the continent during glacial periods (20), the Atlantic rainforest would have not been restricted to small patches as predicted by the FRH because forest cover would have extended onto the continental shelf (Fig. 3A).
The fragmentation of suitable habitats we found during both current and previous (LIG) interglacial periods is unexpected, given that a warmer and wetter climate is more suitable to forest expansion. This putative fragmentation might have been the result of topographic heterogeneity across eastern Brazil. Our species distribution models show that, during glacial periods, ranges are displaced to the north, where the coastal plain is very wide and the topography is monotonously flat, allowing for continuous species ranges. During interglacials, on the other hand, ranges moved south, where the steep coastal mountains (Serra do Mar) drastically replace a narrow coastal plain, leading to a very complex inland topography (Fig. 3B). The seaward side of the Serra do Mar has the highest mean annual rainfall (up to 3,600 mm) of the entire Atlantic Forest range whereas the inland side has typical seasonal climates (1,300–1,600 mm) (21). This heterogeneous topography is associated with distinct rainfall patterns and a mosaic of vegetation types, potentially leading to habitat and range fragmentation.
The Atlantic Forest shows strong latitudinal differentiation, and the geomorphological history of the continental shelf must have played an important, yet ignored, role in this process. A north–south floristic differentiation is caused by distinct temperature and rainfall regimes as mountain ranges are progressively farther from the coast and lower in altitude to the north of the Doce River (21) (Fig. 3B), which is a major biogeographic contact zone for many terrestrial vertebrate groups (22). Furthermore, an analysis of phylogeographic endemism of 25 vertebrate species indicated that this region was recently a location for a subdivision of the Atlantic Forest in two bioclimatic domains (23). Rivers and river valleys are frequently associated with biogeographic breaks in the Atlantic Forest (24), but contrasting results raised some questions about their roles as primary barriers whereas currently hidden geological or climatic barriers have been invoked (25). We argue that the Brazilian continental shelf is a putative barrier because it changes drastically at the Abrolhos Bank, from very narrow to the north to very wide to the south, where significantly more land was exposed during glacial periods (Fig. 3). Also, the steep slopes of the Tubarão Bight and the Espírito Santo Mountains form one of the narrowest lowland areas on the Brazilian coast (Fig. 3B). This region therefore marks major geomorphological and climatic transitions and must have been crucial in shaping the evolutionary history of the adjacent biota (26). Statistical phylogeographic studies on other Atlantic Forest organisms showed mixed results regarding demographic changes during the LGM, but moderate population growth and demographic stability were found for a warbler (27), toads (25), and one spider (28).
Eustatic changes play a major role on the evolution of islands and coastal habitats (29, 30), but such changes have been overlooked in the FRH debate. In addition, the vast continental shelf has been completely ignored in paleomodeling studies of the Atlantic Forest, despite its evolutionary significance in other tropical areas, such as Southeast Asia, where widespread rainforest covered the exposed Sunda shelf during the LGM (31, 32). These facts show that the predominantly glacial Quaternary environment is the norm, characterized by cooler climate, significantly lower sea levels, and slightly reduced precipitation, which was still enough to keep large stretches of rainforest instead of small refugia. This greatly expanded land area, resulting from the exposure of the continental shelf, supported widespread forests (31, 33). The current high sea level condition is therefore the exception, reflecting in reduced forest area in Southeast Asia and the Brazilian Atlantic coast. Our results show that the interplay of sea level and land distribution has been crucial in the biogeographic history of the Atlantic Forest and that forest refuges may have played only a minor role, if any, in this biodiversity hotspot during glacial periods.
Five species of small mammals, three didelphid marsupials [Marmosops incanus (Lund), Monodelphis americana (Müller), and Gracilinanus microtarsus (Wagner)] and two cricetid sigmodontine rodents [Euryoryzomys russatus (Wagner) and Delomys sublineatus (Thomas)], were included in this study. These small mammals are all forest specialist species, for which ecological data indicate a clear preference for forests compared with open habitats, and should respond to both current and past alterations in the distribution of the Atlantic Forest. All five species are widespread and locally abundant in most of the Atlantic Forest (34). To verify and quantify forest fidelity, we used occurrence data from previous studies (35, 36), obtained through standardized samplings (i.e., using the same type, number, and disposition of traps, and the same number of sampling days) in 54 sites in the Atlantic Plateau of the state of São Paulo, southeastern Brazil. The information on species occurrence is therefore comparable among the 54 sites. One third of the sites were located in continuously forested areas within one of the largest tracts of remaining Atlantic Forest. The other 36 sites were located in open, deforested areas used for agriculture (plantations of annual crops) within two 10,000-ha fragmented landscapes located adjacent to the sampled continuously forested areas (18 sites within a fragmented landscape containing only 10% of remaining forest, and 18 sites within a fragmented landscapes containing 50% of remaining forest). We computed the proportion of forest and open area sites where each of the five species occurred (PF and PO, respectively). We then used these proportions to calculate a simple index of forest fidelity (FD = PF − PO) that varies from −1 (preference for open habitats) to +1 (preference for forests), and for which a value close to zero indicates species equally well-distributed in both habitats. All five species were well distributed in continuously forested areas and were absent from most of the open area sites, and fidelity values ranges from 0.47 to 0.97 (Table S1). As expected from their high forest fidelity, these selected forest specialist species negatively respond to current deforestation and forest fragmentation whereas habitat generalist species do not respond or benefit from these processes (35, 37). The obtained values of forest fidelity corroborate several other ecological studies on habitat selection or preference for all five small mammal species (34).
Georeferenced occurrence localities of the five target species were gathered from data from museums, scientific articles, theses, and dissertations (38). The museum data were obtained from the online platforms SpeciesLink (splink.cria.org.br/), Global Biodiversity Information Facility (GBIF) (www.gbif.org), and Mammal Networked Information System (MaNIS) (manisnet.org/). The locations of the points of occurrence were verified using Google Maps (www.google.com.br/maps), Google Earth 7.0 (www.google.com/earth/), and SpeciesLink tools (geoLoc and infoXY), and possible mislabeled coordinates were excluded. Occurrence points for each species were then compared with distribution maps from the most recent taxonomic reviews, and discrepancies were further scrutinized to remove suspicious records. In total, we used 44 points of occurrence for E. russatus, 114 for M. incanus, 29 for M. americana, 17 for D. sublineatus, and 73 for G. microtarsus, scattered throughout the study area to create the suitability maps (Dataset S1). To estimate historical broadleaf evergreen rainforest occurrence, we used 395 points (Dataset S1) randomly extracted from the extent of this vegetation (39). To prevent spatial autocorrelation, only points distant at least 10 km of each other were included. To estimate habitat suitability for each species and habitat, we used 19 climate variables from WorldClim (40) in the 2.5-min resolution (www.worldclim.org/). For each species and habitat, correlated environment variables were rejected, as indicated by principal component analysis (PCA) and Mantel tests conducted in R 3.0.1 (41). Climatic layers were cropped to the Atlantic Forest river basins that encompass the distribution of the target species: Paraná, South Atlantic, Southeastern Atlantic, East Atlantic, and São Francisco. The algorithm chosen for modeling was that of maximum entropy, as implemented in MaxEnt 3.3.3 (42). Models were run with default regularization and 10 replicates subsampled, using 20% of the points for test and 80% for training each replicate. To determine whether the model's discrimination capacity is better than random chance, models were validated by accessing area under the curve (AUC), sensitivity, specificity, and accuracy values averaged across the replicates. The adequate environment conditions for each species and habitat in the present were projected into past layers by creating suitability maps for the last glacial maximum (LGM) period 21,000 y ago (including the continental shelf) and for the last interglacial (LIG) period 120,000 y ago. We used General Circulation Models (GCMs) for LIG and LGM estimated by the Community Climate System Model (CCSM). To estimate broadleaf evergreen rainforest occurrence during the LGM, we also used GCMs estimated by the Model for Interdisciplinary Research on Climate-Earth System Model (MIROC-ESM) and the Max-Planck-Institute Earth System Model running in low resolution grid and paleo mode (MPI-ESM-P). All layers of bioclimatic variables are available at the WorldClim website (www.worldclim.org/). We created binary maps of broadleaf evergreen rainforest by estimating the presence/absence of suitable areas using a 10th percentile training presence logistic threshold.
Genomic DNA was isolated from tissue samples (muscle or liver samples) preserved in ethanol (Dataset S2) using the salt extraction method (43). For each species, we sampled between 13 and 30 localities and between 1 and 15 specimens per locality. Polymerase chain reactions (PCRs) were performed to amplify 801 base pairs (bp) of the Cytochrome b gene (Cyt b), using primers MVZ05 and MVZ16 (44). PCR (12.5 μL) contained 1.25 μL of 10× buffer, 0.5 μL of MgCl2 (50 mM), 0.25 μL of dNTP solution (10 mM per nucleotide), 0.15 μL of each primer (10 mM), 0.15 μL of Taq Platinum (Invitrogen Life Technologies), and 1 μL of DNA (20 ng/μL). The PCR profile was as follows: initial denaturation at 94 °C for 5 min, followed by 39 cycles at 94 °C for 30 s, 48 °C for 45 s, 72 °C for 45 s, and a final extension at 72 °C for 5 min. PCR products were purified using ExoSAP enzymes (GE Healthcare Life Sciences), and the cycle-sequencing reactions were performed with a Big Dye v3.1 kit (Applied Biosystems Inc.), following the manufacturer's protocol. Samples were sequenced in both directions using automated DNA sequencer ABI 3500 (Thermo Fisher Scientific, Applied Biosystems), with the same primers listed above. Electropherograms were inspected, and sequences were aligned using the CLUSTAL W method implemented in MEGA 6 (45). All alignments were inspected and corrected manually. All generated sequences have been deposited in GenBank (Dataset S2).
Newly generated and GenBank sequences were used to estimate summary statistics for each target species. Numbers of haplotypes (h), polymorphic sites (S), haplotype diversity (Hd), and nucleotide diversity (Π) for each species were estimated with DnaSP v.5 (46). Deviation from neutrality was tested through different statistics: D (47), Fs (48), and R2 (49), mismatch distribution, and raggedness (50), which were calculated using DnaSP v.5. Coalescent simulations with 1,000 replicates were applied to determine the P value of each statistics. Significant P values (<0.05) were taken as evidence of demographic expansion (Table S2). These summary statistics were further used in coalescent simulations, to test the fit of each simulated scenario to observed data.
To estimate substitution rates for five taxa to be used in subsequent analyses, we used a fossil-calibrated approach (51). Estimates were performed independently for marsupials and for rodents. Both estimates were implemented in BEAST 2.1.3 (52) assuming a Yule speciation model and using a relaxed molecular clock with lognormal distribution. For marsupials, we used a Cyt b matrix containing previously published sequences (Table S3) and a published phylogeny (53) to set monophyletic constraints. We used two calibration points: the diversification of Didelphidae and the diversification of Didelphimorphia (Table S4). These calibration points were set as minimum and maximum ages, being the minimum age of the oldest unequivocal fossil belonging to clade (54). The maximum age of Didelphidae was based on stratigraphic bounding, which encompasses two stratigraphic units without fossil records of this clade (54). The maximum age of Didelphimorphia was based on phylogenetic bounding, which encompasses the age of the stem of didelphimorph (Peradectes) (54, 55). For rodents, we also used a Cyt b matrix containing previously published sequences (Table S3) and followed a published phylogeny (56) to set monophyletic constraints. We used six calibration points that represent multiple shallow (e.g., genera) and deep (e.g., tribes) calibrations of sigmodontines (56, 57) (Table S4). These calibration points were set as lognormal prior distributions, where fossil ages were set as offsets, and mean and SD were set as 1 (57). We performed two independent runs of 100,000,000 generations each and sampled every 10,000 generations. Analyses were performed in the CIPRES Science Gateway (https://www.phylo.org). Convergence of the Markov Chain Monte Carlo (MCMC) chains was verified in Tracer v.1.6 (58), making sure that all effective sample size (ESS) values were higher than 200. The resulting trees were combined using LogCombiner, and the consensus tree was generated in TreeAnotator, after a burn-in of 10%. Consensus trees were visualized in FigTree 1.4. The substitution rates obtained for each species were further used to estimate time to the most recent common ancestor (tMRCA) and in the coalescent simulations.
The tMRCA was estimated independently for each of five taxa, using all available sequences of each species (Dataset S2) and two close outgroups. Before this estimation, we tested the molecular clock for each dataset by performing maximum likelihood (ML) analyses in MEGA 6. The ML values of the obtained topologies were compared with and without molecular clock constraints also using MEGA 6. We rejected the presence of a strict molecular clock for the five small mammal species and used relaxed molecular clock approaches in the all of the subsequent analyses. The time to the most recent common ancestor (tMRCA) for each species was estimated assuming a Yule speciation model. We used a relaxed molecular clock with lognormal distribution as indicated by ML tests, our estimation of the mean substitution rates as the clock rate, and evolution models as selected using MrModeltest (59) for each species (Table S5). To obtain absolute dates in this analysis, we used our estimation of the mean substitution rates for each species (Table S5). Estimates were implemented in BEAST 2.1.3 and ran in the CIPRES Science Gateway. Runs and MCMC chain convergence, as well as resulting consensus trees, were defined as described above. Mean divergence dates and 95% highest posterior densities (HPDs) were obtained from the deepest node of the species in the tree.
To test putative scenarios of population diversification, we implemented coalescent simulations in the Bayesian Serial SimCoal (BayeSSC) (60, 61), and we then evaluated the fit of the empirical genetic data to distinct historical models. To evaluate fitting of each simulated scenario to the empirical data (summary statistics: h, S, Π, Tajima’s D, and pairwise differences), we applied the statistical framework described in the literature (27, 62⇓⇓–65). We started with simple models as a baseline against which more complex ones were compared (66). First, we simulated a scenario of a single panmictic population without any demographic event throughout time, which was then compared with three population subdivision models, where the ancestral population was separated into two, three, and five subpopulations (Fig. S1). The subdivision date was based on the estimates of tMRCA of each species (Table S5). We then simulated bottlenecks and expansions during the last glacial maximum (LGM) [∼11,000–33,000 y before present (11–33 kybp)], the last interglacial (LIG) (∼109–130 kybp), and the penultimate glacial period, which corresponds to marine isotope stage 6 (MIS6) (∼130–185 kybp) (67, 68). Bottlenecks and expansions were first simulated for each period separately and were then conjugated in both LGM and LIG (Fig. S1). These demographic events were simulated to reduce or increase the population effective size (Ne) to 1–10% and 5–20% of the current size (under a uniform prior). To further test a scenario of broader expansion, we also simulated an increase of 30–50% (under a uniform prior) during LIG followed by an increase of 5–20% (under a uniform prior) during the LGM. Ne was estimated from theta (Θ), according to the following formula: Θ = 2μNe, where μ is the mutation rate previously obtained for each species (Table S5). Θ was estimated in MIGRATE-N 3.6 (69) using transition/transversion ratios for each species estimated in MEGA 6 (Table S5), a static heating scheme, and the Bayesian Inference option with the default search parameters. For population subdivision models, we assumed that Ne from each present population was 1/2, 1/3, and 1/5 of the species Ne. In these scenarios, we assumed a linear stepping-stone model, with migration (m) equal to 0.01. BayeSSC input files for each species included the estimated range of population effective size Ne (under a uniform prior), mutation rates, generation time of 1 y (70), and transition/transversion (T/T) ratios (Table S5). For each model, we ran 1,000 simulations, in which five simulated summary statistics were obtained (h, S, Π, Tajima’s D, and pairwise differences). These simulations were then tested against observed data for each species (71). For assessing the significance of our data, we followed the procedure described in the literature (62, 63), where each summary statistic was compared with values representing the empirical distribution of the statistic from simulation. We then compared the combined empirical likelihood of each summary statistic (Cobs) against a null distribution of C (62, 63) and rejected simulated scenarios when P ≤ 0.05 (Table S6).
We thank J. Justino, L. Zanchetta, J. Dalapicolla, J. Agrizzi, R. Guimarães, R. Duda, and the many undergraduate students for laboratory work. We thank J. L. Patton and A. B. Rylands for helpful discussions and advice; M. Nascimento, J. Guedes, V. Colombi, V. Vale, E. Guerra, R. Duda, and K. Alves for field assistance; and J. Hoppe and F. Barreto for assistance with distribution models. This study was funded by Fundação de Amparo à Pesquisa e Inovação do Espírito Santo (FAPES, Brazil) and Conselho Nacional de Desenvolvimento Científico e Tenológico (CNPq, Brazil).
Author contributions: Y.L.R.L., L.P.C., and R. Pardini designed research; Y.L.R.L., L.P.C., A.C.L., R.G.R., and R. Pardini performed research; Y.L.R.L., L.P.C., A.C.L., R.G.R., H.B.-F., V.F., R. Paresque, M.P., and R. Pardini contributed new reagents/analytic tools; A.C.L., R.G.R., and H.B.-F. analyzed data; and Y.L.R.L., L.P.C., A.C.L., R.G.R., H.B.-F., A.C.B., V.S.Q., and R. Pardini wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
Data deposition: The sequences reported in this paper have been deposited in the GenBank database. For a list of accession numbers, see Dataset S2.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1513062113/-/DCSupplemental.
