Genetic structuring of the coastal herb Arthropodium cirratum (Asparagaceae) is shaped by low gene flow, hybridization and prehistoric translocation

We examined the genetic structuring of rengarenga (Arthropodium cirratum; Asparagaceae), an endemic New Zealand coastal herb, using nuclear microsatellite markers. This species was brought into cultivation by Māori within the last 700–800 years for its edible roots and was transplanted beyond its natural distribution as part of its cultivation. We found very high levels of genetic structuring in the natural populations (FST = 0.84), indicating low levels of gene flow. Reduced genetic diversity was found in the translocated populations, suggesting a large loss of genetic diversity early in the domestication process. The data indicates that rengarenga was brought into cultivation independently at least three times, with the sources of these introductions located within a narrow area encompassing about 250km of coastline. Hybridization was inferred between A. cirratum and the closely related A. bifurcatum, despite A. birfucatum not occurring in the vicinity.


Introduction
New Zealand was the last substantial landmass to be colonised [1], about 700-800 years ago [2]. Pacific Islanders translocated and cultivated many food, medicinal and fibre plants as they colonised the Pacific [3]. Six of these were known to be cultivated in New Zealand by Māori at the time of European contact but their cultivation was marginal owing to New Zealand's relatively cooler climate. Instead Māori began to cultivate several endemic New Zealand plant species for food, fibre and medicine [4]. These cultivated endemic plants are useful for studying the domestication process because their initial cultivation must have been no later than 800 years ago, and they therefore provide a window into the early stages of the domestication process.
Arthropodium cirratum (G.Forst) R.Br. (rengarenga, repihina-papa, maikaika, New Zealand rock lily) (Asparagaceae) is a perennial lily-like herb endemic to New Zealand. It is primarily a coastal species and its natural distribution is thought to be north of around 38˚S (Fig 1) [5]. The fleshy roots of A. cirratum were used as a food source by Māori [6] and its distribution Structure outputs based on microsatellite data from 10 loci. K = 5 was selected as the most likely K value by Structure Harvester but higher levels of K revealed further population subdivision. Outputs for K = 14 and K = 21, which were identified as secondary optima by Structure Harvester are also shown. Individuals are grouped by sampling location and the numbers along the bottom indicate population identification numbers. The two populations marked by stars above the K = 5 plot had the chloroplast haplotype found in A. bifurcatum. The map shows the distribution of sampling sites. Arthropodium bifurcatum sites are identified by (1). For A. cirratum numbers in brackets are the sites' identifying number; see Table 1 for more details. For A. cirratum sites 2-45 are believed to be naturally occurring and sites 46-55 derive from translocations. The position of 38˚S, an important biogeographic boundary in New Zealand, is indicated. The basemap was supplied by Kahuroa.
were included with 1 to 22 individuals per population sampled (Table 1). Many populations were small so large numbers of individuals could not be sampled for these. Six samples of A. bifurcatum, from four sites, were also included; the samples of A. bifurcatum were pooled as one group for all analyses. Approval to collect samples from conservation land was provided by the Department of Conservation (permits WA-23814-FLO, BOP-23814-FLO, TT-23661-FLO and NO-233360-FLO). Permission to collect samples from the Otari Wilton's Bush botanic gardens was provided under permit 145 and the permission of individual land owners was obtained prior to the collection of samples from private land. For samples that did not already have DNA available, a modified CTAB method was used to extract genomic DNA (steps 1, 3-7 from Table 1 in [10]. Twelve microsatellite markers developed for A. cirratum [11] were genotyped. Fluorescent labelling, genotyping and allele scoring followed [11]. Possible scoring errors caused by null alleles, stutter and allelic dropout were assessed with Microchecker v2.2.3 [12]. Matching allele profiles between different individuals were detected using GenAlEx 6.5 [13]. Descriptive statistics including observed (H O ) and expected heterozygosity (H E ), number of alleles (A), number of private alleles (P A ) and inbreeding coefficient (F IS ) were calculated in GenAlEx 6.5 [13].
Loci were tested for deviation from Hardy-Weinberg equilibrium (HWE) with GenAlEx 6.5. Fisher's exact tests of linkage disequilibrium between loci were calculated in Genepop v.4.2 [14]. The Bonferroni correction was used to correct significance values for multiple tests [15]. Genetic differentiation was estimated by calculating pairwise F ST values between all populations with more than five samples with Arlequin v3.5.2.2 [16]. Global F ST for all A. cirratum populations was also calculated with Arlequin. For both global and pairwise F ST values, statistical significance was tested by 10 000 permutations whereby for each permutation individuals were randomly exchanged between populations and a new F ST calculated.
To test for correlations between F ST and geographic distance, we performed Mantel tests between the pairwise F ST values and spatial distance among the nine natural populations with more than five samples in Arlequin v3.5.2.2. Both straight-line distances and minimum coastline distances were tested and significance tested with 1000 permutations. Straight-line geographic distances between populations were calculated using the Geographic Distance Matrix Generator [17] and coastal distances were measured from a map.
Population structure was examined with STRUCTURE v2.3.4 [18,19], without prior grouping assumptions. The number of genetic clusters (K) was set between 1 and 25, with 10 permutations for each. We used the admixture model with correlated allele frequencies and ran 100,000 generations of burn-in followed by 500,000 Markov Chain Monte Carlo (MCMC) iterations. The optimal number of genetic clusters (K) was obtained by calculating the ΔK statistic [20] in STRUCTURE HARVESTER web v.0.6.94 [21]. CLUMPP v.1.1.2 [22] was used to average iterative runs of K and the results visualized graphically with Distruct 1.1 [23].

