Introduction
The glacial cycles of the Pleistocene recurrently brought formerly isolated populations into contact with one another during cooler periods that increased low-elevation distributions (Hewitt 1996, Hewitt 2000, Avise 2009). Following the last glaciation period, cooler-adapted species either shifted ranges to higher latitude or else became relicts in mountaintop refugia (”sky islands”; McCormack, Huang, Knowles, Gillespie, & Clague, 2009); sky islands are analogous to ocean archipelagos, with individual mountain populations isolated due to low-elevation habitat that acts as a barrier to dispersal. Predicted increases in temperatures due to climate change will likely result in increased isolation of these high-altitude species as lowland habitat expands and alpine habitat shifts upslope (Chamberlain et al., 2012; Kupfer & Cairns, 1996; Scridel et al., 2018; Sekercioglu, Schneider, Fay, & Loarie, 2008). Depending on the distance between populations and their ability to disperse, this isolation can lead to a feedback cycle of increased inbreeding depression among populations descending into an ”extinction vortex” (Frankham, 2015).
The inhibited gene flow and localized genetic drift resulting from decreased dispersal in sky islands means alpine species often develop genetic differentiation among isolated populations (McCormack et al., 2009). This isolation can be seen in geographically-based genetic structuring (e.g. Bech, Boissier, Drovetski, & Novoa, 2009; Jackson, Gergel, & Martin, 2015; Lonsinger, Schweizer, Pollinger, Wayne, & Roemer, 2015). By examining genetic structuring in isolated populations of species, areas of low connectivity may be identified (Allendorf & Luikart, 2009). Especially for species at risk, understanding the genetic diversity underlying the spatial structure of populations is important for establishing the appropriate scale and subunits for conservation management (Laikre, 2010). While it has been suggested that species in isolated populations would evolve more efficient dispersal ability (Travis et al., 2013), effective dispersal distance (i.e. that which results in the disperser reproducing and adding to the gene pool) is constrained by population growth rate and distance between habitat patches (Baguette & Schtickzelle, 2006; Van Dyck & Baguette, 2005). However, even if populations are not separated by geographic barriers, they may have low genetic diversity due to strong selective pressure if they occupy highly specialized habitat niches (Orsini, Vanoverbeke, Swillen, Mergeay, & De Meester, 2013).
As would be expected, genetic structuring in sky island endemics has been seen across nearly all taxa, including plants (e.g. Gizaw et al., 2013; Lexer et al., 2013), arthropods (e.g. Knowles & Richards, 2005; Masta, 2000), amphibians (e.g. Vörös et al., 2017), birds (e.g. Bech et al., 2009; Jackson et al., 2015), and mammals (e.g. Browne & Ferree, 2007; Lonsinger et al., 2015). The sky islands of South Africa’s Cape Fold Belt have been overlooked for such studies among birds, despite the possibility of elevation-based species diversification (Verboom, Bergh, Haiden, Hoffmann, & Britton, 2015), and presence of high levels of endemism (proportional to geographic size) within the Fynbos biome across taxa (Cowling, 1992; Picker & Samways, 1996; Sharratt, Picker, & Samways, 2000; Wishart & Day, 2002). Here, we focused on the Cape Rockjumper Chaetops frenatus , an avian sky endemic restricted to the mountain fynbos of southern African sky islands. The contemporary geographical distribution of mountain fynbos results from recent historical shifts in precipitation regimes; it expanded ~96–37 KYA (”KYA”: thousand years ago) and then contracted resulting in the current isolated patches of alpine fynbos across the Cape Fold Belt (Quick et al., 2016). Population declines correlated with warming habitat (Milne, Cunningham, Lee, & Smit, 2015), as well as a small population size, led to C. frenatus being placed as Near Threatened on the IUCN Red List of Endangered Species in 2017 (IUCN, 2017). The possible continued contraction of mountain fynbos habitat due to the expansion of warmer lowland habitat from climate change may be one reason for the decline and isolation ofC. frenatus populations in warmer areas of their habitat.
In this study, we aimed to determine genetic structuring among sky island populations of C. frenatus using a combination of mitochondrial DNA (mtDNA) and nuclear DNA (nDNA) to test ifC. frenatus are experiencing geographically-based genetic isolation. These three markers were chosen as they remain relevant in studies on avian population structure (e.g. Liu et al., 2016; Dantas et al., 2019; Stervander, Ryan, Melo, & Hansson, 2019), as well as general avian taxonomy (Ericson, Klopfstein, Irestedt, Nguyen, & Nylander, 2014). We hypothesized that unsuitable lowland habitat between the geographically distinct mountain ranges would act as a barrier to effective dispersal between populations, resulting in genetic structuring aligned with the topography of the Cape Fold Belt. If we find geographic genetic structuring, suggesting some populations ofC. frenatus are genetically isolated from one another, this may indicate key areas in need of conservation. Alternatively, if selective pressure is at work due to the highly specialized habitat occupied byC. frenatus , we may find low overall genetic diversity and no genetic structuring (Orsini et al., 2013). As both sister-species within family Chaetopidae represent small relict groups at the far end of early radiation of African passerines of high taxonomic interest (Fjeldså & Bowie, 2008), we chose the Drakensberg Rockjumper C. aurantius as our outgroup. We also chose Picathartes gymnocephalus as an additional outgroup because the exact taxonomy and phylogenetic placement between the sister families Picathartidae and Chaetopidae, and within Chaetopidae itself, remains unresolved within the oscines (Ericson, Klopfstein, Irestedt, Nguyen, & Nylander, 2014).
Methods
Study area and sample collection
The study area encompassed much of the Western and Eastern Cape provinces of South Africa, spanning the species distribution range for the C. frenatus . During 2016 and 2017, we collected 73 blood samples of C. frenatus from 13 localities representing eight mountain ranges of the Cape Fold Belt (see Figure 1 for localities as well as samples per mountain range); these eight ranges represent areas with >10 % recorded occurrence rate for C. frenatuson the South African Bird Atlas Project (”SABAP”) 2 (2007 - present; http://sabap2.adu.org.za), and were sampled as extensively as was possible. For two additional localities we attempted to sample C. frenatus based on records from SABAP 1 (1987–1991) as well as from earlier records of SABAP 2 (<2015), but were unable to locate birds at these localities (e.g. Tsitsikamma and Lady Slipper; Figure 1).
Although generally contiguous, many of the mountain ranges in the Cape Fold Belt are interspersed by >20 km of lowland habitat (Lee & Barnard, 2016). This lowland habitat is generally warmer than habitat occupied by C. frenatus , and is often also either drier (e.g. Karoo) or densely vegetated (e.g. Albany thicket), thus representing a potential dispersal barrier for the C. frenatus(McCormack et al., 2009). Sample localities ranged from 32° to 34° S latitude and 18° to 25° E longitude, spanning from ~100 metres above sea level (”masl”) near the southwest tip of South Africa to ~1,800 masl in the northern Cederberg mountains (Figure 1). Despite the low elevation occurrence of C. frenatusnear the southwestern tip of South Africa, this area does in fact consist of the short, Ericaceae- dominated mountain fynbos (Cowling & Heijnis, 2001). We additionally collected 15 blood samples of C. aurantius from two localities at 30° S latitude and 27° to 28° E longitude in the Lesotho Highlands, spanning from ~2,600 masl to ~2,800 masl (Figure 1). Localities for C. aurantius were also selected based on occurrence data from SABAP 2.
We caught wild-living individual C. frenatus and C. aurantius using snap traps baited with tenebrionoid spp. Upon capture, we recorded a GPS point, and fitted birds with unique numbered aluminium rings (SAFRING). Blood samples (~50 μL) were collected by brachial vein puncture, whereby a small pin-prick was made with a 26-gauge needle and the resulting blood droplet collected into a capillary tube and preserved in ~90 % alcohol. The entire sampling process took <5 min of handling time and occurred adjacent to, or within, birds’ territories for release immediately after sampling. Blood samples were stored at 4 °C until further use for DNA extraction.
DNA extraction and primer selection
Genomic DNA was extracted from blood samples in a dedicated lab following standard salt extraction techniques (Bruford, Hanotte, Brookfield, & Burke, 1992) using Proteinase K (QIAGEN, Germany; 10 mg mL−1), modified by use of lysis and elution buffers (buffers ATL and AE respectively; QIAGEN, Germany).
Primers were prepared by Integrated DNA Technologies, USA. Two primer sets were selected from Zuccon and Ericson (2012) for mtDNA markers: “ND2” which encodes the NADH dehydrogenase subunit 2 and “cytb” that encodes the cytochrome b. For the selected nDNA marker, the recombination activating gene 1 —- “RAG1”, we designed a ~550 bp primer online (www.ncbi.nlm.nih.gov/tools/primer-blast/) using sequences of C. frenatus and our additional outgroup P. gymnocephalus(representing the sister family Picathartidae) available from GenBank (www.ncbi.nlm.nih.gov/genbank/; Date: 04-03-2019; Table S1).
DNA amplification and sequencing
For amplification of mtDNA, PCR reaction occurred in a final volume of 27 µL containing approximately 100 ng of genomic DNA, 0.4 µM of each primer, and 12.5 µL of a mixture dNTPs, MgCl2, Taq polymerase and stabilizers (cytb: iTaq Universal SYBR Green Supermix® by Bio-Rad Laboratories Inc.; ND2: TopTaq PCR MasterMix QIAGEN). For amplification of the nDNA the PCR final volume was 33 µL and contained approximately 200 ng of genomic DNA, 0.6 µM of each primer, and 12.5 µL of a mixture of dNTPs, MgCl2, Taq polymerase and stabilizers (iTaq Universal SYBR Green Supermix® by Bio-Rad Laboratories Inc.).
Mitochondrial DNA PCR protocol was as follows: initial denaturing step of 94 °C for 5 min, 40 cycles of denaturing at 94 °C (30 s), annealing at 58 °C (ND2) or 60 °C (cytb; 30 s), and extension at 72°C (40 s), final extension at 72 °C (5 min). Nuclear DNA PCR protocol was as follows: initial denaturing at 94 °C (5 min), 40 cycles of denaturing at 94 °C (30 s), annealing at 51.5 °C (40 s), and extension at 72 °C (40 s), final extension at 72 °C (5 min). Sequence data for P. gymnocephalus were retrieved from GenBank for ND2 (Fuchs, Fjeldså, Bowie, Voelker, & Pasquet, 2006), cytb (Ericson & Johansson, 2003), and RAG1 (Barker, Barrowclough, & Groth, 2002); see Table S1 for accession numbers). Sequencing of the three loci was outsourced to Macrogen Europe (Amsterdam, the Netherlands); sequences were only produced in the forward direction. We visually checked each sequence individually for correct base-calling of individual nucleotides.
Phylogenetic analysis
Sequences were aligned in MEGA (v 7.0.62; Kumar, Stecher, & Tamura, 2016) with the ClustalW method. We tested for saturation in each codon of each marker using DAMBE v. 7 (Xia, 2018), and found no evidence of saturation in any of the codon positions. We estimated uncorrected p-distances both among groups/species (C. frenatus , C. aurantius , and P . gymnocephalus ) and within groups/species (C. frenatus and C. aurantius ) in MEGA to investigate the degree of divergence and examine overall genetic variation. Population structure was examined using minimum-spanning networks (Bandelt, Forster, & Röhl, 1999) for both cytb and ND2 in PopART (Leigh & Bryant, 2015; http://popart.otago.ac.nz).
Genetic differentiation between populations was examined using an Analysis of Molecular Variance (AMOVA; Excoffier, Smouse, & Quattro, 1992) in ARLEQUIN v. 3.5 (Excoffier & Lischer, 2010). Chaetops frenatus and C. aurantius sequences were split into eight and two groups respectively based on localities where samples were obtained (see Figure 1). Fixation indices were used to estimate the proportion of genetic variability found within populations (F ST; i.e. within each mountain range), among populations within groups (F SC; i.e. among mountain ranges) and between groups (F CT; i.e. between C. frenatus and C. aurantius ). Due to the lack of nucleotide differences between individuals and populations of C. frenatus and C. aurantius at the RAG 1 locus, we did not perform an AMOVA on RAG1.
Sequence data were also used to detect signal of population expansion, within each species (C. frenatus and C. aurantius ) as estimated with Tajima’s D (ARLEQUIN v3.5). When D = 0 indicates no evidence of selection, D > 1 indicates population contraction, and D < 0 indicates a post-bottleneck expansion (Tajima, 1989).
We first used jModelTest v. 2.1.10 (Darriba, Taboada, Doallo, & Posada, 2012) for each marker to determine the evolutionary model that best fit the dataset using Bayesian Information Criterion adjusted for low sample size. Results of jModelTest suggested best fits were two rate categories for each locus, with cytb using uniformly distributed proportions of invariable sites (HKY + I; Hasegawa, Kishino, & Yano, 1985), ND2 using four-category gamma distributed invariable sites (HKY + Γ + I; Hasegawa et al., 1985), and RAG1 using equal rates (K80; Kimura, 1980). Phylogenetic inference trees were created online in CIPRES Science Gateway v. 3.3 (Miller, Pfeiffer, & Schwartz, 2010) with data partitions based on each genetic marker (cytb, ND2, and RAG1), and priors set according to results of jModelTest. As we did not have sequences for the same individuals of Pgymnocephalusacross all three markers, we used available sequences for RAG1 in our consensus trees. For Bayesian inference of phylogeny we used the program MrBayes v. 3.2.6 on XSEDE with BEAGLE (Miller et al., 2010) to obtain Markov Chain Monte Carlo (MCMC) approximations of posterior trees; MCMC was run for 20,000,000 generations, sampling trees every 1,000 generations, with the first 2,500,000 generations (2,500 trees) discarded as burn-in. For maximum likelihood (ML) inference, we used Garli v. 2.0.1 on XSEDE (Bazinet, Zwickl, & Cummings, 2014) to obtain approximations of posterior trees based on 100 bootstraps; we ran this analysis twice to ensure the independent ML searches produced the same tree topology. We determined best fit phylogeny of Bayesian vs. ML trees based on greater number of node posterior probabilities of >95 % or >70 % (Bayesian and ML respectively; see Bates et al., 2013).