Frontiers of Biogeography Geographic and ecological segregation in an extinct guild of flightless birds: New Zealand’s moa

The nine currently recognized species of moa (Order – Dinornithiformes; Bonaparte 1853) suffered extinction soon after New Zealand was settled by humans. They were the result of an evolutionary radiation that produced a unique guild of birds – giant, and totally wingless species that evolved in the absence of non-volant mammals. of the nine species of moa, along with information on the geographic, topographic, climatic and edaphic characteristics of sites from which moa remains have been recovered, enabled us to test whether their evolutionary radiation truly was ‘adaptive’, producing ecologically distinct species. Randomization, resampling analyses of moa distributions across North and South Islands revealed highly significant geographic and ecological segregation, with different species tending to occupy different islands, regions within islands, or elevations within regions. Quadratic Discriminant Analyses demonstrated niche segregation at even finer scales, including that based on vegetation‑defined habitats and on local climatic, topographic and edaphic conditions. Moa distributions also appear to have been dynamic over time, shifting in their upper elevational limits as climatic conditions changed and vegetative zones shifted upward during the Holocene Epoch. Our ongoing studies are building on the results presented here to explore the temporal dynamics of moa distributions, assess differential responses of moa species to natural and anthropogenic drivers, and determine how these forces may have combined to cause the extinction of moa just a few centuries ago.


