aCentre for GeoGenetics, Natural History Museum, University of Copenhagen,1350 Copenhagen K, Denmark;bAncient DNA Laboratory, School of Veterinary and Life Sciences, Murdoch University, Perth, WA 6150, Australia;cSchool of Biological Sciences, University of Canterbury, Christchurch 8140, New Zealand;

Significance

In New Zealand, nine species of moa (large, wingless ratite birds) went extinct shortly after Polynesian settlement. In this study, we characterize the gene pools of four moa species during the final 4,000 y of their existence and gain new insights into moa biology and their population sizes. Our analyses show that moa populations were large and viable prior to human arrival in New Zealand, and their demise therefore represents a striking example of human overexploitation of megafauna.

Abstract

The extinction of New Zealand's moa (Aves: Dinornithiformes) followed the arrival of humans in the late 13th century and was the final event of the prehistoric Late Quaternary megafauna extinctions. Determining the state of the moa populations in the pre-extinction period is fundamental to understanding the causes of the event. We sampled 281 moa individuals and combined radiocarbon dating with ancient DNA analyses to help resolve the extinction debate and gain insights into moa biology. The samples, which were predominantly from the last 4,000 years preceding the extinction, represent four sympatric moa species excavated from five adjacent fossil deposits. We characterized the moa assemblage using mitochondrial DNA and nuclear microsatellite markers developed specifically for moa. Although genetic diversity differed significantly among the four species, we found that the millennia preceding the extinction were characterized by a remarkable degree of genetic stability in all species, with no loss of heterozygosity and no shifts in allele frequencies over time. The extinction event itself was too rapid to be manifested in the moa gene pools. Contradicting previous claims of a decline in moa before Polynesian settlement in New Zealand, our findings indicate that the populations were large and stable before suddenly disappearing. This interpretation is supported by approximate Bayesian computation analyses. Our analyses consolidate the disappearance of moa as the most rapid, human-facilitated megafauna extinction documented to date.

The causes of Late Quaternary megafauna extinctions continue to be debated (e.g., refs. 1⇓–3). Climate has been invoked as a major factor driving demographic shifts over evolutionary timescales, but it is undeniable that most recent megafauna extinctions occurred in the presence of humans. However, the role of humans in the extinction process differs among continents and the species studied (4), and it has proven difficult to find evidence for a direct causative link between anthropogenic activity and megafauna loss.

Ancient DNA (aDNA) research has contributed significantly to the extinction debate. DNA extracted from fossil material spanning thousands of years can yield insights into the demographic histories of extirpated populations and extinct species in the period leading up to their loss (5⇓–7). However, most efforts have involved continental-scale data where the large geographic distances hamper fine-scale inferences, limiting our ability to determine the causative agents of discrete demographic events. In contrast, island extinctions offer analytical advantages not afforded by studies of widely distributed species. In island endemics, the absence of gene flow from mainland populations allows us to disregard the spatial dimension in the genetic analyses and focus on extinction dynamics through time.

New Zealand is central to the megafauna extinction debate. It was the last major landmass to be colonized by humans and harbored a diverse assemblage of avian megafauna (8⇓–10). Among them were nine species of moa (11): large, wingless ratite birds ranging in size from the ∼12-kg North Island morph of Euryapteryx curtus to the ∼250-kg females of the two Dinornis species (8). Moa inhabited a variety of habitats across the New Zealand archipelago until their extinction shortly after the arrival of Polynesian settlers, estimated at approximately the late 13th century (8⇓–10, 12). The abundance of well-preserved archaeological sites containing evidence of large-scale exploitation of moa (e.g., ref. 13) brings the controversy of the role of humans in the extinction event into sharp focus.