Microsatellite variability in A. cirratum and A. bifurcatum
In some populations, including natural populations demonstrating variation at other loci, all of the individuals at two of the 12 loci, ArtCir18 and ArtCir38, had fixed heterozygotes. These two loci were excluded from further analyses. Microchecker found no evidence of large allele dropout or stuttering. Possible null alleles were inferred for seven loci, each in a single population. Three of these loci were in the Matapouri Bay population (population 11), which has within-population sub-structuring (see Structure analyses below), likely explaining this result. For the remaining four loci the null allele frequencies was estimated to be low (<0.2) and, because previous research has shown that low frequency null alleles have little influence on the detection of genetic differentiation [24], we retained these loci. Within A. cirratum the 10 loci displayed a total of 87 alleles, with 4 to 20 alleles per locus (mean ± SD: 8.7 ± 4.6). The six A. bifurcatum samples had 23 alleles, with 1 to 4 alleles per locus (mean ± SD: 2.3 ± 1.1). Arthropodium bifurcatum had three private alleles. Within A. cirratum, 12 populations had private alleles with one to three private alleles in each (Table 1).
Following sequential Bonferroni correction, no significant linkage disequilibrium was detected among paired loci comparisons. Significant deviation from HWE was observed for two loci, following sequential Bonferroni correction (ArtCir22 in five populations, including the Matapouri Bay population, and ArtCir 26 in two populations).
The mean F IS for Arthropodium cirratum was high and positive (F IS = 0.72±0.09), indicating inbreeding. At the population level, the F IS values varied from -0.23 to 0.57 (Table 1). The mean F IS for A. bifurcatum was 0.28±0.13.
The cultivated populations (populations 46-55) had very low genetic diversity compared to the natural populations. Across the 10 loci only 21 of the 87 A. cirratum alleles were found in the cultivated populations; all of these were also detected in the natural populations. The cultivated populations have retained only 41.6% of the genetic diversity of the natural populations, as calculated from H E .
Genotype matching revealed that all individuals sampled from the cultivated Tora and Wainuiomata populations (47 and 48) shared an identical genotype but this genotype was not detected from any individuals from natural populations. All the alleles found in the Tora and Wainuiomata populations were also found in the natural populations at Haparapara and Motu (populations 38 and 39), except for the fixed allele at locus ArtCir43.
All individuals from the cultivated Miramar and Paekakariki populations (populations 49 and 50) plus all of the South Island individuals (populations 51-55) had an identical genotype and this matched that found in all individuals from the natural populations at Tohora Pirau (42) and Otiki (45) and four individuals from the Hick's Bay population (43).
The only cultivated population that was not monomorphic across the 10 loci was Kairakau (population 46), whose three individuals showed variation at six loci. None of the three individuals from Kairakau matched any other individuals sampled across all loci.  Table 2). The only non-significant F ST values were between the adjacent populations of Whanarua and Waikawa (populations 40 and 41) and the comparisons between the three cultivated South Island populations (populations 52, 54 and 55). Pairwise F ST values between A. bifurcatum and the A. cirratum populations were also high and ranged from 0.14 to 0.93 ( Table 2). The Structure Harvester analysis of ΔK indicated that the optimal K was 5 (Fig 1, S1 Fig) At higher values of K additional population structuring was evident (Fig 1), with many populations assigned to single clusters with high probability. One population, at Matapouri Bay (population 11), was assigned with high probability to two distinct clusters at values of K!18 (see K = 21, Fig 1). Individuals in these two clusters were fixed for a different allele at ArtCir22. Samples at Matapouri Bay were collected from two adjacent headlands approximately 200m apart, separated by a sandy beach where A. cirratum was not present. The genetic split detected by STRUCTURE corresponds to this geographic separation.
Arthropodium bifurcatum did not form a separate cluster until K = 16. At this value of K the Papanui Point A. cirratum population (34), which shares a chloroplast haplotype with A. bifurcatum, was split between the A. bifurcatum cluster and populations 29-33. Individuals from the second population of A. cirratum with an A. bifurcatum chloroplast haplotype, the Maungaraho population (17), were assigned with high probability to the A. bifurcatum cluster until K ! 16, when they formed their own cluster.
At K = 14 and K = 21 the cultivated populations 49-55 clustered with those from East Cape (42-45). At K = 14 the Kairakau population (46) grouped with Lake Rotoehu (35) but it formed an independent cluster at K = 21. At K = 14 the Tora (47) and Wainuiomata populations (48) were split between clustering with the Lake Okataina population (36) and a second cluster that included samples from Motu (38), Tauwhare Pa (37) and some individuals from Haparapara (39). At K = 21 the Tora and Wainuiomata individuals were assigned with high probability to a distinct cluster, with individuals from Lake Okataina (36) also showing some affinity to this cluster.
A Mantel test indicated that there was a moderate significant correlation between straightline geographic and genetic distance (r = 0.421; P<0.005) for the natural populations with more than 5 individuals. This relationship weakened slightly when coastal distances were used (r = 0.396; P<0.02)

