Calling the Banners: exploring the complex biogeographical history of the Neotropical banner-wing damselflies (Odonata: Polythoridae)

Sorting out the biogeographical history of species in the New World Tropics is challenging. We here examine the roles of evolution and biogeography in driving the distribution and diversification of the damselflies of the family Polythoridae. This family comprises seven genera with a total of 57 species, distributed across much of Central and South America. Through phylogenetic analysis, relaxed-clock molecular dating and biogeographical analysis we find these genera originated ~56 Ma and started diversifying ~45 Ma. As with other neotropical groups, the Polythoridae have a primary origin in the Northern Andes; at least one genus first appears in the Amazon Basin. Diversification rates have been uniform in all genera except one—Polythore—where a significant increase in the late Pliocene (~3 mya) correlates with mountain building. While our molecular clock suggests correlations with some major geographical events, our biogeographical modeling (with BioGeoBEARS and RASP) found little influence of the formation of the Pebas and Acre systems or Andean mountain building, possibly due to the short branch lengths in our time-dated phylogeny.

Introduction group are generalist predators in both their larval and adult forms, the types of aquatic habitats in which their odonate larvae are found vary, with some species using lakes and ponds, while others inhabit streams and rivers. The larvae of Polythoridae prefer moving water; most species are found in small streams in mountain regions, though some additionally exist within slow streams in the Amazon Basin 11 . These species display a diverse range of colors and patterns on both the wing and body within and among genera; previous work by Sánchez Herrera et al. 12,13 has shown that for some genera, wing colors may be polymorphic, raising the question of what factors drive this color diversity.
We explore hypotheses concerning the diversity within Polythoridae, investigating the influence of the changing biogeography on this family, and its ability to explain the distribution and phylogenetic relationships we see today. Using relaxed-clock molecular methods we estimate divergence times within and among the genera of this family. We also investigate their biogeographical patterns of diversification in relation to the mountain building events of the Andes, as well as the extent of marine and freshwater environment incursions in different periods that might have influenced distribution and speciation in these genera. Specifically, we evaluated the following questions: (1) Have the Andes and/or the Amazon Basin acted as the centers of origin and diversification of the family Polythoridae? (2) Are there multiple interchanges among Neotropical regions, particularly between Amazon and Andean regions and the Andean and Central America?
(3) How has Polythoridae diversified (e.g. constant, episodic) through evolutionary time? (4) Are there significant diversification shifts associated with these biogeographical events that can indicate an adaptive radiation for the genera within the family? Overall, this work represents the first continental-scale biogeographical analysis of a neotropical family of damselflies. With their reliance on aquatic habitats, damselflies are likely to respond very differently to the landscape than plants, vertebrates, or terrestrial insects; our analysis allows us to consider the formation of diversity in the modern neotropics through an unexplored perspective.

Phylogeny and Time Divergence
Our phylogenetic reconstruction was generally consistent with the clades previously recovered by Sánchez Herrera et al. 13 (Fig 1) from all the other taxa ~33.14 Ma, however its diversification estimates to be ~17 Ma (24.24 -9.06 Ma, pp = 1) around the end of the early-Miocene (Fig 1). The Miocora clade diversification took place ~ 14.5 Ma (22.5 -6.78 Ma, pp = 1), a few million years after Cora s.s during the mid-Miocene, although it separates from the Euthore and Polythore clade ~ 24.75 Ma (32.18 -17.28 Ma 95% CI, 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, pp = 1), but both start their diversification concurrently around the early Miocene, though Euthore estimates are slightly older than those for Polythore (Fig 1) (Fig 1). Our lineage through time plot suggest that the major accumulation of lineages for Polythoridae occurred during the Miocene epoch (Fig 1). Our lineage through time plot also shows that most of the species diversity within the family was acquired during the later Miocene, Pliocene and Pleistocene epochs (Fig 1).