Early claims of environmental changes or poor adaptive abilities of moa as causes for the extinction (reviewed in ref. 8) have now been largely replaced by the view that direct or indirect human impacts—including hunting, fires, and the introduction of exotic species—were the primary drivers (14⇓⇓⇓–18). Ecological modeling suggests that such human-mediated extinction could have happened within 100 y of Polynesian colonization (10). In contrast, it has been argued, based on limited mitochondrial DNA (mtDNA) data, that moa populations had already collapsed before human arrival, as a consequence of volcanic eruptions or diseases, suggesting that humans were just one of several additive factors responsible for the extinction (19).

To address this issue, we investigated the demographic trajectories of four sympatric moa species in the four millennia leading up to their extinction. We genotyped 281 individuals of Dinornis robustus (Dinornithidae), Euryapteryx curtus, Pachyornis elephantopus, and Emeus crassus (all Emeidae) using mtDNA and six nuclear microsatellite markers developed specifically for moa (20, 21). Moa were recovered from five fossil sites within a 10-km radius in North Canterbury on New Zealand’s South Island (Fig. 1). Through a series of genetic analyses, we have gained insights into moa paleobiology, population sizes, and reproductive success, and we have specifically addressed whether moa populations were in decline before Polynesian settlement of New Zealand.

Of 24 taxon–locus combinations, only two significant deviations from Hardy–Weinberg proportions were observed (α = 0.05) and only one after Bonferroni correction (Table S1). No tests for linkage disequilibrium among loci were significant (P > 0.07) except for one instance in P. elephantopus (Moa_MS2/Moa_MA1; P = 0.0005). The linkage could have resulted from a 106-bp homozygous profile in the Moa_MS2 locus, which was accompanied by a 90-bp homozygous profile in the Moa_MA1 locus in five individuals. This finding was probably an effect of the heterochronous data. Using established criteria (20), we did not find evidence of scoring bias from null alleles, allelic dropout, or stuttering in any of the datasets.

The four species differed in their levels of genetic diversity in both the mtDNA and microsatellite data (Figs. 2 and 3 and Table 2). Differences are visualized in the mtDNA haplotype networks (Fig. 2), where the E. crassus gene pool is largely dominated by two haplotypes separated by a single mutation. Haplotype networks for P. elephantopus and E. curtus are provided in Fig. S1.

mtDNA haplotype networks of D. robustus (n = 87) (A) and E. crassus (n = 81) (B), based on ∼340 bp of mtDNA. The color composition of each haplotype is defined by the individual radiocarbon ages of the fossils. See Table 2 for summary statistics.

Demographic history and genetic diversity. (A) Bayesian skyline plot for D. robustus (n = 87), where the y axis depicts the effective female population size multiplied by generation time. Year zero corresponds to the age of the youngest sample at 602 B.P. (B) Expected heterozygosity (HE) for six microsatellite loci, measured across time in the four moa species (n = 188). Data points represent the mean age and mean HE (with SE) of the moa individuals in 1,000-y time bins.

The microsatellite markers were developed specifically for D. robustus (20, 21), and we therefore cannot eliminate ascertainment bias as being responsible for the lower levels of genetic diversity observed in the other species. However, the congruence between the relative diversity levels of the mtDNA and microsatellite data supports the validity of our data (Table 2). Also, ascertainment cannot explain the much lower genetic diversity in E. crassus in comparison with that in the other emeids (E. curtus and P. elephantopus); these three species are phylogenetically equidistant from D. robustus (11).

Genetic Structuring.

We found strong genetic differentiation among the four moa species (Fig. S2). Only one of the 188 individuals with microsatellite profiles was assigned to the wrong species. Morphologically, this tibiotarsus (AV8470 from Pyramid Valley) could only be a female D. robustus, and the individual exhibited a rare combination of alleles for this species. We found no temporal genetic structuring within species, even when manipulating datasets to comprise only the 10 oldest and 10 youngest individuals for each species (Fig. S2). This finding was supported by a lack of isolation by time using the Mantel test (P values ranging from 0.24 to 0.48) and small insignificant FST values across time bins. In D. robustus, <2.8% of the genetic variation was explained by temporal differentiation, and no pairwise FST values were significant (Table S2).

