Disjunct distributions within a species are of great interest in systematics and biogeography. This separation can function as a barrier to gene flow when the distance among populations exceeds the dispersal capacity of individuals, and depending on the duration of the barrier, it may eventually lead to speciation. Here we describe patterns of geographic differentiation of two disjunct populations of Diglossa brunneiventris separated by approximately 1000 km along the Andes. Diglossa brunneiventris vuilleumieri is isolated in northern Colombia, while Diglossa brunneiventris brunneiventris has a seemingly continuous distribution across Peru, Bolivia, and Chile. We sequenced mitochondrial and nuclear DNA of the two Diglossa brunneiventris subspecies to evaluate whether they form a monophyletic clade, while including the other three species within the carbonaria complex (D. gloriosa, D. humeralis and D. carbonaria). We also constructed ecological niche models for each Diglossa brunneiventris subspecies to compare their climatic niches. We found that when using all available molecular data, the two D. brunneiventris subspecies are not sister lineages. In fact, each subspecies is more closely related to other species in the carbonaria complex. Our niche modeling analyses showed that the subspecies are occupying almost entirely different climatic niches. An additional, and not expected result was that the carbonaria complex might encompass more cryptic species than previously considered. We suggest reevaluating the taxonomic status of these brunneiventris populations, especially the northern subspecies, given its highly restricted range and potential threatened status.

Taxon sampling

We selected two localities where populations of the northern subspecies are known to occur, Páramo del Sol in the northern part of the Western Cordillera, and Serranía de las Baldías in the Central Cordillera. A total of 14 days of field work resulted in the capture of 18 individuals using mist nets. Three individuals were collected as reference specimens (two in the Central Cordillera and one in the Western Cordillera), and were deposited in the Museo de la Universidad de Antioquia in Medellín (MUUA), Colombia. Additionally, we took a small blood sample (ca. 20 µl) from 15 individuals (Serranía de las Baldías, Central Cordillera, N = 10; Páramo del Sol, Western Cordillera, N = 5). Eventually, DNA was extracted from two tissue samples and one blood sample (Table 1). For the southern subspecies we used tissue samples from three individuals collected in three different localities of Peru, all deposited in the Louisiana State University Museum of Zoology (LSUMZ; Table 1). To include all other known taxa (including subspecies) for the carbonaria complex and the three species used as outgroup (D. caerulescens, D. gloriossisima and D. albilatera), we collected blood samples in the field and obtained tissue samples from different museum collections (Table 1). Blood samples were collected using capillary tubes after brachial venipuncture with sterile needles. Blood was preserved in Queen’s lysis buffer (Seutin et al. 1991) in microcentrifuge tubes. Additional sequences of D. b, brunneiventris were available in Genbank for ND2 and CytB and were included only in analyses based on mitochondrial markers (Table 1).

DNA extraction, amplification and sequencing

We extracted DNA from pectoral muscle tissue using two different isolation kits. DNeasy tissue extraction kit (Quiagen, Valencia, CA, USA) and the UltraClean Tissue and Cells DNA isolation kit (MO BIO, Carlsbad, USA). DNA was extracted from blood using the salt extraction method. Two coding mitochondrial genes (NADH dehydrogenase subunit 2; ND2; 942 bp, and cytochrome b; CytB; 944 bp) and two non-coding nuclear introns (ornithine decarboxylase introns 6 and 7, ODC, 538 bp; and. beta Fibrinogen intron 7, βFib, 849 bp) were analyzed in this study (Prytchiko and Moore 1997, Friesen et al. 1999, Sorenson et al. 1999). PCR amplification was done in 30 µl reaction volumes including 2 mM MgCl2, 2 µM dNTPs, 0.5 µM primers, 1X Taq buffer, 1.8 units of Taq polymerase (Thermo Scientific) and ca. 30 ng DNA. Amplification reactions of ND2 also included 0.4% bovine serum albumin (BSA). Thermal cycling conditions included an initial denaturation step at 96 °C for 5 min followed by 35 cycles of denaturation at 96 °C for 30 s, annealing at marker specific temperature (60 °C for βFib, CytB, and ND2, 57 °C for ODC) for 30 s and extension at 72 °C for marker specific times (ODC 30s, βFib 60s, ND2 and CytB 70s), and ended with a final extension step at 72 °C for 10 min. Reactions were purified and sequenced by Macrogen Inc. (Korea). Sequences were manually edited and aligned with Clustal W in BioEdit v 7.0.5.3 (Hall 1999). Sequences of CytB and ND2 were translated in DnaSP to make sure they were the product of amplification from the mitochondrial genome. We confirmed that all sequences obtained from protein-coding genes did not have stop codons and therefore were not likely to be the product of nuclear copies of mitochondrial DNA amplification (e.g., numts). All sequences are new and were deposited in GenBank (accession numbers MW456565-MW456637).

Phylogenetic analyses

