The colonization of the Puna and Atacama Biogeographic Province by sister clades of Psectrascelis (Coleoptera: Tenebrionidae): Synchronous expansion without spatial overlap

We investigated the biogeographic pattern of the species‐rich genus Psectrascelis (Coleoptera: Tenebrionidae) in the Central Andes and the Chilean Atacama Desert as an example for insect evolution in such a geologically and climatically highly dynamic arid region. The main aim was to test two alternative hypotheses about the main drivers of diversification of biotas from deserts in the region: Andean uplift versus Pleistocene fragmentation/dispersal.


| INTRODUC TI ON
Understanding and characterizing thresholds for biological colonization and speciation and their potential causes, such as species diversification in response to climate and geological processes (e.g. Gillespie & Roderick, 2014) or interplay between geographical barriers and species dispersal (e.g. Burrows et al., 2014) is a major topic in biogeography. We are particularly interested in biological evolution and speciation at the transition from aridity to hyperaridity. In deserts, simple food webs and limited biological interactions offer ideal conditions for detailed research into key factors in the past that may have influenced the current distribution and diversification of organisms. Such analyses can then be used to model future trends in regions with increasing aridity, including possible time scales for adaptation and speciation processes.
Our study area is part of the South American Transition Zone (SATZ; Morrone, 2014Morrone, , 2015 and includes the Chilean Atacama Desert, adjacent regions of the Pacific coast of Chile and the Central Andes. The area of the Central Andes considered here, including the Western and Eastern Cordilleras, the Altiplano and the Puna south of the Altiplano, is known as the Puna Biogeographic Province (Puna BP; Morrone, 2014). It is well established that the uplift of the Andes to its current altitude during the Miocene (Garzione et al., 2017) caused significant climatic changes during the Pliocene with strong impacts on the biodiversity throughout South America (Graham, 2009). Since then, the Andes have functioned as an effective eastwest migration barrier for fauna and flora and thus a major driver for the diversification of many plant and animal taxa that have been studied from this region so far (e.g. Antonelli & Sanmartin, 2011;Bohnert et al., 2019;Ceccarelli et al., 2016;Pocco et al., 2013), including insects (e.g. Cigliano, 1989;Roig-Juñent, 2002). The interplay of tectonic activity and changing climate led to the proliferation of arid environments in the southern cone known as Arid Diagonal of South America (ADSA; Abraham et al., 2000;Hinojosa & Villagran, 1997), providing favourable conditions for the dispersal of organisms adapted to aridity (e.g. Ceccarelli et al., 2016;Zúñiga-Reinoso & Predel, 2019). In the eastern part of the ADSA, the Andes uplift was accompanied by hyperaridity in coastal Pacific regions (Jordan et al., 2014), caused in part by the rain shadow effect of the mountain ranges. Species diversification was also affected by global climate changes that overlapped temporally with the Andes uplift. These climate changes have repeatedly interconnected or fragmented suitable habitats. In addition, the emergence of arid habitats favoured adaptive radiations as a result of the availability of new ecological niches (Hughes & Eastwood, 2006;Antonelli & Sanmartín, 2011).
Taking all this into account, the complex geoclimatic structure in our focus area seems to provide a very suitable starting point to test effects of tectonic and/or climate changes on the diversification of organisms under conditions of increasing aridity. For the present study, we used wingless beetles of the family Tenebrionidae (Insecta: Coleoptera), which generally have low dispersal abilities.
The Tenebrionidae belong to the most species-rich animal families with about 20,000 species and, in contrast to most animal groups, show an increasing species diversity with increasing aridity (Koch, 1962). In deserts, which in the past were often colonized in a convergent scenario by different lineages of Tenebrionidae, they are among the ecologically most important consumers. It is still unclear which genetic predisposition of Tenebrionidae favours worldwide desert colonization, but it is known that such colonization is usually followed by subsequent radiations in new and often still unoccupied ecological niches, which is then accompanied by numerous adaptations that affect morphology, physiology and behaviour (Cloudsley-Thompson, 2001). It was postulated for Tenebrionidae in our focus area that the Central Andes uplift acted as a vicariant event at genus and species level (Dominguez et al., 2016;Flores & Pizarro-Araya, 2006). With Psectrascelis, we have selected one of the most speciesrich tenebrionid genera of South America (Flores, 1997;Peña, 1985), which is distributed across the ADSA. Of the 83 species described, 42 species are known from central and northern Chile, making this region particularly biodiverse for Psectrascelis (Flores, 1997;Larrea-Meza et al., 2020;Peña, 1985). So far, 24 species have been documented from the Puna BP, of which only five inhabit the Chilean section (Ferrú & Elgueta, 2011;Peña, 1985Peña, , 1994. Unlike several other tenebrionid taxa, which are already adapted to hyperaridity and colonized the hyperarid core of the Atacama (e.g. Physogaster, Entomochilus, Psammetichus, Cordibates), species of the genus Psectrascelis largely avoid hyperarid areas and their current distribution pattern (Figure 1), therefore nicely illustrate the boundaries of hyperaridity in the Atacama. Psectrascelis' wide distribution as well as its less developed adaptation to hyperarid conditions makes it particularly suitable for analysing the evolution of insects in a geologically and climatically very dynamic arid environment. However, molecular data are largely missing, and the divergence times, the colonization histories and hence the key factors that led to their current distribution and species diversity still need to be verified.
The main objective of our study is, therefore, the examination of the biogeographic history of Psectrascelis to determine the driving forces behind the distribution and diversification of taxa adapted to high altitudes and aridity. We have used tenebrionid species as they form the largest radiation of insects in hyperarid deserts worldwide, and Psectrascelis is considered as taxon with progressive adaptation to hyperaridity. In this context, we tested whether the distribution and molecular diversity of Psectrascelis can be related to historical geoclimatic events. In particular, we tested two alternative hypotheses about the main drivers of diversification of biotas from deserts in the region: Andean uplift versus Pleistocene fragmentation and/ or dispersal.  (1985,1994) and Vidal and Guerrero (2007).

| DNA extraction, amplification and sequencing
DNA was extracted from thoracic muscles using the EZNA® Insect DNA Kit (Omega Bio-tek, Inc., Norcross, U.S.A). First, the mitochondrial gene cytochrome oxidase I (COI) was amplified for 87 specimens. To facilitate the detection of different haplotypes in more conserved genes, we selected the genetically more distant specimens of each Psectrascelis lineage (see Zúñiga-Reinoso & Méndez, 2018) and subsequently amplified the ribosomal RNA 16S and the nuclear gene wingless (Wg) from these samples. The primers used to amplify each gene are listed in Table S1. The thermal polymerase chain reaction (PCR) steps were 94°C for 2 min, followed by 36 cycles at 94°C for 30 s, 56°C for 45 s and 72°C for 1 min, with a final extension at 72°C for 2 min. The PCR products were purified using peqGOLD Cycle-Pure kit (Peqlab Biotechnologie GmbH, Erlangen, Germany), and both forward and reverse products were sequenced at Eurofins Genomics GmbH (Germany) using the Sanger sequencing method. The DNA sequence of each sample was reviewed, and all orthologous sequences were aligned using the Clustal W algorithm implemented in BIOEDIT version 7.0.5.3 (Hall, 1999) and then manually checked for inconsistencies. In addition, Xia's test (Xia et al., 2003) implemented in DAMBE, v5.1.5 (Xia & Xie, 2001) was used to evaluate the saturation of each gene matrix.

| Phylogenetic reconstruction
Phylogenetic reconstructions were performed with the concatenated COI, 16S and Wg genes using Bayesian inference (BI) algorithms. The best model for sequence evolution was selected with the Akaike Information Criterion (AIC) in the program jModelTest 0.1.1 (Posada, 2008). The best models for the sequences were TIM2+I+G for COI and Wg and TIM2+G for 16S; these models were used for each partition in all analyses. The BI was performed with the program Mr. Bayes 3. 2. 6 (Huelsenbeck & Ronquist, 2001) implemented in the server CIPRES Science Gateway 3.3 (Miller et al., 2010). We conducted four independent runs, with a setting of four chains; starting with a random tree, running for 30 million generations and sampling every 1000 trees. The initial 25% of the resulting trees was discarded as burn-in. Once convergence of the four independent runs was confirmed by the average standard deviation of split frequencies and the potential scale reduction factor, results from the runs were combined to obtain a total of 30,004 trees.
Finally, a consensus tree was constructed by a 50% majority rule; the node supports were evaluated using posterior probabilities.

| Divergence times and reconstruction of biogeographic history
For an estimation of the divergence times, we pruned the branches of the BI tree with the 24 specimens whose COI, 16S and Wg genes have been determined. Estimation of divergence times was performed with the software BEAST 1.8.2 (Drummond & Rambaut, 2007). The substitution models selected for the jModelTest in BI are not implemented in BEAST; these models were replaced by the complex models GTR+I+G for COI and Wg and GTR+G for 16S. In addition, lognormal relaxed clock and Yule speciation models were selected, and the analysis was then run for 30 million generations, sampling every 3000 generations. For the calibration of divergence times, the mutation rate for COI was set at 3.54%, according to the estimate for Tenebrionidae by Papadopoulou et al., (2010). Results were checked with TRACER 1.6 (Drummond & Rambaut, 2007) to determine the convergence of the chains and the effective sample size. The sampled trees were combined with TreeAnnotator 1.8.2 (Drummond & Rambaut, 2007) and used to plot the mean of the lineages through time and the respective 95% high posterior density (HPD) in TRACER 1.6. The tree was displayed and edited in FIGTREE 1.4.2 (Rambaut, 2009).
The reconstruction of ancient distribution ranges and speciation events was performed with the software package 'BioGeoBEARS' (Matzke, 2013), which allows a probabilistic deduction of the biogeographic history of a taxon and has the advantage of comparing different biogeographical models. This approach considers both anagenetic events (i.e. dispersion, extinction and shift of range) and cladogenetic events (vicariance, founder events; Matzke, 2014).
BioGeoBEARS was performed with the models DEC, DIVALIKE and BAYAREALIKE. Additionally, we performed the same models with the extra parameter 'j' (i.e. dispersion by jump; Matzke, 2014). All models implemented in the BioGeoBears were run independently, and the best model was selected using AIC. For our analyses, we used the ultrametric tree of BEAST and a presence/absence matrix with the biogeographic region where the species are currently distributed. Primarily, the biogeographic provinces proposed by Morrone (2014Morrone ( , 2015 were used. We have modified the data for the Puna BP according to the classification proposed for Chile based on the distribution of Tenebrionidae (Peña, 1966) and other terrestrial animals (Artigas, 1975).

| RE SULTS
Sampling in the Chilean Atacama Desert and Chilean Puna BP ( Figure 1) yielded all species known from this region (P. intricaticollis, P. laevigata, P. escobari, P. confinis, P. izquierdoi) and added several locations not previously listed for Psectrascelis. The following morphologically determined species were collected from neighbouring regions: (1) P. laevigata (region Lake Titicaca, Peru). This species is one of the few Psectrascelis known to occur north of the focus area.
(2) P. rotundata (Villa Mar, Bolivian Altiplano). This species belongs to several morphologically similar species in the region; most are known from the Argentinian Puna (Kulzer, 1954;Peña, 1994). (3) P. elongata, P. specularis, P. sublaevicollis, and P. pilipes. The latter species belong morphologically to Psectrascelis sensu stricto that become species-rich in the Coquimban BP south of the focus area and include coastal species. In total, 87 Psectrascelis individuals from 37 locations were included in our study (Table S2).

| Phylogenetic reconstruction
The combined matrix of COI, 16S and Wg genes has a length of 2460 bp and contains 572 polymorphic sites and 552 parsimony informative characters (for further details, see Table S3 and S4). No genetic saturation was observed for these genes.

| Divergence times and reconstruction of the biogeographic history
Different selection models for biogeographic estimations were tested (Table S5), and the best fit were obtained with DIVALIKE+J F I G U R E 2 Phylogenetic analysis of Psectrascelis using a concatenated data set of (COI, 16S and Wg) genes. The numbers at the nodes show the posterior probabilities. The coloured boxes of the different lineages correspond to the colour-marked areas in the map, which show the hypothetical distribution of these lineages. For this purpose, data from Peña (1985) were cross-checked with the material investigated in our study. The dotted lines are the average annual rainfall isohyets for 100 (red), 250 (green) and 500 mm (blue) (Fick & Hijmans, 2017) [Colour figure can be viewed at wileyonlinelibrary.com] and DEC+J; both models yielded identical reconstructions. Clades

| DISCUSS ION
The present study provides the first detailed molecular analysis of the diversification and colonization history of insects across the SATZ, biogeographically covering the Puna BP and the Chilean    (Zachos et al., 2001). During the Mid-Pleistocene Transition, the previous dominant periodicity of climate cycles changed from a 41-kyr to a 100-kyr cycle, causing high-amplitude climate oscillations (e.g. Clark et al., 2006;Tziperman & Gildor, 2003;Zachos et al., 2001

| Allopatry and adaptation
Our results suggest that speciation in Psectrascelis was clearly linked to changes in global climate and climate variability in the last 3.5 Ma. Expansion and shrinking of suitable habitats might have regularly caused fragmentation of populations; primarily due to a failure to adapt to the changing environmental conditions.
We hypothesize that such niche conservatism, which resulted in non-adaptive allopatric speciation (Predel et al., 2012;Rundall & Price, 2009;Wiens, 2004), is particularly typical for Psectrascelis clade A. The hypothetical origin of clade A is thought to be on the eastern slopes of the Central Andes in north-western Argentina; this is the proposed eastern retreat of the mrca of Psectrascelis during the major climate cooling after the Mid-Pliocene warm period, demonstrated by glaciation in the Central Andes 3.5 Ma (Haselton et al., 2002;Roberts et al., 2018). Today, clade A species more or less continuously occupy the entire high-altitude region of the Central Andes, except for the dry western slopes of the Western Cordillera, inhabited by the B clade (see below), and the Altiplano. All species of the A clade were collected in areas receiving precipitation of about 100 mm/year or more (Fick & Hijmans, 2017), including the northernmost species of this clade from an area west of Lake Titicaca that receives precipitation of more than 500 mm/year (Figure 2). The assumed non-adaptive radiations in Psectrascelis do not exclude subsequent adaptations or phenotypic changes, but speciations likely happened prior to these changes.

ACK N OWLED G EM ENTS
We thank the following colleagues for the support of the study:

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets generated during the current study are available in the NCBI repository (see table S2 for accession numbers), and the complete genes data set is available in the CRC1211 project repository under the https://doi.org/10.5880/CRC12 11DB.37.

B I OS K E TCH
Alvaro Zúñiga-Reinoso is mainly interested in entomology and evolution, using phylogenetic approaches to understand the systematics, diversification processes and biogeographic scenarios that drive the insect diversity.
Benedikt Ritter works in the field of geochronology, sedimentology and paleoclimate reconstructions with special focus on arid to hyperarid climates/environments, particularly in the Atacama and Namib Desert.