Biogeographical Reconstruction
For the three scenarios (Control, P&AS, and ANDES) we evaluated, the best biogeographical model selected by BioGeoBEARS was DEC+J and the least fitted was the BayArea model (corrected Akaike Information Criterion (AICc) values > 245, Supplementary Table 1). However, following the Ree and Sanmartín 14 critique of the DEC+J model, we also reconstructed the next best fitted biogeographic model based on the AICc, using both BioGeoBEARS and RASP, for all the scenarios (Table 1). Comparing across the different BioGeoBEARS and RASP reconstructions there were no major differences among our three scenarios, however we found several disparities among the ancestral area estimation performed by BioGeoBEARS and RASP within each scenario (first row in Table 1).  Table 1). The Cora s.s. clade highlighted a similar issue across scenarios; here the BioGeoBEARs models consistently support the Northwest Andes (C) as an ancestral area, although depending on the model ( Fig. 2, Table 1) it can be accompanied by the Northeast Andes (D) and/or the Amazon Basin (G) and the Venezuela Highlands (E). For the RASP reconstructions of Cora s.s., the Northwest Andes (C) again is the most commonly predicted ancestral area, followed by the Northeast Andes (D) and/or the Central Andes (F). For the Miocora, Euthore s.l. and Polythore clades specifically, the BioGeoBEARs DIVALIKE model showed the greatest inconsistency in comparison with the other models (Fig. 2, Table 1 Table 1). It is interesting to note that the different scenarios (Control, P&AS, ANDES) within our BioGeoBEARS and RASP analyses, which accounted for many of the major geographical events in South America over the last 50 million years, did not produce fundamentally different results; 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. This may be due to the relatively short branch lengths within the tree for most of the species groups, especially within Polythore.
Conversely, any long branch lengths within the tree may also suggest the potential for undetected extinction events early in the formation of some of these genera. Our S-DEC model implemented in RASP suggest several dispersal (29), vicariance (10) and extinction

Diversification Analyses
The Episodic Birth-Death with multiple time periods was the best model explaining the diversification pattern in this group of damselflies ( Table 2). 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 Supplementary Figure S3). 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 3). When we observed the relative extinction (λ/μ) rate across the tree, we observed that this rate also decreases for this clade (Fig 3). 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.

