Phylogeny and Time Divergence estimation
Our IQTree ML phylogenetic reconstruction, as well as our BEAST time calibrated tree, were completely concordant with the clades previously recovered by Sánchez Herrera et al. [13](Fig 1, Fig S1), though our current dataset includes twelve additional species. The family Polythoridae separated from our Outgroups ~56 Ma (70 - 48 Ma 95% HPD [High Probability Density intervals], pp [posterior probability] = 1). Our time divergence reconstruction suggests that the MRCA (Most Recent Common Ancestor) of the family Polythoridae diverged ~ 42 Ma (48 - 35.11 Ma 95% HPD, pp = 0.992) during the Eocene epoch. The Chalcopteryx –Chalcothore crown clade itself began to diversify later during the late Oligocene ~24 Ma (35.7- 12.52 Ma 95% HPD, pp = 0.999), while the crown group of all of the other Polythoridae genera (Cora s.s., Miocora, Euthore and Polythore clades) began to diversify ~33.14 Ma (41.52 - 25.52 Ma 95% HPD, pp = 1). Chalcothore separated from Chalcopteryx ~ 23.49 Ma (35.73-12.51Ma 95% HPD, pp = 0.999) around the late Oligocene and early Miocene boundary, however the crown group of Chalcopteryx only diversified during the late Miocene ~ 9.62 Ma (16.97 - 4.325 Ma 95% HPD, pp = 0.997). The MRCA of the Cora s.s clade separates from all the other taxa ~33.14 Ma, however its crown group of Cora s.s. appears ~17 Ma (24.24 - 9.06 Ma 95% HPD, pp = 1) around the end of the early-Miocene (Fig 1). The crown group of the Miocora clade appears ~ 14.5 Ma (22.5 - 6.78 Ma 95% HPD, pp = 1), a few million years after Cora s.s during the mid-Miocene, although the MRCA of Miocora (separate from the Euthore and Polythore clade) appears ~ 24.75 Ma (32.18 - 17.28 Ma 95% HPD, pp = 1) during the late-Oligocene or early Miocene epochs (Fig 1). The two sister clades of Euthore and Polythore diverged from each other around the same epoch as Miocora, ~22.28 Ma (29.37 - 15.70 Ma 95% HPD, pp = 1), but both form crown groups concurrently around the early Miocene, though Euthore estimates are slightly older than those for Polythore (Fig 1). Euthore crown group diversified ~18.39 Ma (25.06 - 12.17 Ma 95% HPD, pp = 1) while Polythore seems to be ~16.68 Ma (22.80 - 10.46 Ma 95% HPD, pp = 1). Within the Euthore clade we recovered two sister clades: Euthore sensu stricto (e.i. E. fasciata fasciata, E. fasciata fastigiata, E. fasciata plagiata, E. fassli, E. sp. nov), and what used to be part of Cora (i.e. E. lugubris and E. klenei). Within the Euthore clade the divergence times recovered suggest that Euthore sensu stricto, ~6.62 Ma (10.72 - 3.09 Ma, 95% HPD; pp = 1), is significantly younger that the previous classified Cora species within this clade, ~12.8 Ma (19.04 - 6.13 Ma, 95% HPD; pp = 0.978).
In the highly speciose Polythore clade we recovered the divergence times for all the previous geographical clades reported in Sanchez Herrera et al[13] (Fig 1, Fig S1). The crown group of the Amazonian clade appeared ~8.72 Ma (14.36 - 3.24 Ma 95% HPD, pp = 1), while the Andean clade was ~10.57 Ma (15.43 - 6.41 Ma 95% HPD, pp = 1). Within the Andean clade, the crown groups of the West and East Andean clades appeared ~ 2.98 Ma (5.34 - 1.19 Ma 95% HPD, pp = 1) and ~7.16 Ma (10.78 - 4.18 Ma 95% HPD, pp = 1), respectively. Moreover, estimates for the North (~4.21 Ma, 6.35 - 2.35 Ma 95% HPD, pp = 1) and South Eastern (~3.25 Ma, 5.57 - 1.64 Ma 95% HPD, pp = 1) crown groups suggest they appeared around the Pliocene - Pleistocene epochs (Fig 1).
Biogeographical Reconstruction
We designed and tested five different scenarios—Control (S0), Pebas & Acre Systems (P&AS, S1), ANDES Old uplift (ANDES O, S2), ANDES young uplift (Andes Y, S3), Pebas & Old Andes (P&OA, S4) and Pebas & Young Andes (P&YA, S5)—to reconstruct the Ancestral Areas through multiple biogeographical models (e.g. DIVA, DEC or BayArea) for the Polythoridae. The designated areas include; (A) Central America, (B) Tumbes-Choco-Magdalena, (C) North Western Andes, (D) North Eastern Andes, (E) Venezuela Highlands, (F) Central Andes, (G) Amazon Basin, (H) Guiana Shield, and (I) Brazilian Shield. For the control scenario, we tested if the phylogenetic relationships (i.e. cladogram) alone will recover the Ancestral Areas without imposing any restrictions (e.g. dispersal and adjacency) of the geological events including the Andes uplift and Marine incursions of South America (see details in Fig 2). We selected the best biogeographical models for all scenarios (S0, S1, S2 , etc.) using the AICc criteria from BioGeoBEARS. For most of the tested scenarios, the DEC+J was the best and the least fitted was the BayArea model (Table S2). However, following the Ree and Sanmartín [15] critique of the DEC+J model, we also reconstructed the next best fitted biogeographic model based on the AICc. For the Control (S0) and P&AS (S1) we reconstructed the ancestral areas using the DIVALike model, while for the scenarios involving the Andes (S2, S3, S4 and S5) we selected the DEC model using both BioGeoBEARS and RASP software. The lowest AICc values across all the scenarios was for the P&AS (S1), and it was followed very closely by our control S0 (Table 1). The next lowest AICc values were the ANDES Young (S3) and the P&YA (S5), the others (S2 and S4) seem not be fitting the diversification model. However, there were no major differences in the reconstructed ancestral areas (Table S2) when comparing among all five models in both BioGeoBEARS and RASP. To simplify, here we will describe the most likely ancestral areas estimated for the four scenarios with the lowest AICc values (Fig 3). All reconstructions became problematic at the deeper nodes, so the likelihood values supporting the areas are lower than 70% in nodes 81, 80 and 79 (Fig 3), suggesting uncertainty in the ancestral areas for the MRCA for the family Polythoridae (node 81), for the clade that includes Cora s.s, Miocora, Euthore and Polythore (node 80), and for the MRCA of the Miocora, Euthore and Polythore clades (node 79). For the MRCA (node 78) of Euthore and Polythore, the Northwest Andes (C) is the most commonly predicted ancestral area for all scenarios (Fig 3). The MRCA (node 56) for the crown group of Euthore is consistently found in the Northwest Andes (C) in all scenarios, however for the S0 and S1 the Northeastern Andes (D) is also supported (Fig 3). The MRCA (node 47) for the crown group of Cora s.s. has an ancestral area in the Northwest Andes (C), the Central Andes (F) and Venezuelan highlands (E); the latter was not supported in S1 (Fig 3). The MRCA (node 77) of the Polythore crown group has ancestral areas in the Northwest Andes (C), and the Amazon Basin (G), and for S3 it also supports the Central Andes (F) (Fig 3). The MRCA (node 49) for the Miocora crown group is found in the Tumbes-Chocó-Magdalena Valley (B) and the Northwest Andes (C). Finally, the youngest MRCA (node 43) of the Chalcopteryx-Chalcothore crown group showed as more plausible areas the Amazon Basin (G), and the Venezuelan Highlands (E); S2 also supports the Brazilian Shield (I) (Fig. 3, Table S2). It is interesting to note that the different scenarios (S0-S5) within our BioGeoBEARS and RASP analyses did not produce fundamentally different ancestral areas; inclusion of the formation of the Pebas and Acre wetlands systems, as well as mountain building events, did not alter our results in comparison with a simple control scenario that allowed free dispersal between adjoining regions at different periods of geological time (Fig. 2). Our best selected DIVALIKE model implemented in RASP for the Andes scenario suggests several dispersal (24), vicariance (11) and extinction (3) events within the different genera of Polythoridae (Figure S3); many of these involve movement within and between the different ranges of the Andes, as well as movement into the Amazon, Guiana Shield, Venezuela Highlands, the Tumbes-Chocó-Magdalena Valley, and Central America.
Diversification Analyses
To investigate the patterns of diversification and extinction rate variation through time and across lineages, we chose RevBayes[16, 17] to test the best diversification models and rate shifts, due to the high uncertainty associated with phylogenetic tree estimation. The Episodic Birth-Death with multiple time periods was the best model explaining the diversification pattern in this group of damselflies (Table S4). The estimated net diversification (λ-μ) and speciation (λ) rates show an increase, however the relative extinction (λ/μ) and extinction (μ) rates pattern seem to behave constantly through time (see Figure S5). The Branch specific model detected a shift in diversification (λ-μ) rate corresponding to the Andean Clade in Polythore with the highest shift in the Eastern Clade (Fig 4). When we observed the relative extinction (λ/μ) rate across the tree, we observed that this rate also decreases for this clade (Fig 4). Our estimates show that the other clades within this family have a more constant diversification (λ-μ) and relative extinction (λ/μ) rates across the branches of the estimated tree.