Demographic History.

Bayesian Skyline plots were generated using the software BEAST (22) and the plots from the three emeid species (P. elephantopus, E. curtus, and E. crassus) displayed flat lines with large highest posterior densities (HPDs), likely reflecting a lack of power in the data to model a reliable genealogy, as has been reported for other aDNA megafauna data (4). We found support for a modest population expansion in D. robustus before human settlement of New Zealand. Despite the relatively large HPDs in the skyline plot, we found an increase in population size in ∼13,000 B.P. (Fig. 3 and Table S3), although we were unable to reject the constant population size model as an alternative fit for the data. This increase was supported by the free model of demographic change, which always yielded a posterior distribution with positive values for growth rate. Also, an expand model proved superior to a decline model, which consistently yielded -infinity likelihood values (Table S3).

The results from BEAST were supported by approximate Bayesian computation (ABC) analyses, which encompassed both mtDNA and microsatellite data. Under the free-model approach, the modal value for the demographic change parameter (Nanc/Ncur) was 0.9 (90% HPD interval 0.26–25.84), indicating that the population was unlikely to have changed dramatically in the period 31,700–100 y before the youngest sample (Table 3). The time of onset of demographic change had a wide posterior distribution (Table 3), as expected if the population size were constant or only changed slightly. The effective population size at the time of human arrival was estimated to be 9,200 individuals, but also with wide HPDs (Table 3).

In the model selection approach, we calculated a Bayes factor of 5.25 in favor of the expand model over the decline model. We observed good posterior coverage and high information content of the summary statistics, making us confident that the ABC analyses were informative with regard to model choice and parameter estimation (Table S4).

Despite allowing the effective population size of the mtDNA (Nemt) to remain independent of the autosomal (microsatellite) effective size (Nems) so as to assess whether there was support for an unequal reproductive contribution of the two sexes, Nemt remained close to the theoretically expected value of Nems/4 (posterior mode of Nemt = 3,200; posterior mode of Nems = 15,800). Hence, we fixed Nemt to a quarter of Nems in all reported simulations.

Discussion

No Precolonization Decline in Moa.

Our genetic data have provided a unique source of information regarding community-level megafauna population dynamics preceding the extinction event. We applied a range of methods to analyze mtDNA and microsatellite data and failed to detect any evidence of moa decline, suggesting that the populations were large and viable throughout the Holocene until their sudden loss.

Temporal sampling is required to directly observe changes in genetic variability over time. To date, relatively few aDNA studies have included microsatellite data; those that did have shown that microsatellite analyses can be a powerful tool to detect temporal loss of genetic diversity and changes in allele frequencies (23⇓⇓–26). We observed no loss of genetic diversity in any of the four moa species (Fig. 3). Moreover, we observed a star-like haplotype network in D. robustus, which suggests an increase in population size (27, 28), although the signal may also reflect an artifact of analyzing sequence data sampled across different points in time (29, 30). However, our detailed demographic analyses of D. robustus, using BEAST and ABC, take temporal sampling into account and support a scenario of a slowly increasing population size during the Holocene.

All moa were sampled from within a 10-km radius, and we can therefore disregard geographic structuring patterns and focus exclusively on genetic patterns through time. Our observation of Hardy–Weinberg proportions and no linkage disequilibrium in the moa populations could be interpreted as random mating. However, these samples span several millennia, and, rather, the result reflects that the microsatellite allele frequencies remained approximately constant during the entire second half of the Holocene. This analysis is also supported by the absence of temporal genetic structuring and lack of isolation by time. Our results reflect very low levels of genetic drift, in agreement with a recent mtDNA analysis of the moa genus Pachyornis, which also failed to detect significant demographic shifts during this period (31).