Sequences were aligned using Sequence matrix 1.8 and additionally by eye. We partitioned the sequence data by locus and coding position and used ModelFinder (Kalyaanamoorthy et al. 2017), following the Bayesian Information Criterion, to determine the evolutionary models for each partition. To estimate tree topology and support values for clades, we used Bayesian (MrBayes) and maximum likelihood (IQtree) approaches using the mitochondrial (ND2 and CytB), the nuclear (BFib and ODC), and the complete partitioned datasets. In MrBayes, we ran partitioned analyses using three separate chains for 30 million generations, sampling every 10000, and discarding the first 1000 sampled trees as burnin. In IQtree (Nguyen et al 2015), we also ran partitioned analyses and generated support values using the ultrafast bootstrap algorithm with 10,000 bootstraps (Minh et al. 2013). To further estimate divergence times among sequences we used BEAST v 1.8.3 with the complete dataset (Drummond and Rambaut 2007, Drummond et al. 2012). We allowed separate substitution and clock models for each partition (including codon position) but constrained analyses to a single-tree model for all sequence data. The analysis incorporated a strict clock and a Yule speciation process with a random starting tree. We chose the strict clock option given demonstrated support for a constant evolutionary rate in mitochondrial genes for a time interval of 12 my (Weir and Schluter 2008), and the Yule speciation process (birth process with constant rate) given that our clade of interest is relatively young (< 1 my) and thus we do not expect extinctions to be frequent. Substitution rates for each locus followed the guidelines in other phylogenetic analyses of birds (Lerner et al. 2011, McGuire et al. 2014). Sequences from the carbonaria complex were used as ingroup, and sequences from D. albilatera, D. gloriosissima and D. caerulescens as outgroups. We used the estimated divergence times given in Mauck and Burns (2009) for the most recent common ancestor of the whole group including D. caerulescens (15.1 mya, 95% Highest Posterior Density interval [HPDI]: 10.4-20.5 mya) and for the carbonaria complex (0.5 mya, 95% HPDI: 0.4-0.6 mya) but we did not enforce monophyly for the carbonaria complex since the placement of samples from the northern subspecies of D. brunneiventris was uncertain. The analysis was run for 50,000,000 generations, sampling every 10,000 generations (extra parameter details are provided in Supplemental Material Appendix A). Output log files were examined in Tracer v.1.8.3 (Rambaut et al. 2017). After removing 20% of samples as burn-in, independent runs were combined and a maximum clade-credibility (MCC) tree was constructed using TreeAnnotator v.1.8.3. FigTree v 1.4.3 was used to edit the tree. In addition, we constructed minimum spanning haplotype networks (Bandelt et al. 1999) for CytB and ND2 in PopART (Leigh and Bryant 2015) to visualize evolutionary relationships between the haplotypes, and estimated genetic distances among sample sequences using MEGA X (Kumar et al. 2018).

Environmental niche modeling

We compiled a database of occurrence records from both populations of Diglossa brunneiventris from the Global Biodiversity Information Facility (GBIF) as well as from our own records (Fig. 1). Each record was evaluated for inconsistencies in the locality description, elevation, taxonomy and georeference. Suspicious records (i.e. locality not corresponding to coordinates, inconsistent elevation) were discarded. After all occurrence records were checked for quality, we used them to estimate and compare the climatic niches of each subspecies following the methodology proposed by Broennimann et al. (2012) in which occurrence density surfaces are generated via kernel density estimation based on the environments available and the frequency of occurrences in each accessible environment. By environment, we specifically refer to the 19 bioclimatic variables available from the WorldClim database (Hijmans et al. 2005, v1.4; Supplemental Material Table S1). These variables represent climate variation (mean, seasonality and extremes) for temperature and precipitation at ~ 1 km² spatial resolution and for an interval of 50 years (1950-2000). Many of these climatic variables are highly correlated and the niche comparison is based on a principal component analysis of them. The methodology proposed by Broennimann et al. (2012) estimates climatic preferences by comparing the available or accessible climate (determined by all climates within the accessible area) in relation to the occupied or used climate (determined by all climates occupied by the species). The accessible area, understood as all areas in geographic space that the species is able to access, is referred to as M in ecological niche modeling theory (Soberón and Nakamura 2009, Barve et al. 2011). Since the method is based on occurrence frequencies in environmental space and to prevent correlation among detections, we chose to use only a single record per 1 km² pixel. We determined the accessible environments for each subspecies as those within the ecoregions (Dinerstein et al. 2017) where at least one occurrence of the population was identified (Fig. 1). Finally, to formally evaluate how similar the niches of the two populations of D. brunneiventris were, we quantified their overlap using a similarity metric developed by Schoener (1968), but applied in environmental space as proposed by Broennimann et al. (2012). This metric essentially quantifies the amount of climatic space used commonly by both populations, acknowledging the climates that are accessible for each one. We used two null model approaches developed originally by Warren et al. (2008) to evaluate if similarity is higher than expected under a proposed null model. These two are known as the niche equivalency and the niche similarity tests. The null model of the niche equivalency test randomizes the identity of the population that each locality belongs to and calculates the overlap between  populations (northern and southern subspecies) for a given number of iterations, while the null model of the niche similarity test randomizes the position of the occurrences within the accessible area for one population while holding the other fixed and calculates the overlap between populations (northern and southern subspecies) for a given number of iterations. In either case, the observed overlap is compared to the distribution of overlap values obtained by the given null model. The expectation under the null hypothesis for the niche equivalency (identity) test is that both populations share identical climatic preferences. If rejected, then the subspecies exhibit differences in their climatic preferences. For the niche similarity test, the null hypothesis is that subspecies differences in climatic preferences are due to differences in climate availability. Thus, if the null hypothesis is rejected, climatic preferences can be more or less different than expected under the climates available in each area.