Population Growth Confounds Phylogeographic Inference in Namaqua Sandgrouse

Abstract

The Namaqua sandgrouse, Pterocles namaqua, is a highly nomadic granivore of semiarid to arid habitats. As a result of nomadic movements in response to rainfall, the size of the breeding population in any one area fluctuates dramatically between breeding seasons. This high mobility in response to spatial and temporal abundance of food resources is expected to result in little population genetic structuring. Namaqua sandgrouse also shows a seasonally predictable partial migration between the southeast and northwest regions of South Africa, and a further possible north–south migration between southwestern South Africa and central Namibia. It is unclear whether birds migrating between these regions breed in only one or both regions. If populations breed in only one region of their migratory range, then population genetic structuring is predicted to occur. This study addresses Namaqua sandgrouse movements with the analysis of mitochondrial DNA control region sequences. In general, little population genetic structure was evident, yet strong signals of population growth were detected. Several populations have private alleles, which is in direct contradiction to the spatial genetic pattern expected under high levels of gene flow. We suggest that the inference of high levels of female gene flow could be an artifact of population growth and that additional loci will allow a greater understanding of Namaqua sandgrouse movements.

The inference of the demographic processes of migration, admixture, isolation by distance, bottlenecks, and population growth has become a central challenge in phylogeographic studies (Emerson et al. 2001; Knowles and Maddison 2002). Such demographic processes, if acting in isolation, would each be expected to generate different spatial genetic patterns. However, each of these processes may occur during limited time periods, or multiple processes may occur simultaneously, and thus, the untangling of multiple historical factors contributing to spatial genetic pattern remains a major challenge (Wakeley and Hey 1998; Knowles 2004). Recent developments in coalescent approaches (Kuhner et al. 1998; Bahlo and Griffiths 2000; Beerli and Felsenstein 2001; Nielsen and Wakeley 2001) have, however, allowed the coestimation of demographic parameters of interest, such as population growth, effective population size, migration, time of isolation, among other parameters. In this article, we use such available computational methods to identify demographic processes that have contributed to observed spatial genetic patterns in the partially migratory and nomadic Namaqua sandgrouse (Pterocles namaqua).

Namaqua sandgrouse inhabit arid to semiarid regions of southern Africa (Figure 1). As obligate granivores, they must drink regularly and are known to congregate in flocks during a synchronized flight to the nearest surface water each morning (Ward 1972). These flocks may cover a distance of more than 50 km a day in their movement between water, feeding, and roosting sites (Lloyd et al. 2000). Namaqua sandgrouse are highly nomadic over most of their distribution, exhibiting large interannual fluctuations in abundance in any one area (Lloyd, Little, Crowe, and Simmons 2001). In addition, this species has a partial migration, restricted to the southern African sub-region, between the southwestern and northeastern regions of its distribution (Malan et al. 1994; Lloyd, Little, Crowe, and Simmons 2001). Populations in the southwest (Namaqualand/Karoo) are generally larger in the summer months (September–April) and often disappear during the winter months (May–August), whereas populations in the northeast (southern Kalahari) are larger in the winter months (Malan et al. 1994; Lloyd, Little, Crowe, and Simmons 2001). This suggests a partial migration between these regions, with a large proportion of birds from both regions residing in the northeast and southwest, during winter and summer, respectively. Breeding occurs in both regions, but it is unclear whether birds moving between these regions breed in only one or both regions (Lloyd, Little, Crowe, and Simmons 2001).

Figure 1

Distributions of the 3 subspecies of Namaqua sandgrouse (light gray shading) in southern Africa and population groups (dark gray polygons) defined in the GENETREE analysis. Collecting locality acronyms are as in Methods and Table 1.