Discussion
Overall, our results suggest that Polythoridae originated around the Eocene, however most of the diversification within the family appeared during the Miocene, Pliocene and Pleistocene epochs (Fig 1). During these periods there were a number of major biogeographical events occurring, including the uplift of different regions of the Andes, as well as the formation of the Pebas and Acre Systems. Our biogeographical analyses suggest strong support for an origin-or at least early diversification-in the Northeastern and Northwestern Andes for most genera in Polythoridae. Chalcopteryx clade, the earliest of the modern genera to emerge in our tree, may have had an Amazonian origin (Fig 2, Supplementary Figure S4). Chalcopteryx clade separated during the marine condition retreat period in South America ~ 42 Ma; this basal node is consistent with recent studies that suggest Amazonia as the primary source of Neotropical diversity 15 . Cora sensu stricto has an origin in the Northeastern and Northwestern Andes (Supplementary Figure   S4) during the period when the sub-Andean river system was forming; Miocora appears later in this period with an origin in the Northwestern Andes (Supplementary Figure S4).  Figure S4). This potential for an Amazonian origin in Polythore is likely driven by a basal clade of the genus, containing Polythore aurora, and P. mutata, which is currently found in the Amazon. These arose as a separate lineage in the mid-Miocene, not long after the estimated formation of the Pebas System.
The other Amazonian Polythore ( P. manua and P. spaeteri) appear to have evolved very recently, in the mid-Pleistocene, and are found within the Amazon Basin relatively near the Andes; they emerge from within a clade found mainly in the Central Andes. Taken as a whole an origin in the Northern Andes, with a subsequent distribution to the Amazon and Central Andes, seems the most likely scenario for Polythore.
An origin in the Northern Andes is consistent with the analysis of several other neotropical groups, which suggests that mountain building in these regions was an important impetus generally for neotropical diversity 9,16 . This appears to also be the case with Polythoridae, as the genera not only have estimated ancestral areas associated with the Andes, but the timing of origination and diversification also coincides with periods of Andean uplift (Fig   2). Diversification in many of the genera of Polythoridae took place much later than their origin within the tree, leading to long branch lengths at the base of the tree, and many very short branches near the tips. While Chalcopteryx originates during the period or marine retreat, diversification in this group doesn't take place until the development of the Acre system. Likewise, Euthore, with its origin during the early Pre-Pebas, does not begin to diversify until the late Pre-Pebas. Cora, Miocora and Polythore diversify during the Acre, with most diversity in Polythore developing during the Modern River System era, much of it during the Pleistocene (Fig 2). Many of these later periods of diversification are also associated with periods of intensified Andean uplift, despite our biogeographical modeling not finding definitive support for an Andean influence on species diversity (Fig 2, Supplementary Figure S5).
Our analyses suggest that Polythoridae has been diversifying through an episodic pattern (Table 2). Speciation (λ) and net diversification (λ-μ) show an overall increase through time, while relative extinction (λ⁄μ) and extinction (μ) rates seem to remain somewhat constant through their evolutionary history. However, our branch-specific model of diversification shows a significant shift for the Andean Clade within the genus Polythore which might suggest that at least for this genus the intensified Andean uplift during the Pliocene has been promoting speciation (Fig 3). Genera such as Stenocora (not included in this analysis) and Chalcothore have relatively limited distributions, while Chalcopteryx, Cora, Miocora, Euthore and Polythore are found broadly across different regions in Central and South America. This may be somewhat tied to species diversity within these genera: for example, Polythore and Euthore are species rich (21 and 14 species, respectively) while Chalcothore and Stenocora are each represented by a single species. Still, species richness does not immediately explain the distribution range of all these genera; Chalcopteryx, with a broad distribution across much of South America, comprises 5 species, while Cora s.s, with its habitat similarities to Polythore and Euthore, has only 9 species.
For much of the family Polythoridae, vicariance and dispersal, with species moving into new habitats as they become available, is a likely explanation for the current patterns of diversity. We found multiple interchanges among the Amazon and Andean regions; many of them involved 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. As noted above, Chalcopteryx may have originated in the Amazonian region and diversified through vicariance during the formation of the Acre System. Polythore shows a likely pattern of dispersal from the Andean regions to the Amazon in some clades, as well as movement from the Northern to Central Andes. As observed above, a basal node comprised of a Polythore clade dispersed early to the Amazon, while speciation into the Amazon from the Central Andes occurred more recently. Nodes separating the Northwestern Andean from Northeastern/Central Andean species likely represent a case of vicariance, as does the separating of the Amazonian clade. Some extinction is predicted to have been associated with the 'picta' node in Polythore (Supplementary Figure S2). On the other hand, Euthore moved from the Northwestern to Northeastern Andes and then to the Venezuelan highlands. This genus is currently found at locations from 1000-2000 meters in elevation; it is plausible that they colonized new habitats that became available during mountain building events in the Andes and Venezuelan highlands. There is also some likelihood of extinction in the Euthore 'proper' node. Cora s.s and Miocora diversification was more likely driven by dispersal, as it covers a range of landscape types, and has moved across the Isthmus of Panama, the only group in Polythoridae to do so. Interestingly, the age at which both of these genera diversified matches the age suggested by recent geological studies of the Central American Seaway closure during the middle of the Miocene 3 .
Beyond physical isolation, are there other factors likely involved in the diversification in these groups during the Andean uplift? In groups such as plants, the changing conditions as elevational gradients develop have created ecotones along which species have diversified 8,17 . In butterflies, the changing range of host plants in response to gradients, and shifts to new hosts plants in parallel with these changes, have driven diversity 7,18 .
For the damselflies of Polythoridae, habitat specialization is less obvious. Chalcopteryx are found mainly in lowland and Amazonian stream habitats for example, while other groups, such as Euthore, are found at higher elevations. Within Polythore we see one distinct Amazonian clade, while the majority of the remaining species are found at higher elevations throughout the Andes. As organisms with an obligate aquatic life stage, the streams used by these damselflies do not appear to differ greatly at different elevations.
The assemblage of prey available may be different in and around these streams, but as generalist predators odonates are able to take advantage of whatever prey may be available, as long as it is in sufficient abundance.
The most outstanding group within the Polythoridae from a diversity standpoint is the Andean Clade of Polythore (Fig 1). What might be driving the increase in Andean Polythore diversity? While all these species have relatively similar habitat requirements, they have generally disjunct distributions within the Andes, such that any local stream normally hosts only a single Polythore species, sometimes two species. The wing color diversity characteristic of Polythore is highest within these groups--the central Andean Polythore have wings that include bands of black, orange and yellow patterns, while those of the Northeastern Andes have intense black and white patterning (except for P. concinna, which is orange). This color diversity also appears to be polymorphic within some of these clades 12,13 ; do this color diversity represent an adaptive radiation within Polythore? Some of the color diversity could be explained by sexual selection, with local mate choices driving diverse color patterns in different regions 19,20 . Likewise, wing color may be under other constraints, such as thermal tolerance or selection by predators 21 . Having robust phylogenetic hypotheses will allow for further exploration of the reasons behind this radiation in Polythore, where population-level analysis will be a likely next step to disentangling the complex history of these striking creatures.  Supplementary Table S6. DNA amplification, sequencing, and alignment For the 12 species new to this analysis we extracted DNA from either the legs or ¼ of the pterothorax using a DNeasy Tissue Kit (QIAGEN) from each specimen following the manufacturer's protocol. We amplified three mitochondrial and four nuclear fragments:

Taxon Sampling
Cytochrome Oxidase I (~799bp), NAD dehydrogenase (~548 bp), 16S, (~340bp), Elongation Factor (~900 bp), 28S (~340bp) and 18S (~600bp) (see Supplementary Table   S7 for a list of primers used). All gene fragments were amplified using PCR conditions as described in the associated publications for each pair of primers (Supplementary Table S7) products (15μl final volume for each primer) and the Sanger DNA sequencing. Primer contig assembly, peak chromatogram verification and the generation of per-individual consensus sequences were done using Geneious v 8 22 . All fragments were aligned using MAFFT 23 and then manually aligned in Mesquite 24 . Ribosomal genes were aligned manually with reference to secondary structure using the methods described in Kjer 25 and Kjer et al. 26 . Finally, all genes were concatenated using Mesquite for the overall analyses.

Time Divergence Analysis
A relaxed-clock molecular dating analysis on the partitioned dataset was run using BEAST v 1.8.4 27 . Specifically, we partitioned the gene fragments as follows: (i) We linked the sites and clock models for all mtDNA fragments and (ii) unlinked all nuclear ones from their clock and site models. We implemented the appropriate model selection for each partition: HKY + G4 for all mtDNA, JC for 18S and PMRT, 28S and for EF1 we set the model to GTR + G, models obtained using the model selection tool of IQTree 28 . We used lognormal relaxed clock models for all partitions, under a Yule speciation model. We generated a best ML phylogram with proportional branch lengths using the appropriate partition models using IQTree 28 as the starting tree for the analysis. Table 3 shows the calibrated nodes (4), stem fossils (10) and prior distributions selected for the analyses.
We ran four independent analyses to ensure convergence of the MCMC; convergence was checked using Tracer 1.7 29 . Finally, the independent runs for each treatment were combined using LogCombiner v 1.8.3 27 . The dated ultrametric tree was obtained using TreeAnnotator v 1.8.3 27 and visualized using Figtree v 1.4 (http://tree.bio.ed.ac.uk/software/figtree/) and the R package ggtree 30 . With the best ultrametric tree we obtain a lineage through time plot using the R package APE 31 .

Biogeographical Reconstructions
The ancestral range estimation was performed with the following software: R package Biogeobears v.0.2.1 32 and RASP (Reconstruct Ancestral State in Phylogenies) 33 . Both software packages allowed us to customize dispersal rates matrices and time stratification events, as well as infer areas among the following historical biogeography frameworks: and RASP. Afterwards, we favored the model based on best fit for the empirical information on Polythoridae and method consistency. The best selected reconstructed areas models, for each scenario were mapped over the best time calibrated phylogeny and the directionality of the dispersal and vicariance events was represented using maps for each major clade. 32,38 Time Diversification Analyses To investigate the patterns of diversification and extinction rate variation through time and across lineages, we chose RevBayes 40,41 to test the best diversification models and rate shifts, due to the high uncertainty associated with phylogenetic tree estimation. We Once the best model through time was selected we calculated the speciation (λ), extinction (μ), net diversification (λ-μ) and relative extinction (λ/μ) for the best model; all plots were generated using the RevGadgets package in R 40 . In addition, we used the Branch Specific Diversification Model implemented in RevBayes 40,46 to detect rate shifts across lineages. For all the models implemented in RevBayes 40 we assumed a uniform taxon sampling and an incomplete sampling fraction of 48/57. For the model selection, we ran a total of 5000 MCMC generations with a burnin of 1000 generations, all the parameter outputs were checked in Tracer 29 to assess the ESS, prior and posterior probability, and markov chain proper behavior. We estimated the mlnL of the models sampling the power posterior of each model for 1000 MCMC generations with a burnin of 10%. For the branch specific estimation, we assumed the heterogeneous model and we test several rate categories (N=1(constant), 4 and 10) across the lineages, for each rate we obtained the marginal probabilities and Bayes Factors as explained above. For the best selected model, we wrote a file with all the estimated parameters (λ avg, μ avg, λ-μ avg, λ/μ avg and its 95% confidence ranges) for each branch that can be mapped over our BEAST time divergence tree and was visualized using FigTree (http://tree.bio.ed.ac.uk/software/figtree/). Declarations performed field collection and manuscript drafting. RN performed analysis and drafted components of the manuscript. CS and JW contributed to manuscript preparation.

Competing interests
The author(s) declare no competing interests.

Data Availability
Data are available in GenBank; see Supplementary Table S6 for accession numbers.  Tables   Due to technical limitations, table ___ is only available