Despite comprehensive genetic analyses of four moa species, we found no genetic signatures of a hypothesized Holocene prehuman decline. The previous study describing this population collapse (19) may have been compromised by some questionable assumptions (discussed in ref. 32). Rather, our results indicate that the D. robustus population increased slowly at the onset of the Holocene, which seems reasonable from an ecological perspective. Moa were primarily forest and shrubland dwellers, with some entering herbfields and the subalpine zone (8, 33), and pollen records from the South Island reflect an establishment of postglacial shrubland in 14,000–10,000 y B.P., followed by podocarp forest expansion during approximately 13,600–7,500 y B.P. (34, 35).

Population Size of D. robustus.

Our data suggest that D. robustus had an effective population size (Ne) of ∼9,200 individuals when Polyneisans reached New Zealand. Although the posterior distributions are wide, this is the modal value from the free-model ABC analysis (Table 3), which encompasses both mtDNA and nuclear microsatellite data and is without any directional restrictions on simulated growth rate. The conversion of effective population size to census population size (Nc) is problematic, especially for extinct taxa with limited biological information, and no Ne:Nc ratios have been published for extant ratites. However, our best estimate of Ne:Nc ratio for D. robustus is 0.4 (SI Materials and Methods). With an Ne of 9,200, we therefore estimate a census size of ∼23,000 D. robustus individuals, which is the same order of magnitude as the 14,100 individuals estimated from ecological modeling (SI Materials and Methods and Datasets S2 and S3). Such a large population size of a heavy and slowly maturing species (36) suggests a larger, panmictic South Island population, rather than a local enclave isolated in North Canterbury. This finding suggests that D. robustus individuals could disperse over long distances and sustain population connectivity in a landscape fragmented by rivers, glaciers, and mountain ranges. It seems reasonable to infer that for flightless birds, larger body size translates into higher dispersal rates, resulting in larger effective population sizes and higher levels of genetic diversity. This hypothesis could also explain why the smallest species of the region (E. crassus) had the lowest observed genetic diversity (Table 2).

The fossil record suggests that moa had skewed sex ratios with an excess of females (8, 37). If the sexes did not contribute evenly to the effective population size, we would expect a deviation from a 1:4 ratio between mtDNA and microsatellite (autosomal) effective population sizes. Indeed, we observe a 1:4 ratio in the ABC analysis of D. robustus, suggesting that putative skewed sex ratios did not cause differential reproductive success between the sexes. Our result could indicate mating competition among females and accordance with the observation of pronounced reverse sexual size-dimorphism (38, 39) and hypotheses of female territoriality (37) in moa.

Genetic Diversity in the Moa Community.

We found differing levels of genetic diversity in the four moa species, although each remained constant through time (Fig. 3 and Table 2). In the four millennia preceding extinction, we found only half the genetic diversity in E. crassus compared to the three other moa species, suggesting either a smaller population size or a previous demographic bottleneck. Unlike E. curtus and P. elephantopus, E. crassus was not present in the wetter western and northwestern South Island during the Pleistocene (40); the more limited distribution likely resulted in smaller population sizes. If a bottleneck was the cause of lower levels of genetic diversity in E. crassus, it could reflect population isolation in a habitat refugium during the Otiran glaciation, during ∼74,000–17,900 y B.P. (41), followed by recolonization of the eastern lowlands, including North Canterbury, when suitable habitat increased at the onset of the Holocene (e.g., ref. 11). Despite markedly lower genetic diversity, E. crassus appears to have been thriving in the Eastern forests, where late Holocene fossils have been found in great abundance (42, 43).

Conclusion