Genetic structuring in A. cirratum
Arthropodium cirratum exhibited very strong genetic differentiation between populations within its natural distribution, and some of this structure can be explained by isolation-by-distance. Population structuring was even detected within the Matapouri Bay site, where two headlands 200m apart had different subpopulations. The nuclear microsatellite results confirm that the high level of genetic structuring previously reported [8] is not restricted to the chloroplast genome and indicates that both pollen and seed dispersal is highly restricted in this species. It has been suggested that A. cirratum seeds are dispersed by gravity [25] and/or wind [26]. The species is insect pollinated but also uses delayed autonomous self-pollination, where selfing occurs once the chance to outcross has passed [27]. This pollination strategy may explain the large variation in F IS values among populations; i.e., some populations have sufficient pollinators for outcrossing, while others do not. It is possible that clonal propagation has contributed to the population structuring and low diversity observed within some A. cirratum populations. However, since this species tends to form distinct clumps in different crevices on rock bluffs, the observed patterns cannot solely be explained by clonal propagation and our sampling targeted distinct clumps to avoid collecting clones.
The high structuring at nuclear microsatellites and presence of private alleles in many A. cirratum populations reinforces our previous recommendation [8] that numerous populations need to be conserved to preserve the genetic diversity of this species.

Origin of the cultivated A. cirratum populations
The two chloroplast haplotypes previously detected from the cultivated A. cirratum populations indicated that these populations had been sourced from the Bay of Plenty/East Cape region. However, it was unclear whether there had been multiple introductions from different sites or a single introduction from a single unsampled variable population [8]. The microsatellite data provides additional insight into the origins of the cultivated population and indicates at least three independent pre-European introductions of A. cirratum into cultivation.
The cultivated samples with chloroplast haplotype C (populations 49-55) had identical microsatellite profiles. This genotype was also found at four locations in the natural populations: Onepoto, Otiki, Tohora Pirau and Hick's Bay (populations 42-45). The Tohora Pirau sample has a different chloroplast haplotype, AC [8], thus excluding it as the source of the cultivated plants. The remaining sites, which occur across a 30 km stretch of coastline at East Cape, are the most likely source for populations 49-55.
Cultivated samples with chloroplast haplotype B were more variable at the microsatellite loci than those with haplotype C. The population at Kairakau (46), showed variation at six of the ten microsatellite loci, despite only three individuals sampled. This was considerably more variation than that detected from Tauwhare Pa and Motu (populations 37 and 38), the only natural populations with haplotype B. Those populations only showed variation at two and one of the loci, respectively. The variation found at Kairakau may result from an introduction from an unsampled, possibly extinct, genetically variable population in the vicinity of populations 37 and 38. There is little natural coastal vegetation remaining in the Bay of Plenty [28] and it is likely that A. cirratum was more widespread in this region in the past. Alternatively, the Kairakau population may derive from multiple introductions from different source populations from this region.
The other cultivated populations with haplotype B, Tora (47) and Wainuiomata (48), had identical microsatellite profiles, suggesting a common origin. However, this genotype did not match any of the individuals we sampled from the natural population. STRUCTURE analyses indicated that these cultivated populations are most similar to natural populations 36-41. However, populations 36 and 39-41 do not have haplotype B, excluding them as the source, and populations 37 and 38 are fixed for a different allele at locus ArtCir43 to the Tora and Wainuiomata samples. This suggests the source of Tora and Wainuiomata was an unsampled, possibly extinct, population from the Bay of Plenty, in the vicinity of populations 37 and 38. This source population is unlikely to be the same as the Kairakau population because the Kairakau samples were fixed for different alleles than those from Tora and Wainuiomata at two loci (ArtCir32 and Art43). This also suggests that the Tora and Wainuiomata populations have an independent origin to the plants from Kairakau and are unlikely to derive via steppingstone dispersal from Kairakau. Overall our results indicate at least three independent translocation events out of the Bay of Plenty/East Coast region and suggest that both chloroplast and nuclear data can provide complementary sources of information.