The high mobility of Namaqua sandgrouse is expected to facilitate a high level of gene flow between populations, resulting in little genetic structuring, unless these movements are structured in some manner. Clancey (1979) suggested that, during their nonbreeding season, individuals from the southern Karoo in southwestern South Africa may range up to 1000 km north to overwinter alongside breeding individuals in northern and central Namibia. This behavior could potentially result in the genetic mixing of populations if individuals from the southern Karoo breed again in central Namibia. This study addresses intraspecific mtDNA population structure and gene flow between populations of Namaqua sandgrouse in southern Africa. Our results, in general, contribute to the understanding of population processes in birds inhabiting arid zones of southern Africa and more broadly to the inference of demographic processes from molecular data.

Methods

Study Area and Sampling

The Namaqua sandgrouse, a near endemic to southern Africa, occurs throughout the arid western region of the subcontinent (see Figure 1). We sampled 10 populations during wing shooting hunts from 29 October 1994 to 23 June 2004. The 10 populations sampled (Table 1) were Benfontein (BF), Droëgrond (DG), Graaff Reinet (GR), Kgalagadi Transfrontier Park (K), Papkuil (PK), Springbok (SB), Groblershoop (GH), Soetdoring (SD), Tontelbos (TB), and Tswalu Kalahari Desert Reserve (TS). For the remainder of the article, we refer to these collecting localities as populations. Collecting covered the range of 3 morphologically described subspecies (Figure 1), yet preliminary molecular analysis indicated no corresponding genetically defined lineages (Delport W, unpublished data). Heart muscle was removed from a total of 81 Namaqua sandgrouse from these 10 localities. Tissue samples were stored in a dimethyl sulfoxide (DMSO)–saturated salt solution (5 M NaCl and 20% DMSO), at room temperature, with a ratio of 1 part tissue to 3 parts solution.

Table 1

Sampled populations of Namaqua sandgrouse, sample sizes (n) and genetic diversity statistics calculated from mtDNA control region sequences. p = Number of private alleles, A_d = allelic diversity, and N_d = nucleotide diversity. Alleles correspond to those depicted in Figure 2

Population	Alleles	n	No. alleles	p	Ad	Nd
BF	1, 2, 3, 16	5	4	0	0.90	0.005
DG	1, 2, 3, 5, 19, 20, 21	11	7	2	0.89	0.008
GH	1, 2, 9, 11, 22, 27	7	6	1	0.95	0.009
GR	9, 11, 13, 14, 15, 16, 20	7	7	3	1.00	0.006
K	2, 3, 7, 8, 17, 21, 25, 26	9	8	3	0.97	0.004
PK	1, 3, 10, 11, 23	7	5	2	0.91	0.011
SB	1, 6, 12, 21, 24, 26	9	6	3	0.83	0.016
SD	2, 3, 4, 7, 18, 26, 28	11	7	2	0.87	0.008
TB	1, 2, 18, 29	7	4	1	0.81	0.006
TS	2, 3, 22, 30, 31	8	6	2	0.93	0.007
Overall		81	31	20	0.93	0.008

DNA Extraction, Polymerase Chain Reaction Amplification, and Cycle Sequencing

