The expansion of humans into previously unoccupied parts of the globe is thought to have driven the decline and extinction of numerous vertebrate species. In New Zealand, human settlement in the late thirteenth century AD led to the rapid demise of a distinctive vertebrate fauna, and also a number of 'turnover' events where extinct lineages were subsequently replaced by closely related taxa. The recent genetic detection of an Australian little penguin (Eudyptula novaehollandiae) in southeastern New Zealand may potentially represent an additional ‘cryptic’ invasion. Here we use ancient-DNA (aDNA) analysis and radiocarbon dating of pre-human, archaeological and historical Eudyptula remains to reveal that the arrival of E. novaehollandiae in New Zealand probably occurred between AD 1500 and 1900, following the anthropogenic decline of its sister taxon, the endemic Eudyptula minor. This rapid turnover event, revealed by aDNA, suggests that native species decline can be masked by invasive taxa, and highlights the potential for human-mediated biodiversity shifts.
The advance of human populations across the world is thought to have had substantial impacts on indigenous vertebrate faunas [1,2]. In particular, population declines and extinctions of numerous megafaunal species are believed to have coincided with human arrival in many parts of the world [3,4]. However, it has often proved difficult to distinguish between human activity and other factors such as climatic shifts as potential agents of Pleistocene faunal decline . In recent years, the advent and rapid development of ancient-DNA (aDNA) techniques has contributed substantially to resolving such questions, enhancing our understanding of past population histories of extinct and extant species. Fossil records, together with radiocarbon dating and aDNA analysis, have allowed the reconstruction of demographic trends over millennia, in some cases revealing climate change as a major factor underpinning Pleistocene extinctions [6,7]. By contrast, in regions where human settlement occurred during the relative climatic stability of the Holocene, it seems clear that anthropogenic factors do indeed explain major extinction events [8,9].
New Zealand is an isolated landmass that presents a particularly noteworthy example of anthropogenic biodiversity loss since human arrival around AD 1280 . This relatively recent colonization event represents an informative model for the study of anthropogenic impacts on naive endemic fauna that evolved without predatory terrestrial mammals. For instance, numerous species of large flightless ratites (moa) went extinct within only 200 years of human arrival, owing to overhunting and large-scale habitat destruction [11–13]. Indeed, the archaeological record reveals that the diet of the indigenous Māori population, especially in southern New Zealand, was initially centred on moa and other large vertebrates, including coastal taxa such as sea lions, seals and penguins, prior to a major dietary shift to fish and shellfish after AD 1450 [14–16]. This cultural transition, along with an associated human demographic decline [17,18], is thought to have allowed the recovery of depleted vertebrate populations, especially in southern New Zealand .
Human impacts and environmental change are thought to have facilitated numerous trans-oceanic colonizations of New Zealand by avian species over recent centuries (e.g. silvereyes, stilts, spoonbills, pukeko) . Additionally, aDNA studies of coastal vertebrates have revealed dramatic biological turnover events, involving wholesale extinction and recolonization [21,22]. In particular, extant mainland yellow-eyed penguin (Megadyptes antipodes) and New Zealand sea lion (Phocarctos hookeri) lineages represent recent arrivals from the Subantartic, following the rapid anthropogenic extinctions of their endemic mainland sister lineages. Given the parallel extinction-replacement events recently detected for penguins and sea lions, we hypothesized that New Zealand's archaeological record might retain evidence of additional cryptic biological ‘turnover’ events.
The little penguin genus Eudyptula is endemic to southern Australia and New Zealand . Genetic analyses of mitochondrial DNA revealed two highly divergent (10–14%) lineages within Eudyptula: one endemic to New Zealand and the other comprising birds from Australia and New Zealand's southeast (Otago) [24,25]. More recently, a multilocus genotypic study established that these distinct Australian and New Zealand mtDNA lineages in fact represent two distinct species : Eudyptula minor (endemic to New Zealand) and Eudyptula novaehollandiae (represented throughout southern Australia and along New Zealand's Otago coast). While it has previously been hypothesized that E. novaehollandiae arrived in Otago during the Late Pleistocene, the timing of (and explanation for) this intriguing colonization event remains highly contentious. Specifically, whereas initial genetic studies [24,25] inferred colonization up to 200 000 years ago (ya), multilocus coalescent analyses suggest a much more recent colonization timeframe, between 50 and 1500 ya .
In this study, we use aDNA analysis and radiocarbon dating of pre-human Holocene fossil, archaeological and historical Eudyptula remains to test the hypothesis that the Australian little penguin was absent from New Zealand prior to human arrival. If supported, this recent colonization scenario might represent a previously unrecognized anthropogenic turnover event.
2. Material and methods
(a) Ancient and historic DNA extraction
A total of 146 prehistoric Eudyptula bones were sampled from museum and archaeological collections (electronic supplementary material, table S1), covering the geographical range of extant little penguins in New Zealand (electronic supplementary material, figure S1). Bones were sourced from natural pre-human Holocene fossil deposits and early to middle prehistoric middens (AD 1280–1600). No specimens were available from late prehistoric middens. To avoid sampling an individual multiple times we used only common elements of the same orientation from a single deposit, or sampled bones from different stratigraphic units within a site. Specimens were sampled from sites that have previously yielded aDNA [15,21,22].
aDNA extraction and PCR set-up were carried out in purpose-built aDNA laboratories physically isolated from other molecular laboratories . We followed strict procedures to minimize contamination of samples with exogenous DNA and authenticate aDNA sequences, including negative extraction and PCR controls . No Eudyptula samples had been processed in these laboratories prior to this study. DNA was extracted from 10 to 256 mg (avg. 70 mg) of bone powder following Rohland et al. . For a subset of samples, all extraction steps subsequent to the binding step were carried out as per Rohland & Hofreiter . DNA was also extracted from nine historic samples from New Zealand museum collections with reliable locality data (electronic supplementary material, table S1). Samples were extracted under the same stringent conditions as described for bone samples, using a Qiagen® DNeasy® Blood and Tissue Kit following manufacturer's instructions, incubating overnight at 55°C and a second overnight incubation after adding an additional 20 µl Proteinase K. Feather samples were extracted following Rawlence et al. .
(b) Control region amplification and sequencing
Five primer pairs were designed to cover the partial mitochondrial control region HVRI (393 bp) previously published for this species . Owing to the large number of polymorphic sites, only two adjacent primer pairs successfully amplified fragments; EumiCR3F (5′-TCTGTCCCGTTAWGAGGACTAARC)/EumiCR3R (5′-CCTCCATTRAGWTCAAGTAGYCA) (53 bp–excluding primers) and EumiCR4F (5′-YGGWAGRTACGGRTATGTYTGR)/EumiCR4R(5′-ATAGATAACCTAATCCCTGAARCTGG) (51 bp–excluding primers). Both fragments contain information to allow assignment of specimens to either Eudyptula species (11 informative sites for species assignment out of 24 polymorphic sites overall). Combining EumiCR3F and EumiCR4R allowed amplification of a 107 bp fragment. If amplification of the 107 bp fragment was unsuccessful due to DNA degradation, we attempted to amplify the original shorter fragments with additional M13-tags attached to the primers to improve sequence quality . Forward primers were tagged with M13USP (5′-TGTAAAACGACGGCCAGT) and reverse primers with M13RSP (5′-CAGGAAACAGCTATGACCAT). The small fragments amplified with these two primer pairs were non-overlapping, separated by a 3 bp gap when aligned to the 107 bp fragment.
Each PCR (25 µl) contained 1 M Betaine (Sigma), 4 mM MgCl2 (Life Technologies), 1 × Gold Buffer (Life Technologies), 2.5 mM dNTPs (Bioline), 250 nM of each primer, 2 U of AmpliTaq Gold DNA Polymerase (Life Technologies) and 2 µl DNA. Unsuccessful PCRs were repeated with 4 µl DNA or 4 µl of a 1 : 10 dilution. PCR products were sequenced bidirectionally by the Genetics Analysis Service, University of Otago using BigDye® Terminator Technology and an ABI 3730xl DNA Analyzer (Applied Biosystems®). Long fragments were sequenced using EumiCR3F and EumiCR4R and short fragments were sequences using M13USP and M13RSP. In case of inconsistencies between forward and reverse sequences, caused by post-mortem DNA damage (G-A and C-T transitions) or for sequences with singleton G-A and C-T transitions, PCR and sequencing were repeated at least twice and a majority-rule consensus applied to the independent replicates .
(c) Radiocarbon dating
Nine pre-human Holocene fossil and archaeological specimens from Otago (excavated from undated sites/layers or sites with multiple occupational layers) were radiocarbon dated at Beta Analytic Inc. (USA) using accelerator mass spectrometry (table 1). Dates are reported as radiocarbon ages, based on Libby T1/2 = 5568 years, uncorrected for calendar variation, in years before present (present being AD 1950). Radiocarbon dates were calibrated in OxCal v. 4.2  using the Marine13 calibration curve  and a local ΔR value of −41 ± 39 . Calibrated ages are reported as 95.4% confidence range in years BC/AD. For undated specimens, mean calibrated ages were calculated from published dates for the site/layer the specimen was excavated from. A specimen from Kaikai's Beach (Otago) could not be radiocarbon dated due to insufficient amounts of material remaining after DNA extraction. Because this site contains multiple layers of occupation  and has only one published  but inadmissible radiocarbon date , the age of this specimen could not be determined and was hence excluded from subsequent temporal analyses.
(d) Sequence authenticity and phylogenetic analysis
Contigs were built from short fragments amplified with EumiCR3F/R and EumiCR4F/R and the 3 bp gap coded as missing data. Modern sequences (GenBank accessions KP308908–KP309419) were trimmed to 107 bp and aligned to aDNA sequences using Sequencher® v. 5.2.4 (Gene Codes).
Six sequences (Eumi_a26, a49, a54, a63, a95, hi8) contained ambiguous bases after PCR amplification, consistent with post-mortem DNA damage. After multiple replications ambiguous bases were resolved applying a majority-rule consensus, with no ambiguities remaining in the final dataset. Only seven private G-A and C-T changes, not shared with any other samples, were observed (Eumi_a6, a25, a33, a77, a94, a122, hi5). These private changes were not derived from ambiguities and multiple sequencing consistently resulted in identical sequences, which are thus regarded as rare haplotypes. Overall observed post-mortem DNA damage prior to replication was low and is unlikely to substantially impact on our estimates of genetic diversity of prehistoric and historic samples. Contiguous sequences with 3 bp of missing data (from non-overlapping fragments) were conservatively assigned to haplotypes (manually in Sequencher and by building a UPGMA tree with absolute distances in PAUP v. 4.0 and assigning the sequence to its closest neighbour) resulting in a total of 32 haplotypes for the combined modern and ancient dataset.
We employed the R script TempNet  to reconstruct parsimony networks for the heterochronous sequence data. Haplotype diversity for different time intervals was calculated in Arlequin v. 3.5 . Sequences with missing data were replaced with their assigned haplotypes as TempNet and Arlequin overestimate diversity in the presence of missing or ambiguous data.
(e) Demographic analysis
We inferred demographic histories of E. minor and E. novaehollandiae in Otago using a Bayesian serial coalescent approach as implemented in Bayesian Serial SimCoal (BayeSSC [42,43]). Two different approaches were used for demographic inference. First, to investigate changes in E. novaehollandiae alone, three alternative scenarios were modelled, which represent the three standard minimal models for a single population's history (stasis, gradual change and abrupt change) consistent with an invasion (electronic supplementary material, figure S2a): a constant population (AUS Model 0), a population exponentially increasing since the founding event (AUS Model 1), and a population constant in size since splitting from an ancestral population (AUS Model 2). AUS model scenarios were modelled using 53 bp fragments available for all samples. Second, owing to the short sequence length and resulting lack of power for single population analysis, data from both species combined was used to approximate the time of invasion and rate of expansion of E. novaehollandiae/decline of E. minor (electronic supplementary material, figure S2b). Here we used two scenarios: a single invasion that has maintained a constant proportion of the population since colonization (Combined Model 1), or an invasion followed by an expansion sometime after colonization (Combined Model 2). Both combined model scenarios were evaluated and modelled initially using a sequence length of 53 bp (to include specimens that amplified for EumiCR3F/CR3R only), then using the 107 bp dataset (for which fewer samples were available).
Sequences were assigned to two time bins for analysis: fossil and early prehistoric period (<AD 1450), and middle prehistoric period to present day (>AD 1450; see the electronic supplementary material, table S1 for specimen ages). For all analyses, we employed the K80 mutation model with a transition : transversion bias of 12.14 as determined by the lowest Bayesian information criterion score in jModelTest v. 2.1.7 . We used three different constant mutation rates: 0.53, 0.96 and 1.43 substitutions site−1 Myr−1 (representing the mean, and lower and upper bounds of the 95% highest posterior density interval for estimates of evolutionary rates in Adélie penguins; ) and a generation time of 8.75 years (representing a mid-point based on life-history data [46,47] and the equation: t = a + 1/(1 − s) , where a is average age at maturity (2.5 yr), and s is average adult annual survival (0.82–0.84)). The measures of diversity used for the rejection step were the pairwise differences within the ancient and modern Australian (and New Zealand) lineages. We also repeated analyses using additional summary statistics: Tajima's D and inter-pairwise differences; haplotype diversity and Fst; Tajima's D, Fst and mean inter-group diversity , which in some cases narrowed or broadened the posterior density distributions but had negligible effect on maximum-likelihood estimates (electronic supplementary material, figure S3 and table S3).
Posterior probability density was determined using approximate Bayesian computation (ABC): out of one million simulations per scenario, the 0.5% of prior values which generated results with the greatest similarity to observed values were weighted by Epinechnikov kernelization [49–51]. Overall model likelihood was determined by re-running simulations at the maximum-likelihood estimate of the posterior distributions, and evaluating the resulting Akaike information criterion (AIC) value under a hierarchical framework .
(a) Species assignment of fossil, archaeological and historic Eudyptula specimens
We successfully extracted and sequenced DNA from 128 Eudyptula specimens, including 119 of the 146 specimens obtained from Holocene fossil and prehistoric midden deposits, and all nine historical specimens. A small proportion of these samples (16 of 128; 12.5%) could only be partially sequenced (CR3F/3R) but nevertheless yielded sufficient phylogenetic information for species assignment.
Temporal genetic analyses reveal a striking turnover event in New Zealand's southeastern region (Otago), with the near complete replacement of E. minor by E. novaehollandiae. Specifically, all 119 Holocene fossil and prehistoric midden specimens (older than AD 1600) were phylogenetically assigned to New Zealand endemic E. minor (figure 1). By contrast, recent historic and modern samples from Otago are dominated by E. novaehollandiae (figure 2; electronic supplementary material, figure S4). The earliest aDNA evidence for the historic presence of Australian little penguins in Otago was found from a museum skin collected before 1915. Two other specimens from 1905 and 1906 were assigned to E. minor. By contrast, five of the six museum skin samples from birds collected in 1969 carried Australian haplotypes.
Ancient genetic analysis revealed substantial prehistoric haplotype diversity in New Zealand Eudyptula (H = 0.82, figure 1). Across New Zealand as a whole, modern haplotype diversity is greater than observed for prehistoric samples (when E. novaehollandiae samples are included in this analysis), whereas no such temporal difference in diversity is detected within the New Zealand species alone. Six haplotypes unique to the prehistoric period, four of which are from Otago locations, were not found in the modern dataset, despite similar geographical coverage and a threefold larger sample size for the modern temporal layer. Within the Otago region, our analysis detected 12 prehistoric New Zealand haplotypes (figure 2), of which only two are represented in the modern layer (the two most common New Zealand haplotypes). These Otago E. minor data suggest a haplotype diversity reduction of 35%.
(b) Radiocarbon dating
Results from radiocarbon dating of nine prehistoric Eudyptula samples, for which no proxy dates could be inferred from layer and site information, are presented in table 1. We report the oldest-known dates for fossil Eudyptula (approx. 34 000 years BP; Old Rifle Butts, Cape Wanbrow ), with these Pleistocene E. minor specimens being genetically similar to modern samples. We further report the first radiocarbon date for the late prehistoric Otago site Huriawa (Karitane). The bones dated from this site unexpectedly yielded ages of 1499 BP (1631–1366 BP, 95.4% confidence range) and 1112 BP (1250–1004 BP, 95.4% confidence range) indicating that these are most likely pre-human specimens subsequently incorporated into an archaeological context. All other dates fit well within the context of their respective sites (electronic supplementary material, table S2).
(c) Demographic changes
Using a Bayesian serial coalescent approach, we investigated the demographic history of the Australian little penguin population in Otago comparing three standard scenarios. Of these, a model of constant population size (AUS Model 0) was most supported (by 11 AIC units and 0.9970 support probability; electronic supplementary material, table S4) with an effective population size estimate of approximately 1700 (300–4800, 95% credible intervals; electronic supplementary material, table S3). In order to understand the impact of the invasion process on the invading as well as the resident species, we modelled an additional two scenarios including sequences of AUS and NZ little penguins from Otago. Here a two-step model (Combined Model 2) was supported (by 11 AIC units and 0.9955 support probability; electronic supplementary material, table S4) in which an initial invasion by E. novaehollandiae (step one) is followed by a population expansion in Otago some time later (step two). Colonization of E. novaehollandiae under this model occurred most likely 25 generations ago (15–60; AD 1800, AD 1500–1890) with an effective population size of colonizers of approximately 2900 (140–10 700) individuals. This estimated colonization timeframe broadly coincides with a 40% (10–90%) population reduction for the endemic E. minor (Ne before reduction approximately 2200, 170–12 600). However, there was insufficient power in the data to determine an exact chronological order for these events (i.e. E. novaehollandiae colonization before or after E. minor decline) given the short sequence length analysed in this study. The magnitude of the initial demographic decline in E. minor was difficult to determine due to its signal being obscured by a more substantial reduction of 95% (60–98%) approximately nine generations ago (2–10; AD 1940, AD 1930–2000). Parameter estimates are presented for models using pairwise inter- and intra-specific differences and a mutation rate of 0.53 substitutions site−1 Myr−1. Estimates were largely consistent across different summary statistics and mutation rates. Estimates for all five tested models are provided in the electronic supplementary material, table S3 and model support is presented in the electronic supplementary material, table S4.
aDNA analysis and radiocarbon dating of prehistoric Eudyptula specimens provide clear evidence that the Australian little penguin (E. novaehollandiae) was absent from New Zealand in pre-human times. Bayesian modelling suggests that colonization of E. novaehollandiae occurred within only a few centuries of human arrival in New Zealand and broadly coincided with population declines of the New Zealand endemic penguin E. minor. This newly identified turnover event also highlights that extinction events can potentially facilitate major biogeographic shifts , and adds to an emerging body of evidence for dramatic turnover events in response to anthropogenic impacts [15,21,22].
The underlying causes for vertebrate extinctions in different parts of the world have been much debated , with many such extinctions linked to climatic fluctuations . In New Zealand, however, species turnover occurred during a period of relative climatic stability. Moreover, the now extinct Otago population of E. minor did not represent the southern geographical range-limit for the species, with populations further south (e.g. Stewart Island; figure 1) apparently persisting from pre-human times to the present. Indeed, radiocarbon dating and genetic data together imply that E. minor occupied this region of southern New Zealand since at least the Late Pleistocene (34 kya) until human arrival . This study thus adds to a growing appreciation of human-driven impacts (distinct from climate-driven change), ultimately leading to species turnover [16,19,21,22].
In contrast to previous mtDNA analyses that inferred ancient secondary contact for Eudyptula taxa , the results of the current study support a very recent colonization, younger than 500 ya (between AD 1500 and 1900). This recent timeframe is also supported by coalescent analysis of microsatellite genotypes . Additionally, the current Bayesian modelling suggests an initial decline of E. minor around the time of E. novaehollandiae colonization. While the precise chronology is unknown, it seems likely that the colonization of E. novaehollandiae occurred in the aftermath of the decline of E. minor. Release from density-dependent exclusion has been similarly suggested for other coastal taxa in New Zealand , and seems to provide a plausible mechanism for rapid biological turnover in this region.
In addition to density-dependent factors , competitive (fitness-related) factors may help to explain the dynamics of this biological turnover event. For instance, the ability of E. novaehollandiae to double brood (raise two clutches per season; ) may have given it an advantage over the endemic E. minor and facilitated rapid population expansion within a few hundred years. In a biased-drift model, E. novaehollandiae completely displaces E. minor within 200 years in approximately 75% of simulations, when only 30% of females in the population double brood (data not shown). Moreover, human cultural and demographic transitions in the region during the sixteenth and seventeenth centuries AD may have further facilitated population expansion of E. novaehollandiae in southern New Zealand [17,18].
Comparison of natural and archaeological deposits provides further evidence for a substantial decline in the E. minor population in Otago after human arrival. Little penguin bones occur in 60% of early prehistoric midden deposits (ca. AD 1280–1450), in 50% of those from the middle period (AD 1450–1650), but are completely absent from sites analysed from the late prehistoric period (AD 1650–1800) . It has been suggested that Māori hunted some bird taxa only on an opportunistic basis (as they were encountered), in contrast to targeting specific species. In this way, the abundance of species remains in middens may be indicative of their natural abundance , suggesting that E. minor, once one of the most common coastal vertebrate species, declined within 200 or 300 years of human arrival. In addition to direct exploitation by humans, introduced predators (Polynesian rat and dog) might have also had substantial impact on Eudyptula along with other burrow-nesting seabirds .
Population census data indicate that E. novaehollandiae has indeed undergone a substantial population expansion since its arrival in New Zealand. While the genetic data did not allow a reliable approximation of the size of the contemporary E. novaehollandiae population, a survey of little penguin abundance and breeding distribution for the whole coastal Otago region  estimated the total population to be around 10 000 individuals. The founding size for the Otago population of E. novaehollandiae is unknown. Based on the biology of this species, the genetically derived estimate of approximately 3000 founding birds (suggested by ABC) seems likely to be a substantial overestimate. Nevertheless, as little penguins in Australia have been shown to associate together at sea, hunting in groups  and forming rafts before coming ashore, the possibility of a substantial number of founders should not be discounted. The currently large population size compared with the (probably overestimated) founding population size thus supports a population expansion of E. novaehollandiae since arrival in New Zealand.
In conclusion, our aDNA analyses of little penguins directly document a compelling example of rapid faunal shift coinciding with human expansion into New Zealand. These data highlight that the decline of a native species can be masked by the cryptic invasion of a related species . In contrast to previously documented wholesale extinction events, however, the current study characterizes a regional species turnover event, with two formerly allopatric lineages now occurring in secondary contact. This study thus presents an intriguing case of increased regional biodiversity associated with human impacts, emphasizing that anthropogenic processes can facilitate rapid changes in biogeographic distributions.
No ethics approval was necessary to work with Eudyptula museum specimens.
All aDNA sequences have been deposited in GenBank under accessions KT225747–KT225874.
J.M.W. and S.G. conceptualized the research. S.G. carried out the laboratory research. S.G. and C.N.K.A. performed analyses. S.G., J.M.W. and C.N.K.A contributed to writing the manuscript. N.J.R., I.W.G.S. and R.P.S. provided samples, advice on analyses and archaeological context and critically revised the manuscript. All authors gave final approval for publication.
The authors have no competing interests.
This research was financially supported by the Allan Wilson Centre, Royal Society of New Zealand Marsden Contract UoO1112, and by a University of Otago Doctoral Scholarship.
We thank Auckland Museum, Canterbury Museum, Jill Hamel, Museum of New Zealand Te Papa Tongarewa, Otago Museum, and the University of Otago's Department of Archaeology and Anthropology, for providing samples. Ken Miller supplied maps for figures 1, 2 and S1. We also thank Lisa Matisoo-Smith for access to the Otago University Anatomy aDNA laboratory.
- Received December 1, 2015.
- Accepted January 11, 2016.
- © 2016 The Author(s)
Published by the Royal Society. All rights reserved.