This study has highlighted the paleobiological insights that can be gained from ancient population genetics beyond the traditional analyses of mtDNA. We profiled 281 individual moa of four species in the 4,000 y preceding their extinction by combining radiocarbon dating, mtDNA sequencing, and nuclear microsatellite genotyping. We observed differing levels of genetic diversity between species in the Holocene moa community. For D. robustus, we found an equal reproductive output between sexes and estimated the population size and demography leading up to the extinction event. Interestingly, the moa extinction process did not leave any genetic traces in our data, very likely because it was too short for increased genetic drift to have an effect on the gene pools. Our results do not support a collapse in any of the moa populations in the millennia preceding Polynesian settlement of New Zealand (19). Rather, our detailed analysis of D. robustus indicated that this moa species increased in numbers during the Holocene. When humans arrived in New Zealand they encountered a large and perhaps still increasing D. robustus population with an estimated effective size of 9,200 individuals. From the archaeological record, we know that moa were hunted intensively and that D. robustus disappeared along with eight other moa species within just one or two centuries following human arrival (10). Together, these findings point strongly toward human contact as the only factor responsible for the extinction.

Materials and Methods

Sampling, Extraction, Identification, and Age.

A total of 290 moa fossils from five adjacent Holocene fossil deposits (Fig. 1) were sampled according to established protocols (21, 44). To avoid including more than one sample from each individual, only left tibiotarsi were sampled, except for one right femur (Canterbury Museum, catalogue no. AV 41188) that was not associated with any of the tibiotarsi (based on consideration of fossil site, bone size, and preservation). DNA extractions from 200-mg aliquots of bone powder were performed in a dedicated aDNA facility (Murdoch University) by using a silica-column-based method (21). Genetic sex identifications were taken from previous work (37, 44), as were 158 of the 217 calibrated radiocarbon ages (45); the experimental procedures are detailed in these references. The dates presented in this study are median calibrated ages (years B.P.) using the SHCAL04 curve in OxCal (Version 4.1; Oxford Radiocarbon Accelerator Unit). An overview of the material is presented in Table 1, and individual data are in Dataset S1.

PCR and Authentication.

Two primer sets (185F/294R and 262F/441R) were used to amplify a 337- to 341-bp fragment (excluding primers) of the moa mtDNA control region. This region has proven informative for assessing intra- and interspecific genetic differentiation in moa (11); PCR conditions are described in previous studies (37, 38). PCR products were sequenced in both directions, and samples with DNA sequences that continually yielded ambiguous base calls (9 of 290 samples) were excluded from further analyses, yielding a total of 281 samples. Samples with mutations, appearing less than three times in the overall sequence alignments, were sequenced and observed at least twice from independent PCRs before being accepted. Information on the six microsatellite loci and their PCR conditions are provided in previous work (20, 21). To overcome the challenges introduced by allelic dropout, we followed the strict guidelines set out in ref. 20. We restricted the microsatellite profiling attempts to the 217 radiocarbon-dated individuals.

Summary Statistics.

The mtDNA sequences were aligned in Geneious (Version 4.8.3; ref. 46) and imported into DNASP (Version 5.10; ref. 47) to calculate genetic diversity. Microsatellite genetic diversity was calculated in GenAlEx (Version 6.4; ref. 48). Deviations from Hardy–Weinberg proportions were quantified as FIS (49), and significance was tested with the exact test in GENEPOP (Version 4.0.10; ref. 50). Deviations from linkage equilibrium were also tested with GENEPOP. We used MICRO-CHECKER (Version 2.2.3; ref. 51) to investigate the presence of null alleles and allelic dropout.

Genetic Structure.

To visualize the genetic structure and diversity of the mtDNA sequences in each of the four species, we generated median-joining haplotype networks with NETWORK (Version 4.5; Fluxus-engineering). We used STRUCTURE (Version 2.3.3; ref. 52) on the microsatellite data to perform several analyses to detect genetic structuring within and among species (SI Materials and Methods). The fixation index (FST) was used to quantify intraspecific genetic differentiation based on the microsatellite data. Because all individuals were sampled within a radius of 10 km, geographic differentiation was not relevant to our study, and we instead investigated genetic differentiation across time. Using the age of each sample, we constructed temporal groups by pooling individuals into 1,000 calendar years. FST values between all pairs of groups were calculated (49), and the significance of the differentiation was assessed by a permutation test using FSTAT (Version 2.9.3; ref. 53). A Mantel test is commonly applied to assess genetic isolation by distance, but it should be equally suitable for detection of genetic ”isolation by time.” Based on microsatellite data, matrices of temporal and genetic distances between all pairs of individuals (within each species separately) were generated and tested for correlations by using the Mantel test in GenAlEx (Version 6; ref. 48). The temporal distances were recorded as the number of calendar years that separated two individuals. We assessed the P value of each correlation with 1,000 randomizations.