Genomic DNA was extracted from all samples using a standard phenol-chloroform DNA extraction procedure. Approximately 50 ng of genomic DNA was used as template in a polymerase chain reaction (PCR) (Saiki et al. 1988). Reactions were performed in a total volume of 50 μl containing, in addition to the genomic DNA, 2.5 mM MgCl₂, 1× reaction buffer, 0.2 mM of each of 4 nucleotides, 1.5 U of Super-therm DNA polymerase (Southern Cross Biotechnology, Cape Town, South Africa) and 12.5 picamol of each primer. The primers L16075 (5′ ACATTGGTCTTGTAAACC 3′—Delport et al. 2002) and ND6H (5′ GCAGGTTGAGGTTTGTTGT 3′—Grant 1999) were used to amplify the complete mitochondrial control region of the 81 individuals. PCR conditions were as follows: denaturing at 94 °C for 2 min; followed by 35 cycles of denaturing at 94 °C for 30 s, primer annealing at 58 °C for 30 s, and elongation at 72 °C for 45 s; and finally an extended elongation period of 10 min at 72 °C. The PCR products were purified using a High Pure PCR Product Purification kit (Roche Diagnostics, Basel, Switzerland), according to the manufacturer's instructions. Purified PCR products were cycle sequenced using BigDye (Applied Biosystems, Foster City, CA), under conditions specified by the manufactured. Primers, H522 (5′GCC TGA CCG AGG AAC CAG AG 3′—modified from Quinn and Wilson 1993) and THR-EXT (5′ CCC CTT CCC CCC CCG GCT 3′—Grant 1999), were used to yield a 318-bp sequence, from the 5′ end of the control region, of both the heavy and light strands for each of the 81 individuals. Grant (1999) extracted purified mtDNA and sequenced the entire control region for several sandgrouse species. The sequences generated in this study from genomic DNA were comparable to those generated in the aforementioned study, and thus, we could be certain that we had not sequenced nuclear copies. All sequences used in this study were deposited in GenBank under the accession numbers AY181949–AY181975 and DQ025759–DQ025766.

Sequence Analysis

Sequences were proofread in Sequence Navigator (Applied Biosystems) and aligned using CLUSTAL X version 1.74 (Thompson et al. 1997). Analysis was conducted to gain an understanding of the demographic processes that could have produced the observed data. First, we tested for the presence of selection in the data (Tajima 1989; Fu and Li 1993). Outgroup-based tests of selection (D*/F* of Fu and Li) utilized chestnut-bellied sandgrouse (Pterocles exustus). Significance was evaluated from approximations to the underlying null distribution (Tajima's D) or from simulations of empirical distributions for theta in the interval [2, 20] (Fu and Li 1993). Second, we tested for the existence of population structure with 1) permutations of the chi-square test (Roff and Bentzen 1989) using custom software, 2) the calculation of pairwise F_ST's, using Arlequin (Schneider et al. 2000) and their deviation from the random association of alleles among localities (10 000 permutations), and 3) the construction of allele networks using TCS (Clement et al. 2000). This preliminary assessment of population structure was used to determine the next appropriate analysis.

Because little population structure was detected, we assumed a panmictic model and coestimated population growth and theta using FLUCTUATE (Kuhner et al. 1998). The Markov Chain Monte Carlo search strategy comprised 10 short chains and 2 long chains both with sampling increments of 20, yet of chain lengths 200 and 20 000, respectively. Replicate runs with alternate random seeds produced similar results, thus ensuring convergence on the correct parameter estimates. Because we assumed a panmictic population model, we used GENETREE (Bahlo and Griffiths 2000) to estimate symmetric migration rates between groups of populations. GENETREE was chosen over MIGRATE-n (Beerli and Felsenstein 2001) because the infinitely many alleles model of GENETREE allows for faster estimation times and less homoplasic inconsistencies, especially, in light of the high levels of migration between the specified groups. Typically, MIGRATE-n outperforms GENETREE in the estimation of migration when populations are structured, yet the inverse is true for populations that approach panmixia (Beerli P, personal communication). Populations were grouped by proximity for the GENETREE analysis (Figure 1). Because GENETREE calculates the likelihood surface in the vicinity of the true parameter estimates, it is necessary to provide upper and lower bounds for estimates on both theta and population growth. Here we used estimates of population growth and theta from FLUCTUATE to define the region in which to estimate the theta-migration-growth likelihood surface in GENETREE (theta: 0–20, g: 0–1000, 1 000 000 repetitions). This process of first estimating theta and growth parameters in FLUCTUATE, and subsequently, defining upper and lower bounds on the likelihood surface in GENETREE, allowed for faster estimation times and smoother likelihood surfaces from the GENETREE analysis.

Results

