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 P . gymnocephalusacross 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).