Demographic History.

We used two overall methods to estimate demographic history: BEAST and ABC. We analyzed the 217 radiocarbon-dated mtDNA sequences in BEAST (Version 1.6.1; ref. 22). JMODELTEST (Version 0.1.1; ref. 54) and the Akaike Information Criterion (55) were used to estimate the most likely substitution model, favoring a HKY+I+G model for all four species. Prior values on parameters associated with this substitution model were estimated in JMODELTEST and incorporated in BEAST by using wide prior distributions. We analyzed each species dataset using the Bayesian skyline model, which allows population sizes to fluctuate freely through time as governed by the data. D. robustus represented the largest of the datasets (most radiocarbon dates; Table 1) and is the species for which a Holocene decline has been claimed (19). We therefore analyzed this species in more detail and tested an additional four demographic models (Table S3). These models were a constant population size model and three different single-change-point models, assuming a constant population size replaced by an exponential change in population size. These three models were (i) a free model, allowing both exponential growth and decline; (ii) an expand model allowing only growth, as expected when more favorable moa habitat became available after the end of the most recent glaciation in New Zealand, 74,000–17.900 y B.P. (41); and (iii) a decline model, allowing only population decline, following the scenario suggested previously (19). Parameter details are provided in SI Materials and Methods and Table S3. BEAST output files were analyzed in TRACER (Version 1.5; ref. 56) after removing the first 10% of the trees as burn-in. The different models were compared with marginal likelihoods after 1,000 bootstraps.

We used ABC simulations to further elucidate the demographic history of D. robustus. Comparable with the BEAST analyses, we used two different approaches: the free model approach, under which the population size in the most recent phase of the demographic history was allowed to change by a factor of 100 (decline or expansion), and the model selection approach, in which we compared two discrete scenarios of historical population dynamics with declining or expanding population size (Table 3) and compared their respective fit to our observed data using 19 summary statistics (Table S5). See SI Materials and Methods for details on the ABC analyses.

Acknowledgments

We thank the Museum of New Zealand, Te Papa Tongarewa (A. J. D. Tennyson), Canterbury Museum (P. Scofield), and the American Museum of Natural History (C. Mehling). We also thank Eske Willerslev, Malene Møhl, James Haile, and Ross Barnett for assistance with the sequencing and sampling; Jayne Houston, Emma McLay, Helen Hunt, and Frances Brigg for technical assistance; Daniel Wegmann, Simon Ho, Christian Anderson, Alexei Drummond, and Che Si Wu for advice on the analyses; and Marcel Giesen (Bell Hill Vineyard) and the Hodgen family (Pyramid Valley) for their support of our research. This work was supported by Marsden Fund of the Royal Society of New Zealand Contracts 06-PAL-001-EEB (to Palaecol Research Ltd.; R.N.H.) and 09-UOO-164 (to Otago University; C.J.). M.E.A. is supported by European Research Council Marie Curie Actions Grant Agreement 300554. M.B. is supported by Australian Research Council Future Fellowship FT0991741.

Footnotes

↵1To whom correspondence may be addressed. E-mail: morten.allentoft{at}gmail.com or michael.bunce{at}curtin.edu.au.

Blood-sucking sand flies from disparate global regions have a predilection for feeding on the marijuana plant (Cannabis sativa), and the findings hint at a potential avenue for controlling sand flies, which can transmit leishmaniasis.