Mitochondrial control region sequences of 318 bp were obtained from the 81 individuals. Mean nucleotide frequencies (A = 31.8, C = 28.1, G = 14.8, and T = 25.3) were similar to that of other 5′ avian control regions sequenced (Baker and Marshall 1997). Of the 318 positions sequenced, 22 were variable and defined 31 alleles, of which 20 were private alleles, only occurring in their respective populations (Table 1). Allelic diversity was high (0.926) relative to nucleotide diversity (0.008), an indication of either a selective sweep or rapid population growth (Slatkin and Hudson 1991). Neutrality tests indicated that the mitochondrial region was selectively neutral (Tajima's D = −1.35, D* of Fu and Li = −1.65, F* of Fu and Li = −1.84, D of Fu and Li = −1.02, F of Fu and Li = −1.36, and all P > 0.05). Preliminary tests of population structure (Roff and Bentzen 1989) indicated that the observed distribution of alleles among populations was not significantly different from random for populations (χ² = 22.81, P = 0.27) and for the Namaqualand/Karoo and Kalahari regions (χ² = 37.3, P = 0.075). Little geographic structuring is evident in the allele network (Figure 2), yet many populations had unique alleles (Table 1). Pairwise F_ST's calculated in Arlequin again showed little evidence of geographic structure; however, some population pairs had F_ST values significantly different from zero and show significant evidence of population structure (Table 2).

Figure 2

Namaqua sandgrouse (Pterocles namaqua) 95% confidence allele tree constructed using TCS (Clement et al. 2000). Thirty-one mtDNA alleles were identified from 81 individuals from 10 populations. The square box of allele 2 indicates the putative ancestral allele as identified in TCS (Clement et al. 2000). The distribution of alleles among localities is defined in Table 1.

Table 2

Preliminary tests of population structure as derived from F_ST values (lower matrix). Significant deviations (upper matrix) from the random association of alleles at localities, as determined by 10 000 permutations, are indicated in bold type

	BF	DG	GR	PK	K	SB	SD	TB	GH	TS
BF		0.836	0.620	0.158	0.911	0.084	0.763	0.426	0.225	0.999
DG	−0.080		0.455	0.205	0.567	0.096	0.372	0.607	0.294	0.617
GR	−0.050	−0.010		0.848	0.592	0.545	0.263	0.927	0.999	0.256
PK	0.080	0.040	−0.070		0.119	0.939	0.045	0.693	0.694	0.047
K	−0.120	−0.020	−0.030	0.070		0.091	0.854	0.527	0.416	0.936
SB	0.120	0.070	−0.020	−0.090	0.080		0.055	0.607	0.695	0.028
SD	−0.540	−0.010	0.020	0.120	−0.050	0.120		0.378	0.397	0.829
TB	−0.040	−0.050	−0.070	−0.030	−0.030	−0.020	−0.010		0.649	0.331
GH	0.030	0.010	−0.070	−0.070	−0.010	−0.050	0.010	−0.080		0.140
TS	0.140	−0.030	0.020	0.160	−0.050	0.180	−0.050	0.020	0.070

We coestimated theta and population growth using FLUCTUATE. The maximum likelihood estimates of theta and population growth were 0.1484 and 587.6, respectively (Figure 3a). The larger estimate of theta under the maximum likelihood model (0.1484), compared with the unbiased estimate (average number of pairwise differences or nucleotide diversity 0.008), is consistent with longer times to coalescence given the observed population growth. The growth parameter in this case simply indicates that growth is positive. However, if inferences of the timing of growth are to be made, an estimate of the mutation rate is required. We calculated changes in population size, given the estimated growth parameter and upper and lower bounds on mtDNA mutation rates. Because much variability exists in mitochondrial control region mutation rate estimates (Zink and Blackwell 1998; Ruokonen and Kvist 2002; Zink and Weckstein 2003), we applied upper and lower bounds of 10⁻⁵ and 10⁻⁸, respectively, that encompass both protein-coding and control region rates of evolution. Given the bounds on mutation rate, the population has either undergone a recent rapid radiation (within the last 2500 generations) or a steady increase in population size over approximately 200 000 generations. Which of these hypotheses is most likely is difficult to infer without independent information on mutation rate.

Figure 3

Results of the maximum likelihood estimation of demographic parameters. (a) Log likelihood surface of theta and growth estimated in FLUCTUATE. (b) Coestimation of migration and theta using GENETREE. M is the migration matrix multiplication parameter set in GENETREE. Because the migration matrix was set as symmetrical and equivalent with a migration rate of 1 between all pairs of subpopulations, the M in this case represents 2Nm_ij, where N = population size and m = migration rate per generation between population i and j. Confidence limits are indicated as contours of the likelihood surface of −2, −5, −10, −20, −30, −50, −100, −150, −200 and −300 log likelihood units. The −2 (innermost) log likelihood contour represents an approximate 95% confidence limit.

Following the FLUCTUATE analysis, we used GENETREE to infer migration rates between groups of populations. GENETREE requires the data to be fully consistent with an infinitely many sites model (Bahlo and Griffiths 2000), which has the advantages of reducing both homoplasy and computation time. Eight alleles and 2 polymorphic sites not consistent with the infinitely many sites model were removed from the data set following suggestions in the GENETREE manual (Griffiths 2000). The coalescent-based tests of population structure using GENETREE showed high levels of gene flow between all populations (Figure 3b), where an overall level of symmetrical gene flow between all pairs of populations is represented. The GENETREE estimated theta (log theta = 4.5) is substantially larger than for that estimated using FLUCTUATE (log theta = −1.9), when assuming a panmictic population (Figure 3). Although the estimated migration rates between subpopulations in GENETREE are high and suggest panmixia, the specification of population structure in the GENETREE analysis clearly has an effect on the maximum likelihood estimate of theta (Figure 3). Theta is expected to be larger in subdivided populations because coalescent genealogies tend to be longer when populations are structured (Hudson 1990), and the alleles available for coalescence are fewer than what would be available given a panmictic population structure.

Discussion

Demographics and Population Structure

Analysis of population structure showed little evidence of the association of alleles with geographic localities. Neutrality tests indicate the absence of selection, which is furthermore supported by several studies that show mtDNA control region polymorphisms are usually neutral (Fry and Zink 1998; Milot et al. 2000; Griswold and Baker 2002). Although little geographic structure is evident, statistical permutations of F_ST indicated that at least some of the populations are significantly differentiated. This result may be in part due to high frequencies of alleles that are unique to particular populations. The significant differences detected with summary statistics have north to south (SD−PK) and southwest to northeast (PK−TS and SB−TS) orientations. Of these orientations, the north to south divergences could be the result of isolation by distance and the geographic isolation of the SD population. Differences between southwest (PK, SB) and northeast (TS) populations would be consistent with the partial migration accompanied by little mixing of flocks at the overwintering sites. Piertney et al. (2000) in their analysis of red grouse (Lagopus lagopus scoticus) mtDNA population structure similarly identified common alleles that were shared among the majority of populations, yet most other alleles occurred in at least 2 of the 12 populations. Similarly, Ball et al. (1988) found a wide geographic dispersion of alleles in red-winged blackbirds (Agelaius phoeniceus). Both these studies infer panmixia because it is unlikely that the same alleles have arisen in different populations by parallel mutations. A southern African migratory species, the red-billed quelea (Quelea quelea), also shows a lack of geographic population structure, in spite of an apparent migratory divide (Dallimer et al. 2003) in central southern Africa. The authors attribute the lack of genetic structure to the potential for interbreeding between populations on either side of the migratory divide during a second breeding attempt after returning to the interior of southern Africa.

It is likely that a high level of gene flow in Namaqua sandgrouse is causing the populations to approach panmixia. However, Namaqua sandgrouse differs in one respect from these former studies that inferred panmixia. The majority of Namaqua sandgrouse alleles (20 of 31) occur in only a single population. These 20 private alleles occur in 9 of the 10 populations, with BF as the only population without private alleles. This high degree of private alleles within populations should define the existence of population structure; however, the population growth signal overpowers its detection. The more recent an expansion event the more difficult it will to be detect population subdivision. Assuming bounds on mutation rate of 10⁻⁵ and 10⁻⁸, we have calculated that the expansion may have been rapid and fairly recent (2500 generations ago) or slow and old starting approximately 200 000 generations ago. Ray et al. (2003) have used simulations to make inferences of the process of range expansion and its effect on population genetic data. Specifically, they have shown that patterns of expansion observed in populations are not only the result of population size, and the age of the expansion event, but also depend on the level of gene flow between neighboring populations at the time of expansion. Range expansions with large levels of gene flow between populations tend to have most of their coalescent events occurring close to the time of expansion and, therefore, produce star-like genealogies. In contrast, range expansions with low levels of gene flow between populations have most coalescent events occurring before migration events. The crucial result of the former, where spatial expansions occurs rapidly due to high levels of gene flow between populations, is that recent demographic expansions with high levels of initial gene flow lead to patterns of molecular diversity that are equivalent to populations that have always been exchanging a large number of individuals with their neighbors (Ray et al. 2003). This would account for the extremely high levels of gene flow inferred using the maximum likelihood methods in this study, even though some populations have private alleles and significant differences detected with summary statistics.

The causes of such an expansion are uncertain, yet Grant and Leslie (1993) provide support for similar trends in mtDNA patterns from southern African arid zone taxa. Low levels of mtDNA diversity and high levels of inferred gene flow between demes are attributed to population cycles that are the result of extreme climatic variation in the region (Grant and Leslie 1993), where organellar diversity is lost more rapidly than nuclear DNA under severe population bottlenecks (Wilson et al. 1985). Grant and Leslie (1993), however, further demonstrate using simulation modeling that metapopulation structure combined with only low to moderate population bottlenecks can achieve similar patterns in loss of organellar diversity. The population structure of arid zone birds, such as Namaqua sandgrouse, combined with intermediate population size fluctuations, and a life history susceptible to climatic effects could account for the observed patterns. Namaqua sandgrouse has a nest that is a simple scrape in the ground (Lloyd et al. 2000) and is highly vulnerable to the elements. Furthermore, breeding success is highly variable with productivity estimates of 0.07–0.15 young per pair per year to 0.26–0.58 over 4 years in the Namaqualand, but 1.13–4.64 young per pair per year over a 14-month breeding season in the Kgalagadi National Park (Lloyd, Little, and Crowe 2001). Southern Africa is characterized by episodic droughts, which potentially could influence the reproductive success of Namaqua sandgrouse and cause population crashes. These factors most likely contribute to the observed patterns of mtDNA diversity in Namaqua sandgrouse. Higher resolution molecular markers have the potential to provide more power for the detection of population structure in Namaqua sandgrouse. In particular, multilocus approaches for the detection of cryptic population structure have been developed (Pritchard et al. 2000) and would certainly assist in the understanding of population movements within Namaqua sandgrouse.

We thank Peter Bassingthwaighte, Rory Kroon, and Andrew Stainthorpe for collecting the samples from SD, GR, and TS, respectively. Furthermore, we would like to thank all the wing shooters who assisted in the collection of material for this research. The research reported in this article was supported by National Research Foundation grants to T.M.C (grant number: 2053582) and P.B (grant number: 2053653). The opinions and views presented in this article are, however, not necessarily those of the National Research Foundation. We thank 3 anonymous reviewers and subject editor, Robert C. Fleischer, for suggestions that allowed us to improve the manuscript.

Author notes

Corresponding Editor: Robert Fleischer