Population diversity and relatedness in Sugarbirds (Promeropidae: Promerops spp.)

Sample collection and DNA extraction

All samples were collected during prior ecological studies (Calf, Downs & Cherry, 2001; Calf, Downs & Cherry, 2003). Permits were obtained from the relevant provincial conservation agency (Free State Department of Environmental Affairs and Tourism and CapeNature Permits respectively to Kathleen Calf and Gordon Scholtz, 1998–1999) with mist netting of birds licensed (KMCT) under the South African national bird ringing scheme (http://safring.adu.org.za, permit number HK/P1/02095/001). Blood samples were obtained by brachial vein puncture and capillary capture (approximately 20 µl per bird) from a population of P. gurneyi (n = 48) in Golden Gate Highlands National Park, Free State, in 1998, and from one individual in Sani Pass, KwaZulu-Natal, in 2004. Samples were similarly obtained from a population of P. cafer (n = 63) in Helderberg Nature Reserve, Western Cape, in 1999, where sampling was targeted around nests. Samples were stored in 1mL of a modified PBS solution (phosphate buffered saline blood storage buffer: 3 mM KCl, 8 mM Na₂HPO₄, 2 mM KH₂PO₄, 0.14 M NaCl, 6 mM EDTA, 0.2% NaN₃ (w/v)) which was kept on ice after sampling and stored long-term at −20 °C. This buffer includes EDTA to prevent enzymatic digestion of nucleic acids, while the Sodium Azide, NaN₃, is an antimicrobial.

Samples from other species, used to determine the taxonomic range of ND2 primer amplification, were available from previous studies in our laboratory. These included mountain pipit, Anthus hoeschi, blue crane, Anthropoides paradiseus, Angolan cave chat, Cossypha ansorgei, Cape parrot, Poicephalus robustus, Namaqua sandgrouse, Pterocles namaqua, and eastern clapper lark, Mirafra fasciolata. Tissue samples of chicken, Gallus gallus domesticus and Ostrich, Struthio camelus, were obtained from supermarket meat, and marsh owl, Asio capensis, from a roadkilled individual sampled from a carcass found beside the R50 road at Delmas (Long, Lat: 028.681, −26.141).

We extracted genomic DNA by salting out. Digests were prepared by adding 1 mg blood to 500 µL DNA lysis buffer (100 mM NaCl, 50 mM Tris.HCl, 100 mM EDTA, 1% SDS w/v) supplemented with 0.2 mg proteinase-K, before incubation at 57 °C overnight. Following digestion, a volume of 20 µL RNase A (25 mg/mL) was added, before incubation at 37 °C for 60 min. Proteins and other cellular contents were precipitated by addition of 180 µL 5 M NaCl (final concentration 1.3 M) followed by agitation and centrifugation. DNA was then extracted by combining the supernatant with ice-cold isopropyl alcohol (1:1), washed with 70% ethanol, and resuspended in $\frac{1}{2}$ × TE buffer (5 mM tris, 0.5 mM EDTA). We confirmed successful extraction using agarose gel electrophoresis (1% agarose, 1 × TAE buffer, 100 V, 15 min), assessed DNA quantity and quality using a NanoDrop™ spectrophotometer, and prepared working stocks of 50 ng/µL.

Microsatellite genotyping

We used multiplex PCR to amplify the six microsatellite loci developed by Feldheim, McFarlane & Bowie (2006). Reactions contained 5 µL Platinum^® PCR Multiplex Mix (Applied Biosystems, Foster City, CA, USA), 25 ng template DNA, 0.015 M of each primer (Table 1), made up to a final volume of 10 µL using ultra-pure H₂O. Reactions were amplified in an Applied Biosystems™ 2720 thermocycler under the following conditions: long denaturation at 95 °C for 10 min followed by 33 cycles of denaturation at 95 °C for 30 s, annealing at 56 °C for 1 min 30 s, and extension at 72 °C for 30 s, with a final extension step at 60 °C for 30 min.

We used agarose gel electrophoresis to confirm successful PCR amplification (2% agarose, 1 × TAE buffer, 100 V, 15 min), and conducted fragment analysis on an ABI3500xl Genetic Analyzer with a Liz-500 size standard (Applied Biosystems) at the DNA Sanger Sequencing Facility, University of Pretoria. We analysed electropherogram results using GeneMarker^® v1.95 (SoftGenetics, State College, PA, USA) and designed separate genotyping panels for each Sugarbird species. We checked the consistency of fragment assignment by replicate amplification and scoring of 20% of samples.

ND2 primer design

Sorenson et al. (1999) developed a range of universal avian mitochondrial primers, however, most of these make use of degenerate sites which reduce their utility for sequencing. We used the Primer3 v.2.3.7 (Untergasser et al., 2012) plugin in Geneious^® vR10.2.2 to design universal avian ND2 primers against Taeniopygia guttata (NC_007897.1) with comparative alignment across 16 other bird species (Table 2). Primer L3977 is within the Methionine tRNA gene (positions 3958–3977 in T. guttata). Primer H5191 is within the Asparagine tRNA gene (positions complementing 5211–5191 in T. guttata), resulting in an expected product length of 1254 bp (1213 bp target sequence, including 23 bp of the Methionine tRNA, the complete 1,041 bp of ND2, 70 bp Tryptophan tRNA gene, 69 bp complementary to the Alanine tRNA and 10 bp non-genic nucleotides). We authenticated our primers by PCR amplification across 11 species, spanning the avian tree of life (Barker et al., 2004; Jetz et al., 2012; Prum et al., 2017) and used NCBI BLAST (Altchul et al., 1990) to confirm successful gene-targeting for a subset of individuals following Sanger sequencing. By designing primers in the conserved tRNAs flanking ND2, we were able to amplify the entire ND2 coding sequence across a broad range of species using a single protocol.

Mitochondrial gene sequencing

Mitochondrial genes were amplified in 10 µL PCR reactions comprised of 1X PCR buffer, 1.5 mM MgCl₂, 0.2 mM dNTPs, 0.4 µM forward and reverse primer (Table 1), 0.5 U/µL Supertherm^® Taq polymerase, and 2 µL of template DNA, made up to final volume with ultra-pure H₂O. Reactions were placed into an Applied Biosystems™ 2720 thermocycler for amplification under the following conditions: long denaturation at 94 °C for 5 min, followed by 33 cycles of denaturation at 94 °C for 30 s, annealing at 56 °C for 30 s, and extension at 72 °C for 45 s, followed by 2 cycles of denaturation at 94 °C for 30 s and extension at 72 °C for 1 min 30 s. Successful amplification was confirmed using agarose gel electrophoresis (2% agarose, 100 V, 15 min), and DNA precipitation was performed to remove unincorporated dNTPs and primers in preparation for sequencing.

Cycle-sequencing reactions each comprised of 1 µL BigDye^® Terminator v3.1 Ready Reaction Mix, 1 µL BigDye^® Sequencing Buffer, 3.2 µM forward primer (Table 1), and 2 µL template DNA, made up to a final volume of 10 µL with ultra-pure H₂O. Reactions were placed into an Applied Biosystems™ 2720 thermocycler for amplification under the following conditions: long denaturation at 96 °C for 3 min, followed by 30 cycles of denaturation at 96 °C, and annealing and extension at 60 °C for 4 min. Following DNA precipitation, cycle-sequencing products were sequenced on either ABI3730xl or ABI3500xl DNA Analysers at the Sanger Sequencing Facility, University of Pretoria. The process was then repeated using the reverse primer to confirm sequences and increase coverage.

Data analysis

We used Cervus 3.0.7 (Kalinowski, Taper & Marshall, 2007) to measure allele frequencies at the six microsatellite loci, Genepop v4.6 (Rousset, 2008) to test for deviations from Hardy-Weinberg and linkage disequilibrium among loci, and Micro-checker v2.2.3 (Van Oosterhout et al., 2004) to identify and correct for null-alleles. We conducted parentage analysis using Cervus 3.0.7 (Kalinowski, Taper & Marshall, 2007) in order to identify parent–offspring trios. We conducted parent-pair simulations for 10,000 offspring under the conservative assumption that 25% of potential fathers had been sampled after trial runs with different parameters, allowing two mismatches to account for the Z-linked locus (Pro86) as well as any null-alleles or misscoring.

We used Coancestry V1.0.1.7 (Wang, 2011) to calculate pairwise relatedness estimates (Queller & Goodnight, 1989) and inbreeding coefficients (Lynch & Ritland, 1999) for both sugarbird populations, and used R-Studio V1.1.383 (Ross et al., 1996) to graph the frequency distribution of pairwise relatedness estimates. We then estimated the effective size of both populations using the linkage-disequilibrium method (Waples & Do, 2010) implemented in NEEstimator V2.01 (Do et al., 2014) with jack-knifing across loci. Given that P. cafer exhibits high rates of EPP and given that rates of EPP have not yet been measured in P. gurneyi, we estimated effective population sizes assuming both random and monogamous mating.

Sampling of P. cafer was targeted around nests, and thus our P. cafer dataset included a higher proportion of close relatives (parent–offspring, siblings, and half-sibling). In order to investigate whether sampling methods may affect our estimates of effective population size, we used Friends and Family (De Jager et al., 2017) to produce a reduced P. cafer dataset with a lower overall mean relatedness. This analysis identified groups of potentially related individuals (relatedness > 0.25), and randomly removed individuals from each group until sample size was equal for both datasets.

We used MEGA 7.0.26 (Kumar, Stecher & Tamura, 2016) to align and trim our mitochondrial sequence data and to calculate nucleotide diversity. We then converted our sequence data to NEXUS format using PGDSpider v2.1.1.2 (Lischer & Excoffier, 2012), and used Popart v1.7 (Bandelt, Forster & Röhl, 1999) to construct a median-joining ND2 haplotype network. In addition to the sequence data that we produced, complete ND2 coding sequences of both P. cafer (accession number DQ125990) and P. gurneyi (accession number GU16832.1) were downloaded from GenBank and included in our analyses.

Microsatellite diversity

All six microsatellite loci were highly polymorphic in both P. cafer and P. gurneyi (Table 3). There were no significant deviations from Hardy-Weinberg equilibrium within loci but there was a general tendency towards a slight deficiency of heterozygotes across loci (10 of 12 comparisons). Numbers of alleles and heterozygosity were slightly but not significantly higher in P. cafer than P. gurneyi (paired difference t-test t = 1.2, df = 5, N.S.). We suspected null-alleles at the marker Pro66 due to excess homozygosity and failed amplification of this locus in one sample, and so we used Micro-checker to obtain corrected genotypes at this locus.

Additionally, 1–2 base pair shifts were observed at several Pro66 alleles (4 bp repeat motif), these scored consistently across replicate PCRs and parent–offspring comparisons and were considered as separate alleles in analyses. We also observed a slight overlap in allele size ranges between markers Pro25 and Pro86 in the P. cafer panel, both of which are labelled with the same fluorophore. Fortunately, most individuals with alleles in the overlapping region were identified as females, which are hemizygous for the Z-linked Pro86, and so manual scoring was possible.

ND2 primer authentication

As shown in Table 2, the region to which the forward primer L3977 binds is highly conserved across all 17-species included in our alignment, with only one variable site. The region to which H5191 binds is less conserved, with six variable sites across 17 species, however, all 11 of the bird species included produced clean bands following PCR amplification with the avian ND2 primers L3977 and H5191. Of the 5-additional species selected for Sanger sequencing (A. hoeschi, A. paradiseus, C. ansorgei, S. camelus, and P. robustus), each sample yielded over 1100bp of high quality sequence data (GenBank accession numbers; MG972851, MG972852, MG972853, MG972854, and MG972855).

Mitochondrial diversity

We obtained 1,150 bp, including the complete ND2 coding sequence, from 10 P. cafer (GenBank accession numbers; MG972856 –MG972865) and 15 P. gurneyi (GenBank accession numbers; MG972866 –MG972879). These were combined for phylogenetic analysis with the single existing ND2 sequences of each species from GenBank, each of which was a singleton haplotype. MtDNA diversity matches the pattern from microsatellites, with higher nucleotide and haplotype diversity observed in P. cafer (π = 0.0024; H = 0.84) compared to P. gurneyi (π = 0.0012; H = 0.76), despite a smaller sample of the former. We observed eight haplotypes in P. cafer and 6 haplotypes in P. gurneyi (Fig. 2).

Inbreeding and relatedness

Mean relatedness and inbreeding coefficients were low for both species. Bell-shaped curves were obtained for both relatedness frequency distributions (Fig. 3), with mean relatedness being slightly lower in P. gurneyi (mean  = − 0.022; variance = 0.040) than in P. cafer (mean  = − 0.016; variance = 0.030) with wide variance. Mean inbreeding coefficients were slightly higher for P. gurneyi (mean F = 0.028; variance = 0.026) than for P. cafer (mean F = 0.024; variance = 0.013). We also used Cervus (Kalinowski, Taper & Marshall, 2007) to confirm four parent–offspring pairs in our P. gurneyi dataset and six parent–offspring pairs in our P. cafer dataset, although we suspect that the proportion of second order relatives is significantly higher in our P. cafer dataset due to sampling of some half-sibling fledgelings.

Effective population size

The effective population size of P. cafer in Helderberg, Western Cape, was estimated to be 99 individuals (95% CI [66–182]) assuming random mating and 198 individuals (95% CI [132–356]) assuming monogamous mating. The effective population size of P. gurneyi in Golden Gate Highlands National Park, Free State, was estimated to be 133 individuals (95% CI [55–983]) assuming random mating and 223 individuals (95% CI [111–1,627]) assuming monogamous mating. The reduced P. cafer dataset, filtered for relatives, yielded somewhat higher estimates of effective population size, of 157 individuals (95% CI [87–582]) assuming random mating, and 316 individuals (95% CI [175–1,144]) assuming monogamous mating.

Differences in genetic diversity

All measures of diversity were higher in P. cafer than in P. gurneyi. Several factors could contribute to this disparity. Firstly, our sample of P. cafer was somewhat larger than that of P. gurneyi. This is unlikely to be a substantial contributor to interspecific differences in diversity as substantial samples (>48 individuals) were analysed from each species, each from a single population at the centre of the species range. Some of our measures, such as heterozygosity and nucleotide diversity, are robust to variation in sample size (Nei, 1987; Pruett & Winker, 2008; Hale, Burg & Steeves, 2012). Although the proportion of close relatives was low in both datasets, this was higher in P. cafer, which would tend to reduce observed diversity in that species, contrary to the observed differences. Additionally, these differences in diversity remained after randomly subsampling P. cafer to compare equal sized samples.

Secondly, there may potentially be some degree of ascertainment bias—we used microsatellite markers which were developed in P. cafer and selected for high allelic diversity (Feldheim, McFarlane & Bowie, 2006); therefore, there may be lower allelic diversity at these markers in P. gurneyi due to increased mutational stability (Huang et al., 2016). All microsatellite markers showed high diversity in both species, with extensive overlap in allele size ranges and a low and similar frequency of null alleles. Ascertainment bias would not affect differences in mtDNA diversity.

Diversity may be influenced by the promiscuous breeding system in P. cafer, which shows a high level of extra-pair paternity despite long-term social monogamy (Henderson, 1998; McFarlane et al., 2009. Promerops gurneyi also shows long-term social monogamy, however the frequency of extra-pair paternity may be lower in this species, where territories are more dispersed and reduced sexual dimorphism implies lower levels of sexual selection than in the extremely dimorphic P. cafer. We believe that the best explanation for higher diversity in P. cafer may be higher abundance and higher connectivity across the species’ distribution, in the Proteaceae rich fynbos biome.

Inbreeding and relatedness

Our panel only included six microsatellite markers, but all were highly polymorphic and effective in assigning parentage. The frequency distribution of relatedness estimates yielded bell-shaped curves for both populations (Fig. 3), suggesting that our sampling was not biased by related individuals. Mean relatedness was low for both populations, and slightly higher for P. cafer (−0.016) compared to P. gurneyi (−0.022). Mean inbreeding coefficients were low in both populations, with little variation in the degree of inbreeding among individuals. We detected slightly higher rates of inbreeding in P. gurneyi (mean F = 0.028) compared to P. cafer (mean F = 0.024), which may be attributed to higher overall diversity in P. cafer. In sum, we conclude that it is highly unlikely that inbreeding poses a threat to the persistence of either P. gurneyi in the Maloti-Drakensberg or P. cafer in the south-western Cape. However, it remains to be seen whether the same can be said of the several smaller disjunct regional populations of P. gurneyi.

Effective population size

Sugarbirds form socially-monogamous breeding pairs which often re-unite in subsequent breeding seasons, however, while P. cafer exhibits high rates of EPP (>70%), the breeding system of P. gurneyi remains poorly understood (Henderson, 1999; Calf, Downs & Cherry, 2001; Hockey, Dean & Ryan, 2005). We predict that due to competition for territories and sexual selection the mating system of Sugarbirds is not random, nor is it purely monogamous. The linkage disequilibrium method implemented in NEEstimator V2.01 requires specification of either a random or monogamous mating system, we estimated the effective size of both populations assuming either random or monogamous mating.

The effective size estimated for each population was considerably larger than our sample size, suggesting that both populations were derived from relatively large and stable historical populations. The effective population size of P. cafer was estimated as 99 individuals (95% CI [66–182]) under the assumption of random mating, and 198 individuals (95% CI [132–356]) under the assumption of monogamous mating. The effective population size of P. gurneyi was estimated as 133 individuals (95% CI [55–983]) under the assumption of random mating, and 223 individuals (95% CI [111–1,627]) under the assumption of monogamous mating. Estimates of effective size based on linkage are upwardly biased by sample size, with confidence intervals strongly influenced by the number of loci analysed. Estimates of effective population size were larger for P. gurneyi than for P. cafer, perhaps due to the slightly higher proportion of relatives in the P. cafer dataset. There are wide and broadly overlapping confidence intervals around all these estimates. These figures may be understood as indices of local abundance and relatedness, with values larger than sample size suggesting that additional sampling is required to attain stable estimates of long-term effective population size. In contrast to sampling of P. gurneyi, sampling of P. cafer was targeted around nests, and therefore contained a slightly higher proportion of relatives. Filtering the P. cafer dataset for relatives yielded a much higher estimate of effective population size, with 95% upper confidence intervals comparable to those estimated for P. gurneyi. This demonstrates that slight differences in sampling strategies can affect interspecific comparisons of effective population size, specifically when using methods based on linkage-disequilibrium. In the future it may be informative to use coalescent-based Bayesian approaches, which are independent of sample size but generally give wide confidence intervals for microsatellite data and may demand a better understanding of population structure and gene-flow within species.