Introduction
Adaptive radiations describe the ecological and evolutionary diversification of monophyletic lineages, largely driven by divergent selection from competition among closely related species (Lomolino et al. 2017). They are among nature's most intriguing and informative phenomena because, among other reasons, they demonstrate the complex yet compelling interplay of ecological and evolutionary forces. Some evolutionary radiations may, however, appear non-adaptive -yielding lineages with only negligible segregation of functional niches and ecological character (see Rundell and Price 2009 Evolutionary diversification of New Zealand's now extinct moa (Fig. 1) was not nearly as extensive in numbers as some of the more spectacular cases (e.g., Hawaiian honeycreepers and lobeliads (Callicrate et al. 2014, Givnish et al. 2009); Madagascar's lemurs and other Malagasy animals and plants (Wirta et al. 2008, Reddy et al. 2012, Herrera 2017 or the cichlids of East Africa's Rift Valley Lakes (Salzburger et al. 2014, McGee et al. 2016). Nonetheless, it produced one of the most anomalous assemblages of native vertebrates known to science -giant wingless birds that converged on the large herbivore niches typically filled on continents by ungulates. The ancestors of this lineage were volant birds that colonized the islands around 58 Ma (Mitchell et al. 2014). The lineage likely underwent alternating periods of diversification and extinction, surviving the Oligocene marine transgression (~ 22 Ma) and its consequent reduction of New Zealand's land area, with moa persisting as a remnant, and possibly monospecific lineage (Bunce et al. 2009; but see also Tennyson et al. 2010, which suggests two moa species may have persisted at this time). The current diversity of moa, thus, appears to have arisen relatively recently, being closely associated with uplift of the Southern Alps (5 -8.5 Ma; Fig. 1), which simultaneously created intra-island dispersal barriers and increased diversity of habitats and potential niches (Bunce et al. 2009).
Biologists have speculated on the ecological associations of moa since the earliest descriptions of the species by Sir Richard Owen in 1839 (see review by Worthy 1990). The earliest inferences were based on Figure 1. Comparisons of the distributions and body sizes (mass and body length) of all nine currently recognized species of moas (modified after Bunce et al. 2009; mass estimates after Latham et al. 2020).
Frontiers of Biogeography 2021, 13.4, e53416 © the authors, CC-BY 4.0 license 3 morphology, suggesting that the mass and structure of their beaks, necks and cervical vertebrae indicated moa were broadly distributed across the "rhizophagous habitats … dislodging the farinaceous roots of the ferns that grow in abundance over the soil of New Zealand" (Owen 1879). Later investigations would include, along with morphological characteristics of the species, analyses of gizzard contents and reconstructions of environmental characteristics generating alternating opinions that moa were browsers, grazers or frugivores -inferring from this general habitat preferences of the guild for shrublands, grasslands or forests, respectively (see Worthy 1990). Anderson (1983) concluded that, based on his analyses of environmental distributions of moa in relation to past climatic, edaphic and vegetative conditions, that the species overlapped to the degree that there was little reason to distinguish among the distributions of the taxa or genera. Worthy (1990), however, took issue with this -citing Atkinson and Greenwood's (1989) assertion that the high degree of geographic overlap among species was facilitated by differences in beak morphology and diet (i.e., segregation of their functional vs. spatial niches; see Young et al. 2012).
Assessing all data available at the beginning of the new millennium, Worthy and Holdaway (2002) analyzed the distributions and environmental associations of what were then believed to be 11 species, concluding that three or four species often co-occurred in any region, but that "the dominant species varied from place to place" (Worthy and Holdaway 2002:192). That is, despite much geographic overlap among the species, each was thought to have a preferred habitat -the three principal habitats being broadly defined as either the upland zone, lowland wet forest zone or lowland dry climate zone. Most recently, Wood et al. (2020:13) inferred habitat associations from diets of moa by compiling and analyzing plant remains from moa gizzard contents and coprolites, concluding that "the nine different moa species had distinct habitat and food requirements" (see also Wood et al. 2013).
Here we further investigate this assertion of habitat segregation by utilizing the most extensive compilation of data on moa distributions yet analyzed to directly assess the ecological relations of this unique guild of species. Recent advances in ancient DNA and other genetic analyses enable us to identify all species of moa from subfossil remains (e.g., bones, egg-shells, and coprolites), while other technological advances allow reconstructions of past climates and other environmental conditions throughout the late-Pleistocene and into the Holocene and recent times. Together, these emerging frontiers of paleobiogeography and paleoecology now provide valuable opportunities to assess whether evolutionary diversification of moa was indeed adaptive -producing an assemblage of species whose distributions became segregated at geographic, and more local climatic and habitat scales. Here we test the hypothesis that the nine recent species of moa comprised a guild of ecologically distinct herbivores. We then provide a preview and prospective of our ongoing research on the temporal dynamics in the ranges and ecological affinities of moa by testing the prediction that the upper elevation limits of these species increased as climatic and vegetative zones shifted following the last glacial maximum.

Development of the moa distributional and environmental databases
Fossil and subfossil records for the analysis of moa interspecific segregation were accessed from published radiocarbon ( 14 C) dating studies (663 specimens from Perry et al. 2014a). As per prior reports (Perry et al. 2014a), these radiocarbon dates were calibrated using the 'BChron' R package (Haslett and Parnell 2008) and the SHCal13 calibration curve (Hogg 2013). All calibrated dates are reported as years prior to 1950 AD. We discarded 83 specimens for which there was no reliable taxonomic information, then qualityassessed each fossil (Barnosky and Lindsey 2010), retaining 528 specimens for our analysis (Table S1). These data were then used to describe the incidence of records for each species across all 116 fossil sites with at least one species of moa (hereafter, 'moa sites') (Table S2). We have used published data reporting the fossil records of moa in New Zealand, calibrated using the SHCal13 curve. The average deviation between this and the updated SHCal20 curve is 557 ± 50 years for the moa fossils used in this study (Fig. S5). Deviations in more recent periods (over the last 4-thousand yr) were more limited and, more generally, these deviations are unbiased and therefore unlikely to have contributed to the shifts in upper elevational limits we observed. Given the multitude of earlier studies that used the SHCal13 curve, however, we recommend a formal sensitivity analysis to reconcile any potential differences and biases introduced when comparing results in future analysis using the SHCal20 curve to those using the older curve.
Paleoclimate data were simulated as a 30-year average at an annual (1 year) step from 21 kBP to the present using PaleoView v1.5.1 (Fordham et al. 2017), resulting in biologically meaningful estimates of past climatic shifts (Fordham et al. 2018). Modern climate data were extracted from the University of East Anglia Climatic Research Unit (CRU) data set (Harris et al. 2020) as a 30-year average centered on 1985  at 0.50° resolution. To capture important orthographic elements in New Zealand's climate, we downscaled this baseline (modern) data using a high spatial (0.005°; approx. 500 m 2 ) resolution 30-year average climatology produced by the New Zealand National Institute of Water and Atmospheric Research (NIWA;1970-1999 resampled to approximately 0.25° (20 km 2 ) resolution.
Because data in PaleoView are not available after 1989 (see below), we adjusted the baseline so that it was centered on 1970 using CRU data and then corrected for paleoclimatic changes throughout the last 21000 years. In both cases this was done using the change factor method for downscaling climate data (Wilby and Wigley 1997). We generated six continuous estimates of climate variables between 21 kBP to 1989 AD: the total annual average rainfall (mL), the annual average temperature, and the average monthly maximum and minimum temperatures (°C) for the coolest (July) and warmest (January) months.
The monthly minimum and maximum temperatures were used to calculate net primary productivity (NPP) using the empirical Miami model (Lieth 1973 Where T a is temperature and P is precipitation. We calculated two measures reflecting topography at each moa site -ruggedness and steepness. Average ruggedness of each grid cell was calculated from a digital elevation model available on the New Zealand Land Resource Information Systems Portal (https:// lris.scinfo.org.nz/layer/48131-nzdem-north-island-25-metre/ and https://lris.scinfo.org.nz/layer/48127nzdem-south-island-25-metre/) as the largest inter-cell difference between each pixel and its surrounding cells. We calculated steepness as the proportion of each 0.25° cell with a gradient greater than 20°. We also included three measures of edaphic character at moa sites, including the content (g.kg -1 ) of clay, sand and silt in the topsoil layer (5-15cm deep), downscaled from their native 10 arcsec resolution (Hengl et al. 2017). Table S3 shows geographic, climatic, topographic and edaphic variables at each of 116 sites where moa were recorded.
The potential vegetation types (PVTs) during the late-Holocene prior to human settlement of New Zealand at each moa site were determined based on data from Leathwick (2001), and that available on the New Zealand Land Resource Information Systems Portal (https://lris.scinfo.org.nz/layer/48279-new-zealandpotential-vegetation-grid-version/). These vegetative reconstructions are based on extensive surveys of forest and open-habitat ecosystems, complemented by information on meteorological data across New Zealand. Additional environmental variables, including solar radiation, soil and atmospheric water deficit, soil leaching, slope, and soil parent material and drainage, were included in regression analyses to determine the combination of environmental variables that best predicted vegetative patterns (Leathwick 2001, Leathwick et al. 2003). The final product was a series of GIS layers describing the PVTs across New Zealand (i.e., during the late-Holocene but prior to arrival of humans and their commensals in the 13 th Century; Wilmshurst et al. 2008). We then assigned each of the moa sites to a PVT based on its geographic coordinates in order to assess the vegetative associations of moa species.

Statistical Analyses
We used Resampling STATS for EXCEL (2019 Statistics. com LLC, www.resample.com; after Simon 1997) to test whether the observed patterns of moa co-occurrence among sites was significantly different than would be observed if the species were randomly distributed with respect to each other. We first constructed a matrix describing the frequency of records for each species across all 116 sites with at least one species (columns and rows, respectively; Table S2). Next, we counted the number of sites where species occurred separately, and the number of sites where at least two species co-occurred. We then randomly shuffled occurrence frequencies for each species among sites (shuffling data within species columns), while retaining the total incidence of records for each species as in the observed data. We calculated for this randomized matrix the count of sites with just one species, and the count of sites with at least two species co-occurring. We repeated these randomization procedures 1000 times and analyzed the probabilities that the observed results could have resulted from random distributions of the nine species of moa among the 116 sites. The proportion of simulations with results as or more extreme as the observed data was calculated.
We conducted Quadratic Discriminant Analysis (QDA; see Tharwat 2016) in Excel using XLSTAT (2021, Addinsoft LLC, www.xlstat.com) to investigate and visualize the potential segregation among moa species (dependent variable) across geographic, climatic, topographic and edaphic dimensions ( Table  S3). In comparison to Linear Discriminant Analysis, QDA in XLSTAT is less influenced by non-normality, heteroscedasticity of covariance matrices, and collinearity among independent variables (see Clark et al. 2007). We performed QDA for three alternative sets of the data including the three most common species, the five most common species, and then for data including all nine species. While the latter subset of the data included species occurring on a limited number of sites, it did allow some speculation on the environmental affinities of all moa, including these less frequently recorded species.
We used Quantile Regression Analyses (QRA) in Excel using XLSTAT (2021, Addinsoft LLC, www.xlstat.com) to investigate the temporal dynamics in elevational limits of moa following the last glacial maximum (i.e., including just those records after 15,000 radiocarbon yr BP). We conducted a QRA for all nine species of moa combined, and separate QRAs for each of the three most common species (N = 184 for all species combined, and 42, 36 and 32, respectively for the three most common species -Pachyornis elephantopus, Dinornis robustus, and Euryapteryx curtus). Inspection of residuals and bivariate plots of elevation-on-age revealed two outliers for records of Eurapteryx curtus, so these were removed from the analyses (N = 30). We restricted analyses of elevation-on-age to records dated using post-1980 methods, and regressed elevation of each moa record against its radiocarbon date to calculate the 95% upper quantile (upper elevational limits of moa distributions) in each QRA.
PVTs were utilized to calculate niche breadths (NB) and niche overlap among moa species by comparing the proportion of land area for each PVT to the proportion of occurrences for each species within Where q i and p i equal the proportions of available resources (land area) and occurrence records of this species in PVT i (for I = 2 to 24).

Results
Resampling randomization analyses revealed highly significant, negative associations among moa species (Fig. S1). None of the 1000 randomized simulations of species distributions among the 116 moa sites yielded results as extreme as the observed data. That is, none of the simulations yielded results where the count of sites with just one species of moa approached the observed 95 single-species sites, and none yielded a count of sites where at least two species co-occurred as rarely as the observed 21 multi-species sites.
The highly significant negative association among moa distributions resulted from segregation among geographic and environmental (climatic, topographic and edaphic) dimensions. Moa were spatially segregated among islands, regions within islands, and across elevations (Fig. 2). Quadratic Discriminant Analyses also revealed finer scale, environmental segregation among moa ranges (Fig. 3). The results of the latter analyses, thus, provide a relatively high-resolution assessment of moa distributions ( Table 1) and one that complements inferences drawn from previous studies (Bunce et al. 2009, after Worthy andHoldaway 2002), particularly that based on dietary analyses of coprolites .
Focusing first on the five most common species, the QDA ordination charts of Fig. 3 reveal Dinornis robustus to have been distributed across a broad range of habitats, elevations and climates across the South Island. In contrast, Anomalopteryx didiformis was most common in sites at low to mid-elevations with moderate temperatures, precipitation and NPP on both islands, while Megalapteryx didinus was a mid-to high-elevation species inhabiting cool and wet habitats in steep and rugged terrain, with silty soils of the South Island. Euryapterx curtus was distributed across a range of elevations in habitats with clay or sandy soils and moderate levels of precipitation and NPP on both islands. Pachyornis elephantopus also occurred across a range of low to mid-elevations across the South Island, more often in habitats with sandy soils, and moderate to cool climates with limited precipitation and NPP.
Despite their limited records, it may still be informative to inspect results for the infrequently recorded species of moa (Fig. S2). The distributional requirements of Pachyornis australis on the South Island seems most similar to that of Megalapteryx didinus, while Emeus crassus appears to exhibit an environmental association similar to that of Pachyornis elephantopus. Pachyornis geranoides was intermediate in its environmental affinities, with its six records tending to be located in low elevation sites with moderately warm climates on the North Island. Dinornis novaezealandiae may have had the most distinct habitat affinities among these rare species, which is not surprising given it was exclusively limited to the North Island.
Consistent with the above patterns derived from QDA of moa on environmental variables, reconstructions of potential (natural) vegetation (PVT) across New Zealand during the late-Holocene (but prior to human colonization; Fig. 4; Fig. S3) revealed multi-species patterns of distribution qualitatively consistent with those above. Collectively, the entire guild of moa inhabited a broad range of habitats, while select species or pairs of species tended to discriminate among vegetation zones they inhabited. Similar to the patterns illustrated in Fig. S2, the PVT affinities of Megalapteryx didinus and Pachyornis australis were similar, both peaking in Rimu-miro/kamahired beech-hard beech forest (PVT 12). In contrast, Euryapteryx curtus exhibited its highest incidence in Matai-kahikatea-totara forest (PVT 4), Pachyornis geranoides and Dinornis novaezelandiae were most common in Kahikatea-pukatea-tawa forest (PVT 3), while Pachyornis elephantopus exhibited its highest incidence in low forest woodland and scrubland below treeline (PVT 22) and in Matai-kahikatea-totara forest (PVT 4).
Niche breadths (NB) of moa across PVT-defined habitats were relatively narrow, ranging from 0.22 to 0.54 among the nine species (theoretical maximum = 1.0; Table 2). Contrary to our expectations (see for example Wood et al. 2013), however, NB of moa was not directly related to body mass, with Anomalopteryx didiformis (44 kg) exhibiting the highest NB and Dinornis novaezealandiae (138 kg) exhibiting one of the lowest observed NBs in this study (0.54 versus 0.34, respectively).
Consistent with the results based on quadratic discriminant analyses of species and environmental variables, niche overlap (NO) of moa species across the 24 PVT-defined habitats was significantly lower (species more segregated by habitat) than that expected if these species were randomly distributed with respect to each other (Table 2). Although highly variable depending on the species pair (ranging from 0.000 to 0.639), the mean value of observed NO indices (0.265) was less than half the mean value of 1000 randomized distribution matrices (0.667; P < 0.001). Niche overlap tended to be low (higher segregation) between congeneric pairs of species (NO = 0.152 for Dinornis robustus -Dinornis novaezealandiae, and ranged from 0.017 to 0.060 for pairings of the three Pachyornis species).     Table 2. Niche breadths (NB) and niche overlap matrix of moas based on comparisons of species occurrences across 24 habitats (PVTs). These indices are electivity measures and range from 0 to 1 for perfectly specialized versus perfectly generalized niche breadths, and for complete segregation versus complete niche overlap (N = total number of occurrences of this species across all sites and PVTs; P-values represent proportion of 1000 randomized species distribution matrices yielding lower NB values). None of the 1000 randomized species matrices yielded mean niche overlaps that were less than the observed mean value (0.265; mean niche overlap of the randomized matrices = 0.667; ECOSIM -Gotelli, N.J. and A.M. Ellison 2012. http://www.uvm.edu/~ngotelli/EcoSim/EcoSim.html).

Discussion and Conclusion
Evolution of moa represents one of the most compelling examples of 'reversals in natural selection' (sensu Lomolino 2009) yet described for insular birds -they became flightless in the absence of non-volant mammals (see Trewick 1997, Clark 1964, Dekker 1989, and occupied the niches filled by large, herbivorous mammals (i.e., ungulates) on continents. Although more limited in terminal diversity than other, classic cases of evolutionary divergence of insular linages (e.g., those referenced in the Introduction), the diversification of moa may have been the largest avian radiation in New Zealand and clearly was adaptive. That is, as a group, moa filled the "large herbivore niche", but the species also clearly segregated their ecological niches within this broader niche.
As with adaptive radiations of other insular lineages, evolutionary diversification of moa was shaped by both ecological displacement and by character release as well. It was driven by competitive displacement, not just from other moa, but from over 200 other, relatively small avian species native to New Zealand (see Worthy andHoldaway 2002, Davies 2003). This transformation toward gigantism may have been amplified in its later stages by predatory pressures from Haast's eagle (Harpagornis moorei -the largest eagle known to have existed; Braithwaite 1992), and by release into the niche space of large herbivorous vertebrates following the extinction of dinosaurs and in the absence of non-volant mammals (see Mitchell et al. 2014).
Whether considering the results for just the more common species, or those for all nine species combined, adaptive radiation of moa produced a guild of species that segregated their distributions across geographic and more fine-scale, environmental dimensions of the geographic template and available habitats (PVTs) of New Zealand. The ecological segregation among species may have been even more demonstrable than our results have indicated because our inferences were based on records of occurrences pooled across a very broad and climatically variable temporal scale (the past 20,000 years). It is likely that the geographic ranges and habitat and niche breadths of these species were narrower during particular periods of distinctive climatic and ecological conditions.
While moa represented the largest terrestrial birds in the Pacific, ratites are a globally distributed group of birds that have a tendency towards gigantism and flightlessness (Phillips et al. 2010). Typically, they are not strongly diversified. In the Indian Ocean, however, Madagascar also supported a guild of flightless, terrestrial giants -the elephant birds (Aepyornithidae), whose isolation, origins and evolutionary trajectory Frontiers of Biogeography 2021, 13.4, e53416 © the authors, CC-BY 4.0 license 10 closely parallel those of the moa in New Zealand (Hansford and Turvey 2018). Although apparently fewer in number (four vs. nine currently recognized species of elephant birds vs. moa, respectively), evolutionary diversification of elephant birds similarly evidenced allopatric divergence across elevations and ecoregions, with niche segregation appearing to be most pronounced between congeneric species (Aepyornis hildebrandti and A. maximus; Hansford and Turvey 2018). We suggest that the techniques applied here, coupled with appropriate dating of the numerous egg shells of the elephant birds across the full extent of Madagascar's ecoregions, may reveal a similarly intriguing story of insular diversification, ecological segregation and adaptive radiation. While our ongoing research includes niche modelling of the temporal as well as spatial dynamics of moa ranges, a preliminary analysis here of elevational dynamics of moa distributions during the late-Pleistocene and the Holocene is instructive. The macroecological pattern of Fig. 5 reveals a constrained relationship in the elevational limits of moa since the last glacial maximum. We thus extend Rawlence et al.'s (2012) description of pronounced, post-glacial elevational shifts of Pachyornis australis to three other species (Dinornis robustus, Pachyornis elephantopus and Euryapteryx curtus) and to the moa guild, in general. The increase in the upper elevational limits of moa is consistent with concurrent dynamics in climatic conditions and upward shifts in habitats during the Holocene Epoch. We can only speculate that the most recent shifts in elevational limits of these species (beginning around 4,700 BP) may have been a response to a shift toward hotter and drier summers that occurred in New Zealand around this time with concurrent shifts in vegetative cover (van den Bos et al. 2018), this perhaps later amplified by anthropogenic factors (the introduction of Polynesian rats [Rattus exulans] and dogs [Canis familiaris], and establishment and subsequent expansion of ecologically significant populations of Māori who hunted moa and substantially modified native habitats with fire; Perry et al. 2014b).
A more rigorous assessment of the causes of these elevational dynamics and, more generally, the patterns of range collapse and ultimate extinction of this guild of giant, wingless birds awaits further analyses on the interdependence among temporal-spatial dynamics of climate, vegetation, Polynesian rats, Māori and moa. We are currently incorporating data from these paleoarchives into process-based (theory and data-driven) simulation models, which run at fine temporal and spatial scales and across large geographical extents, opening windows into moa dynamics during the late-Quaternary. A central goal of this modelling is to apply the lessons of past extinctions for conserving rare and endangered assemblages of today (see Fordham et al. 2011Fordham et al. , 2016Fordham et al. , 2020.

Acknowledgments
We are indebted to the countless scientists and other individuals who have contributed to information on the distributions of this intriguing guild of species, along with those estimating and compiling information on the environmental characteristics of moa distributions and their environments over the past 20,000 years. This research was supported by an Australian Research Council Discovery Project (DP180102392) awarded to DAF and MVL.

Author Contributions
MVL conceived this study and, together with ST, conducted the statistical analyses. DAF provided guidance and financial support for the studies, overall, and with ST compiled the datasets utilized in these analyses. J. Wood and J. Wilmshurst provided key insights into the ecological dynamics of moa and the environments across New Zealand, and all four authors contributed to the writing and editing of this paper.

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Figure S1. Results of resampling analyses assessing the probabilities of results of moa species co-occurrence as extreme as the observed patterns. Figure S2. Habitat and niche segregation among all nine species of moa. Figure S3. Description of potential vegetation types. Figure S4. Loadings of environmental variables on the first two dimensions of factor scores. Figure S5. Deviations between the SHCal13 curve and the updated SHCal20 curve. Table S1. Fossil database for all species of moa. Table S2. Incidence of species records across 116 sites where at least one species of moa was recorded. Table S3. Geographic, climatic, Net Primary Productivity (NPP), topographic and edaphic characteristics at 116 moa sites used in this study. Table S4. Description of Potential Vegetation Types reproduced from Leathwick et al.(2012).