Loss of diversity in cultivated A. cirratum
All of the cultivated A. cirratum populations, except for Kairakau, were fixed across all microsatellite loci, indicating the loss of diversity through a 'domestication bottleneck'. Overall only 41.6% of the nuclear diversity detected in the natural A. cirratum populations was maintained through the domestication bottleneck, as measure by H E . This value is low compared with other perennial crops [29], where an average of 91.4% diversity was retained and is also lower than the average genetic diversity retained though domestication bottlenecks by annual crops (59.9% [30]). As pointed out by Shepherd et al. [8] a number of factors likely contribute to the high loss of genetic diversity through the A. cirratum domestication bottleneck including the high level of population structuring within the natural distribution of A. cirratum, the narrow area from which cultivated plants were sourced and the physical isolation of the translocated populations from the natural populations, which has prevented gene flow from introducing additional wild diversity into the cultivated populations.
The lack of variation within most of the cultivated populations (populations 47-55) may indicate that A. cirratum was cultivated vegetatively through the division of plants. Alternatively the observed low diversity in these populations may be a consequence of low diversity in the source populations, selfing or loss of diversity through the domestication bottleneck (see below) and/or population bottlenecks following the cessation of cultivation. Next-generation sequencing methods that allow screening of many more DNA markers may be able to discriminate between these hypotheses.

A. cirratum and A. bifurcatum
Two populations of A. cirratum had previously been shown to share a chloroplast haplotype with A. bifurcatum. The nuclear microsatellites also indicate a mixed ancestry for these populations. In the STRUCTURE analyses the Papanui Point population (34) was assigned to two genetic clusters: one containing A. bifurcatum and the other comprising geographically adjacent A. cirratum populations. The Maungaraho population (17) clustered with A. bifurcatum in the STRUCTURE analyses until it formed its own cluster at higher K. Pairwise F ST values also indicated a closer relationship between A. cirratum from Maungaraho and A. bifurcatum than any other comparisons (including within A. cirratum).
There are two possible explanations for the patterns observed: shared ancestral variation and interspecific hybridization, and they are often difficult to distinguish because they produce similar patterns of allele sharing [31]. The small population sizes in A. cirratum favours a hybridization hypothesis (shared polymorphism is more likely in large effective population sizes [32]). The assignment of the Papanui Point plants to a geographically adjacent A. cirratum cluster (as well as to the A. bifurcatum cluster) in STRUCTURE also points to an origin involving hybridization rather than shared ancestral polymorphism. For the observed pattern to have arisen through shared ancestral variation, A. bifurcatum must have diverged subsequent to the development of the genetic structuring within A. cirratum. However, the divergent chloroplast haplotype in A. bifurcatum [8] suggests this is not the case. The fact that the nearest A. bifurcatum to Papanui Point is over 200 km away could infer a chance long distance dispersal event in A. bifurcatum or indicate that this species had a wider distribution in the past. Heenan et al. [5] had noted sporadic northern North Island mainland occurrences of A. bifurcatum. Since that publication our own fieldwork has extended the mainland range of this species, discovering highly fragmented A. bifurcatum populations-sometimes comprising single plants still occur scattered as far south as Whangamata and Ngatutura Point. Otherwise A. bifurcatum is a common plant on predator-free offshore islands of northern North Island [5]. In those places it is a feature of the low often scrubby vegetation associated with seabird nests and roosts (i.e. the 'ornithocoprophilous ecosystem' [32][33][34]. As this ecosystem was once present through all of the main islands of New Zealand, it seems likely that the hybridism indicated at Maungaraho Rock and Papanui Point provides further evidence of the historical loss of the flora associated with these seabirds from the mainland. Supporting information S1 Fig. Results of implementing the Evanno method for detecting the number of K groups that best fit the Arthropodium microsatellite data. According to the ΔK, K = 5 represents the optimal structure partition in our dataset, with secondary optima at K = 14 and K = 21. (PDF) S1 Table. Genotypes at 12 microsatellite loci for Arthropodium cirratum and A. bifurcatum. (XLSX)