Loss and gain of sexual reproduction in the same stick insect.

The outcome of competition between different reproductive strategies within a single species can be used to infer selective advantage of the winning strategy. Where multiple populations have independently lost or gained sexual reproduction it is possible to investigate whether the advantage is contingent on local conditions. In the New Zealand stick insect Clitarchus hookeri, three populations are distinguished by recent change in reproductive strategy and we determine their likely origins. One parthenogenetic population has established in the United Kingdom and we provide evidence that sexual reproduction has been lost in this population. We identify the sexual population from which the parthenogenetic population was derived, but show that the UK females have a post-mating barrier to fertilisation. We also demonstrate that two sexual populations have recently arisen in New Zealand within the natural range of the mtDNA lineage that otherwise characterizes parthenogenesis in this species. We infer independent origins of males at these two locations using microsatellite genotypes. In one population, a mixture of local and nonlocal alleles suggested males were the result of invasion. Males in another population were most probably the result of loss of an X chromosome that produced a male phenotype in situ. Two successful switches in reproductive strategy suggest local competitive advantage for outcrossing over parthenogenetic reproduction. Clitarchus hookeri provides remarkable evidence of repeated and rapid changes in reproductive strategy, with competitive outcomes dependent on local conditions.

INTRODUCTION

Sexual reproduction is a common strategy among multicellular organisms, despite its two‐fold cost relative to asexual reproduction (Maynard Smith, 1978). This apparent paradox has stimulated research exploring the potential benefits of sexual reproduction that would explain its prevalence (Barton & Charlesworth, 1998; Burke & Bonduriansky, 2017; Otto, 2009). Where evolutionary lineages that reproduce in different ways are in competition, valuable evidence about the selective advantages of different strategies can be gleaned (van der Kooi & Schwander, 2014; Neiman, Meirmans, Schwander, & Meirmans, 2018). For example, the gradual upstream range expansion of sexual Potamopyrgus antipodarum snails recorded over 20 years (Wallace, 1992), suggests short‐term local competitive advantage of sex in this species that is supported by evidence from parasite infection rates (Gibson, Xu, & Lively, 2016; Jokela, Dybdahl, & Lively, 2009). Unfortunately, such systems are ephemeral and thus their availability can be a limiting factor in the study of why sex is so common (Barton & Charlesworth, 1998; Maynard Smith, 1978; Mirzaghaderi & Hörandl, 2016; Otto, 2009). Although it might be expected that a single reproductive strategy within a species would prevail in all circumstances, local factors can determine local adaptive advantage (Tabata, Ichiki, Tanaka, & Kageyama, 2016). Our study with temporal sampling provides evidence that local effects are important to the outcome of competition between reproductive strategies within a single stick insect species.

The maintenance of sexual reproduction in nature has been addressed by contrasting the traits of sexual and related asexual lineages, and ideal study systems have natural replication, with independently derived populations displaying contrasting reproductive strategies (see reviews; Neiman et al., 2018; Neiman & Schwander, 2011). Examples of this circumstance include the snail Potamopyrgus antipodarum whose different clonal lineages are derived from the same sexual species (Gibson et al., 2016; Jokela et al., 2009), and Timema stick insects that contain multiple species pairs of sexual/asexual lineages (Bast et al., 2018; Schwander & Crespi, 2009; Schwander, Henry, & Crespi, 2011). An optimal study system would provide replication but lack the confounding effects of ploidy variation and hybridisation that exist in many study systems (Kearney, 2005). The stick insect species Clitarchus hookeri meets these optimal criteria; replicates of sexual and asexual populations exist, and all individuals assessed so far are diploid (Myers, Trewick, & Morgan‐Richards, 2013). However, the range of mechanisms employed by asexually reproducing organisms can influence our ability to isolate the benefits of sex (Neiman & Schwander, 2011). In stick insects, parthenogenetic reproduction can be apomictic (without recombination; Schwander & Crespi, 2009), or automictic (with recombination; Marescalchi & Scali, 2003) and thus reduction of heterozygosity in subsequent generations. The outcome of automixis is similar to selfing in plants and can cause inbreeding depression (Barrett, 2002). It is not yet known whether parthenogenetic C. hookeri reproduces with recombination and resulting homozygosity or without.

Stick insects (Phasmids) have provided many examples of recently derived asexual lineages (Scali, 2009; Schwander & Crespi, 2009), with the switch from sexual reproduction to parthenogenesis having occurred independently numerous times within this family (Ghiselli, Milani, Scali, & Passamonti, 2007; Schwander et al., 2011). Some lineages contain facultative parthenogens (Mantovani & Scali, 1992) that provide the opportunity to investigate competition models of the two reproductive strategies. Facultative parthenogenetic populations have the potential to switch from asexual to sexual reproduction if males are available, and we can distinguish three scenarios in which asexual lineages revert to sexual reproduction: swamping, introgression, and male genesis. Swamping could follow colonization by sexual males and females, that might involve long‐distance dispersal into an empty habitat patch. The genetic signature of such a population would be distinct from adjacent parthenogenetic populations but have the same alleles and mtDNA haplotypes as found in sexual populations further afield. Introgression involves establishment of males from outside the area. These males would successfully fertilize eggs of local females and so restore sexual reproduction. This scenario would mix local alleles and alleles distinctive to sexual populations, resulting in a population with high allelic diversity relative to surrounding all‐female populations, but retain mtDNA from the local parthenogenetic lineage. Alternatively, successful natural male genesis could occur via in situ nondisjunction of the sex chromosome. Loss of X chromosomes during egg production in some parthenogenetic insects can produce locally derived males (e.g., Brock, Lee, Morgan‐Richards, & Trewick, 2018; Pijnacker & Ferwerda, 1980; Scali, 2009). If these males are fertile and restore sexual reproduction to the population, it will contain only alleles and mtDNA haplotypes from the local lineage, but observed heterozygosity will rise from levels when males were rare (if parthenogenetic reproduction is automictic).

In New Zealand, most northern populations of the common diploid stick insect C. hookeri have an even sex ratio and all reproduction is assumed to be sexual (Langton‐Myers, Holwell, & Buckley, 2019; Morgan‐Richards, Trewick, & Stringer, 2010; Myers, Buckley, & Holwell, 2015; Myers, Holwell, & Buckley, 2017), but males are absent from many southern locations (Figure 1). A widespread parthenogenetic lineage appears to have resulted from range expansion by a single mtDNA lineage (Buckley, Marske, & Attanayake, 2010; Morgan‐Richards et al., 2010). Throughout the southern and eastern North Island and eastern South Island of New Zealand (Figure 2), reproduction is inferred to be entirely parthenogenetic (Wu, Twort, Crowhurst, Newcomb, & Buckley, 2017); males have been recorded only sporadically and at very low density (Buckley et al., 2010; Morgan‐Richards et al., 2010). However, two sexual populations of C. hookeri that are surrounded by all‐female populations are described in the present study. Both of these sexual populations are within the geographic range of the parthenogenetic mtDNA lineage (Figure 1), but feed on plants that had been transported and planted by humans. Surprisingly, C. hookeri is also found in the United Kingdom (UK), where an all‐female population was founded by accidental human‐mediated dispersal (Figure 1).

Figure 1

Three questions arise from the distribution of male and female Clitarchus hookeri stick insects. Parthenogenetic populations surround sexual populations at Otaki and Wilton, so from where did the males come? Human‐mediated dispersal has resulted in a parthenogenetic population in the UK, but from where in New Zealand did this all‐female lineage originate? Blue (sexual) and pink (asexual) shading indicates the natural extent of sexual and asexual populations in New Zealand [Colour figure can be viewed at http://wileyonlinelibrary.com]

Figure 2

Similarity of mtDNA sequences indicates the probable origins of Clitachus hookeri stick insect populations. (a) Sample locations in New Zealand, coloured to highlight two intraspecific mtDNA lineages as shown in b (pink parthenogenetic; blue from Taranaki). (b) Phylogenetic relationships of full mtDNA diversity sampled in New Zealand inferred with Maximum likelihood (COI 1,389 bp), colours as in networks and map, scale represents nucleotide site distances. (c) Parthenogenetic lineage haplotype network (medium joining, COI 1,307 bp) in pink, with the haplotypes of male individuals indicated. (d) Haplotypes from Taranaki lineage (blue) includes UK specimens (median joining, COI 1,307 bp). Circle sizes in c & d are proportional to sample size, and number of nucleotide substitutions that differentiate haplotypes are indicated on branches [Colour figure can be viewed at http://wileyonlinelibrary.com]

Captive breeding experiments show that, if they have mated, C. hookeri females from sexual populations produce all offspring via sex. In contrast, even when provided with males, females from the parthenogenetic lineage experience a barrier to fertilisation. Despite retaining sexual signals (Nakano, Morgan‐Richards, Godfrey, & Clavijo McCormick, 2019) and mating, only about 5% of offspring from these females resulted from sexual reproduction (Morgan‐Richards et al., 2010). This resistance to fertilisation might prevent a successful transition to sexual reproduction via males colonizing the range of the parthenogenetic lineage unless there exists an indirect sons' effect (Kawatsu, 2015; Morgan‐Richards et al., 2010). Here, we use both temporal and spatial sampling of wild populations of the New Zealand stick insect C. hookeri to study three putative transitions in reproductive strategy and attempt to determine whether the outcome of competition between different reproductive strategies is contingent on local conditions. We provide evidence that C. hookeri populations have undergone multiple sexual‐asexual transitions. We utilize mating experiments to test the permanency of asexual populations. We compare genetic variation (mtDNA haplotypes and nuclear microsatellite loci) to determine which New Zealand C. hookeri population is the closest relative of the asexual population established in the UK, and use this relationship to infer the UK population's reproductive history. Genetic variation is used to determine the likely origin of the two sexual C. hookeri populations in New Zealand that are surrounded by all female populations. We distinguish between three possible scenarios to explain the origin of these two sexual populations: swamping, introgression, or male genesis.

MATERIALS AND METHODS

Clitarchus hookeri

The common New Zealand mānuka stick insect (the “smooth stick insect” in the UK) is frequently found on the Myrtaceae species Leptospermum scoparium (mānuka) and Kunzea sp. (kānuka), collectively known as tea‐tree. Adult female mānuka stick insects are about 9 cm long; the males are shorter (~7 cm) and thinner (Figure 1). We collected adult and juvenile Clitarchus hookeri from natural sites around New Zealand to obtain a representative sample of their sexual and geographic diversity (Figure 2, Table S1). In the UK, populations of stick insects have established from human‐mediated dispersal (Jewell & Brock, 2002). Clitarchus hookeri was first recorded in the Tresco Abbey Gardens in 1940, where numerous plant species from New Zealand were planted in 1911 and in the 1930s (Brock, 1987; Godley, 1997). Nine C. hookeri individuals were collected by Paul Brock from the Isles of Scilly UK (five Tresco; four St Mary's; one Bryher Island) and held in captivity in the UK to obtain eggs.

Breeding experiments

To assess the ability of the all‐female UK lineage to return to sexual reproduction in captivity we imported UK C. hookeri eggs to New Zealand to raise and mate with New Zealand male C. hookeri. To achieve this Paul Brock collected 922 eggs from nine C. hookeri (referred to here as generation‐1) from the Isles of Scilly UK in Northern Hemisphere autumn 2013 (August, September, October). On import to New Zealand, these were held in a physical containment facility (PC level 2) for 12 months. We monitored the eggs as they hatched between May 2014 and February 2015 (Southern Hemisphere winter, spring, summer; generation‐2) and raised some of the nymphs, which were housed and fed individually, to maturity. Moving from the Northern Hemisphere to the Southern Hemisphere resulted in many Scilly Isle individuals maturing before adult males were available in New Zealand; this constrained the origin and number of potential mates we provided. Ten of the adult generation‐2 females were each provided with a single male C. hookeri collected from either Karapiro or Kerikeri in New Zealand (January 2015). We observed whether copulation took place and collected eggs that were laid (Table S2). Generation‐3 hatched from these eggs between September 2015 and January 2016 (Table S2). We sexed all nymphs using external morphological differences (Morgan‐Richards et al., 2010; Stringer, 1970), and raised a portion of (all‐female) generation‐3 individuals to maturity. This generation, which was more closely synchronised with wild New Zealand populations, matured later in the season and allowed us to pair more than one adult male per female. We provided 13 adult generation‐3 females with males from the Wilton and Taranaki populations (2–4 males with each female). As with generation‐2, we observed copulation (February and March 2016), collected eggs, and sexed all nymphs that hatched (generation‐4 hatched between September and December 2016). We used the frequency of males that hatched from the eggs to estimate the proportion of offspring resulting from sexual reproduction, assuming that half of all offspring resulting from fusion of sperm and egg (sex) would be male. We genotyped a subset of nymphs to confirm that males were the result of sex rather than nondisjunction of one of their X chromosomes and to test our assumption that half the offspring resulting from sexual reproduction were female (Table S2). Nondisjunction during egg formation is rare but has been recorded among unmated C. hookeri females in captivity, which produce fewer than one son for every 300 daughters (Morgan‐Richards et al., 2010; Salmon, 1955).

Natural populations

We collected stick insects from six locations within the New Zealand range of the parthenogenetic mtDNA lineage (Morgan‐Richards et al., 2010): four locations where males have not been recorded and two locations where males were locally abundant (Figure 3; Table S1). These collections included one sexual population on private property where C. hookeri are seasonally common and feed on various non‐native host plants including plum (Prunus domestica; Lupin Road, Otaki). Here, we collected adult insects in autumn (April 2003) and nymphs in spring (October 2013). At Wilton (Otari‐Wilton's Bush in Wellington), we located a C. hookeri population where males were uncommon in the 1960s (Salmon, 1991) and 2000s (Morgan‐Richards et al., 2010) but were at similar frequencies to females in 2016. In 2003 stick insects were collected from native plants (kānuka) within the Otari‐Wilton's Bush native forest remnant (= Wilton). The recently planted mānuka and kānuka on which the 2016 sample were feeding had been germinated from local seed and raised at a plant nursery in Taupo (300 km north of the site; personal communication Eleanor Burton, Wellington City Council).

Figure 3

Genotypes of Clitachus hookeri stick insects suggest that the UK population was derived from a sexual New Zealand population, that Otaki males result from population mixing, and Wilton males result from in situ male genesis (nondisjunction of X chromosome). Bayesian assignments of individual stick insect genotypes (10 microsatellite loci) to three genetic clusters (optimal model K = 3) is shown by three colours. Collecting locations indicated on the map are coloured to indicate dominant genotype assignment cluster

DNA extraction, amplification and sequencing

We dissected leg muscle from fresh, frozen, or alcohol‐preserved C. hookeri, and extracted genomic DNA using a salting‐out method (Sunnuck & Hales, 1996) or the Qiagen DNeasy tissue kit. We amplified and sequenced a mitochondrial genome fragment comprising the 3′ end of cytochrome oxidase I (COI), tRNA‐Leucine, and cytochrome oxidase II (COII), using a combination of primers C1‐J‐2195 and TK‐N‐3785, L2‐N‐3014 and TL2‐J‐3034 (Simon et al., 1994), and standard conditions described elsewhere (Morgan‐Richards & Trewick, 2005). We aligned and visually checked sequences using geneious v9 (Kearse et al., 2012), and excluded the tRNA‐Leu gene between COI and COII as this was missing from some specimens due to its use in PCR priming. We translated coding sequences to check for stop codons, frame shifts, and amino‐acid substitutions that might indicate nuclear copies. We included C. hookeri sequences from previous studies (Genbank: GU299870.1‐GU299966.1, AY940431.1, EU492994.1‐EU492973.1; Buckley et al., 2010; Morgan‐Richards & Trewick, 2005; Morgan‐Richards et al., 2010).

For 111 individuals across 15 population samples, we characterized nuclear DNA variation at 10 independent microsatellite loci (Myers et al., 2017). Multiplex PCR used thermal cycling conditions of 95°C for 5 min, followed by 28 cycles of 95°C for 30 s, annealing temperature of 60°C for 1 min 30 s, and 72°C for 30 s and one additional cycle of 60°C for 30 min. Reactions were resolved on an ABI Prism 3100 Genetic Analyzer (PE Applied Biosystems), with allele size determined using the Microsat plugin in Geneious with Genescan Liz‐500 (Applied Biosystems) as an internal size standard. To determine error rate we re‐genotyped six individuals at random and scored alleles blind.

This insect species has two large metacentric X chromosomes in females, and one X chromosome in males (Morgan‐Richards & Trewick, 2005). We inferred that one locus (Ch29‐14; Myers et al., 2017) was sex‐linked, as males always had a single allele. We confirmed that males inherited their mother's allele but not their father's by genotyping sons.

Population genetic analyses

Mitochondrial DNA

We examined mtDNA COI variation for 125 individuals (46 new sequences) across 41 population samples, and uploaded new unique haplotype sequences to GenBank (accession numbers MK532396, MK606153‐MK606171). The aligned data set is available on Dryad (https://doi.org/10.5061/dryad.b7t80m5) and at http://www.evolves.massey.ac.nz. For the full data set (1,389 bp), we inferred relationships among haplotypes using Maximum Likelihood with PHYML, implementing a GTR + I + G model of DNA evolution (Tavaré, 1986). We used subsets of the mtDNA haplotype data (1,307 bp) to infer median‐joining networks (Bandelt, Forster, & Rohl, 1999) using POPART (Leigh & Bryant, 2015).

Lineage age estimation

One mitochondrial lineage is associated with the all‐female populations that are geographically restrained to eastern and southern locations of New Zealand (Morgan‐Richards et al., 2010). The distribution of this lineage, may be associated with range and population expansion (Buckley et al., 2010). To infer the approximate time of the most recent common ancestor (MRCA) of this lineage, we implemented a molecular clock analysis in beast v2.5.0 (Bouckaert et al., 2014), using a 1,350 bp alignment of mtDNA sequence of COI and COII from 101 C. hookeri individuals. We evaluated nucleotide substitution models using jmodel test v0.1.1 (Posada, 2008) and found GTR + G to be a suitable model. No species specific rates of molecular evolution exists for C. hookeri, nor are there fossils for calibration, and using inferred timing of range expansion associated with climate cycling for calibration would introduce circularity to our analysis. Instead we used five different insect DNA substitution rates to calibrate our analyses. From the literature we obtained two rates based on insect mitochondrial interspecific divergence dates (Clarke, Levin, Kavanaugh, & Reimchen, 2001; Papadopoulou, Anastasiou, & Vogler, 2010), and three rates based on observed intraspecific nucleotide mutations (Gratton, Konopiński, & Sbordoni, 2008; Haag‐Liautard et al., 2008; Ney, Frederick, & Schul, 2018; see Table 1). A relaxed, uncorrelated, lognormal molecular clock model (Drummond, Ho, Phillips, & Rambaut, 2006) was applied, in order to evaluate the clock‐like behaviour of the data. The posterior probability of nonzero rate variance was close to zero (mean of rate variance ~0.0014), and therefore a simpler strict‐clock model was used in subsequent analyses. Because our data are from a subdivided species with a signature of population growth for some lineages (Buckley et al., 2010), we used two models for the tree prior; Coalescent constant population size, and Coalescent Bayesian Skyline (Drummond, Rambaut, Shapiro, & Pybus, 2005). An alternative model, Coalescent Extended Bayesian Skyline, assumes a single population but allows for any number of changes in population size; however, with this more parameter rich model, MCMC simulations failed to converge on a stationary distribution of posterior probabilities (e.g., ESS posterior = 166). Unless stated otherwise, default parameters were used. Markov Chain Monte Carlo (MCMC) simulations had chain lengths of 10 million, sampling every 1,000 generations. We assessed convergence through visual inspection of posterior statistics in tracer v1.5 (Rambaut & Drummond, 2007).

Table 1

The age of the most recent common ancestor (MRCA) of the parthenogenetic lineage of Clitarchus hookeri was estimated using five insect mitochondrial DNA mutation rates derived from different taxa and calibrations

Rate type	Taxa	mtDNA genes	Age of calibration (years)	Substitutions per site per million years	Reference	Age of parthenogen MRCA
						Constant pop size (million years)	SE of mean	Bayesian skyline (million years)	SE of mean
Interspecific	Darkling beetles	COI	10.5 million	0.0177	Papadopoulou et al. (2010)	0.1701	1.57E‐03	0.1537	1.55E‐03
Interspecific	Ground beetles	ND2, COI, COII, Cyt‐B	150,000	0.0285	Clarke et al. (2001)	0.1055	1.15E‐03	0.0964	1.01E‐03
Intraspecific	Fruit‐flies	~60% of mtDNA	200 generations	0.062a	Haag‐Liautard et al. (2008)	0.0493	4.44E‐04	0.044	4.89E‐04
Intraspecific	Katydids	COI	10,740	0.0792	Ney et al. (2018)	0.0381	3.85E‐04	0.0355	4.37E‐04
Intraspecific	Swallowtail Butterflies	COI	10–11,000	0.0960	Gratton et al. (2008)	0.0315	3.66E‐04	0.0281	3.72E‐04

The stick insect data set had 101 haplotypes and 1,350 bp of COI‐COII mtDNA sequence. Molecular clock analyses used a strict clock and either a constant population size or a Bayesian skyline model.

6.2 × 10–8 per site per fly generation (as study insects have ~1 generation per year we converted to 0.062 substitutions per site per million years).

Microsatellite loci

To estimate observed and expected heterozygosity per population sample (H O, H E) and average number of alleles per locus (N a) we used arlequin v3.5 (Excoffier & Lischer, 2010). We scored the sex‐linked locus as either homozygous for all males or as a single allele with missing data for all males. The results of these two coding methods did not alter population inferences from downstream analyses, so we present only the results from coding locus Ch29‐14 as sex‐linked. As these parthenogenetic stick insects could potentially reproduce via automixis or apomixis, we tested our population samples for evidence of deviations from Hardy‐Weinberg expectations (exact tests; autosomal loci only). For those samples that showed significant departures from Hardy‐Weinberg expectations, we tested for heterozygote deficiency using global Hardy‐Weinberg exact tests for multiple samples (U test), which are suitable for small samples and large numbers of alleles (Rousset & Raymond, 1995). We estimated pairwise F ST for each population sample and used an exact probability test of population differentiation (Raymond & Rousset, 1995) with genepop v4.2 (Rousset, 2008). To estimate population structure and assign individuals to clusters based on their genotypes at 10 nuclear loci, we used Bayesian model‐based analyses implemented in structure v2.3.4 (Falush, Stephens, & Pritchard, 2007; Pritchard, Stephens, & Donnelly, 2000). Again, we scored the sex‐linked locus as either homozygous for all males or as a single allele with missing data. As these yielded the same optimal K we present the results from coding locus Ch29‐14 as sex‐linked. For each structure run we conducted analyses with K ranging from one to seven, using an admixture model with correlated allele frequencies (Pritchard et al., 2000). Each run used 500,000 MCMC iterations following a 50,000 iteration burnin, based on recommendations provided by Gilbert et al. (2012). For each K set, we conducted ten replicate analyses. To identify the number of population clusters (K) with the best fit to our data, we used structure harverster (Earl & VonHoldt, 2012) to identify optimal K based on the posterior probability of the data for a given K, and the Delta‐K (Evanno, Regnaut, & Goudet, 2005).

RESULTS

Transition to asexual reproduction

Breeding experiments

Stick insect eggs were imported into New Zealand from the UK (n = 922). These eggs were the product of nine females (generation‐1) collected from the parthenogenetic population of Clitarchus hookeri on the Isles of Scilly. All nymphs that hatched (generation‐2) were female (n = 677). Ten adult females of generation‐2 laid eggs after mating with conspecific New Zealand males (n = 201). These eggs hatched into 82 daughters but no sons (generation‐3). A total of 13 adult females of this third generation laid 492 eggs after mating, which hatched 440 daughters and seven sons. Although we had data from two generations, the females were assumed to be genetically identical and were pooled. From all offspring (laid after mating) over both generations, we estimated that approximately 2.6% of nymphs were the product of sexual reproduction (number of sons × 2; to account for daughters also produced by sex; Table S2).

We genotyped 30 offspring from three generation‐3 mothers, selecting male and female nymphs who had hatched within a week of each other because sexual and parthenogenetic embryos develop at different rates (Morgan‐Richards et al., 2010). Of the 30 offspring, 21 were females with genotypes identical to their mother. Three daughters were not identical to their mother, suggesting they resulted from fusion of gametes (egg and sperm). These daughters were heterozygous at five or six loci (depending on their father's genotype) where their parents differed. All six genotyped sons were heterozygous at these same loci, with the exception of locus Ch29‐14 (Myers et al., 2017), where all males genotyped had only the maternal allele observed in the Isles of Scilly, UK population sample. Examination of our genotyping data for 142 individuals indicates that locus Ch29‐14 is sex‐linked. The genotype information confirmed that the males were produced via the fusion of sperm and egg (i.e., sex), rather than nondisjunction of one of their X chromosomes. Within our nonrandom sample of 30 nymphs hatched in the same week from three mothers, 21 (70%) were the result of asexual reproduction (Table S2).

Mitochondrial diversity

Mitochondrial DNA (cytochrome oxidase I & II) was sequenced and aligned with existing data (n = 125; DNA alignments available from Dryad: https://doi.org/10.5061/dryad.b7t80m5), representing 101 distinct C. hookeri haplotypes in the 1,389 bp alignment. All specimens of C. hookeri collected from the UK had the same haplotype. Haplotype genealogy (Figure 2) revealed that the UK haplotype joined a cluster of 17 haplotypes from specimens collected in Taranaki, New Zealand. Despite intensified sampling in this region, no New Zealand C. hookeri sampled had the UK haplotype that differed by three substitutions from the most similar Taranaki haplotype (0.23%; Figure 2d). Male specimens were collected from three locations contributing to the Taranaki haplotype cluster, where they were as common as females (n = 28 males, 23 females). Two Taranaki samples contained only females (Tarata n = 3; Rotorangi n = 6).

All stick insects collected from southern North Island, New Zealand had haplotypes that are part of the parthenogenetic mtDNA lineage previously identified (pink in Figure 2). This included all specimens sequenced from the sexual populations at Otaki (n = 5) and Wilton (n = 14), as expected from their location, but not expected from their reproductive mode (Table 2). We estimated the age of the most recent common ancestor of the mtDNA lineage associated with parthenogenetic reproduction by using rates of molecular evolution inferred for insects (Table 1). Our oldest estimate of about 170,000 years ago is based on an interspecific rate of nucleotide substitution (Papadopoulou et al., 2010) and is almost certainly too old due to the time dependency of the molecular clock (Ho, Phillips, Cooper, & Drummond, 2005; Molak & Ho, 2015). Estimates that suggest the most recent common ancestor of the parthenogenetic mtDNA lineage was alive about 38,000–28,000 years ago are based on intraspecific rates for the same mtDNA gene region, derived from a lineage sister to the Phasmids (Orthoptera) with similar generation times (Table 1).

Table 2

Higher genetic variation detected in sexual Clitarchus hookeri stick insect population samples compared to parthenogenetic population samples (number of alleles [N a], observed heterozygosity [H O], and expected heterozygosity [H E]) using 10 microsatellite loci, * denotes where genotype proportions deviated from expectations based on the Hardy‐Weinberg assumptions (using all nine autosomal loci with Fisher's exact test; p < .05)

Location	Region	Males common	mtDNA (n)	mtDNA parthenogenetic clade	Microsatellites					Inferred reproductive strategy
					n	Cluster (K = 3)	N a	H O	H E
Opanuku	Auckland	Yes	4	No	3	1	2.3	0.56	0.60	Sexual
Stony Bay	Coromandel	Yes	7	No	8	1	4.4	0.71	0.68	Sexual
Karapiro	Waikato	Yes	8	No	13	1	2.8	0.38	0.49	Sexual
East Cape	East Cape	No	5	No	3	1	1.7	0.07	0.60	Parthenogenetic
Urenui, Tarata, Rotorangi, Gordon	Taranaki	Yes	20	No	18	2	2.0	NA	NA	Sexual
Tresco, Isles of Scilly	UK	No	9	No	10	2	1.0	0	0	Parthenogenetic
Gisborne	Gisborne	No	5	Yes	6	3	1.3	0.06	0.38	Parthenogenetic
Turitea	Manawatu	No	3	Yes	8	3	1.2	0*	0.40	Automictic parthenogenetic
Otaki 2003	Wellington	Yes	3	Yes	5	2, 3	2.3	0.64	0.54	Sexual
Otaki 2013	Wellington	Yes	2	Yes	6	2, 3	2.8	0.52	0.54	Sexual
Otaki combined					11		3.1	0.54	0.53
Wilton 2003	Wellington	No	6	Yes	4	3	1.3	0.08	0.56	Parthenogenetic
Wilton 2016	Wellington	Yes	7	Yes	16	3	1.4	0.25	0.35	Sexual
Wilton combined					20		1.5	0.25	0.36
Manaroa	Marlborough Sounds	No	2	Yes	4	3	1.0	0	0	Parthenogenetic
Peel Forest	Canterbury	No	2	Yes	6	3	1.1	0	0.30	Parthenogenetic

The stick insects were sampled from 14 New Zealand and one UK location.

Population genetic structure

Ten polymorphic nuclear microsatellite loci were amplified and scored for 111 individuals collected from 15 locations (data available on Dryad: https://doi.org/10.5061/dryad.b7t80m5). Repeat amplification and scoring of individuals resolved the same genotypes, suggesting a low error rate (<0.016). Numbers of alleles per locus ranged from two to 11. Expected and observed heterozygosity levels were highest in sexual population samples (Stony Bay) and zero in some asexual population samples (e.g., Isles of Scilly UK; Manaroa NZ). The departure of genotype proportions from Hardy‐Weinberg expectations was significant (p < .05) in one population sample with a deficit of heterozygotes, as expected of an automictic parthenogen (Turitea; Tables 2 and S3). Observed heterozygosity was significantly higher in the six population samples with males (mean = 0.432) than the six all‐female populations (mean = 0.022; t test p = .0008). The sexual population sample from Wilton (2016) had fewer alleles and lower heterozygosity compared to sexual population samples from further north (Figure 3; only four polymorphic loci, no deviations from Hardy‐Weinberg expectations; chi‐squared 11.49; probability 0.17; Tables 2 and S3). The Otaki population samples (2003 and 2013) had eight polymorphic loci and genotype proportions met Hardy‐Weinberg expectations, as did all northern sexual population samples (Tables 2 and S3).

Our population samples were genetically differentiated from one another except for just seven (of 91) pairwise comparisons (exact G test; p < .05). The highest pairwise F ST estimate was between the UK and Manaroa samples (Table S4). Model‐based assignment of genotypes to clusters met an optimal fit to the genetic data with three population clusters (K = 3). Under this model, individuals collected from the same location had high probabilities of being assigned to the same cluster (Figure 3). Samples from northern New Zealand sexual populations grouped together (Opanuku, Stony Bay, Karapiro; Figure 3). This cluster of northern sexual populations also included the three individuals from East Cape, where males have not been detected and mtDNA haplotypes are sister to the parthenogenetic lineage. The UK Isles of Scilly specimens grouped with high assignment probability with the Taranaki samples, and the south‐eastern New Zealand populations formed the third group (Figure 3). All 11 individuals (2003 + 2013) from the sexual Otaki population had low assignment probabilities (0.56–0.69), due to the presence of alleles typical of the south‐eastern asexual lineage and Taranaki (except for two individuals who had alleles typical of the northern lineage; Figure 3). By contrast, all the individuals sampled from Wilton in 2016 clustered with the south‐eastern asexual lineage, with high assignment probability (0.71–0.99). The small sample (n = 4) from Wilton 2003 contained some allelic variation (three loci with two alleles each) but only a single individual (at a single locus) was heterozygous.

DISCUSSION

The capacity of Phasmids to reap the benefits of asexual reproduction are evident in the transition from sexual to parthenogenetic reproduction inferred for multiple independent lineages around the world (Scali, 2009; Schwander & Crespi, 2009). For the New Zealand species Clitarchus hookeri we estimated the most recent common ancestor of the widespread mtDNA lineage associated with parthenogenetic reproduction to have existed about 38,000–28,000 years ago. This coincides with the last glacial maximum (Rother et al., 2014) and is consistent with southward range expansion of this lineage when land connection allowed passage between North and South Islands (Trewick & Bland, 2012). The parthenogenetic mtDNA lineage is sister to another in East Cape, North Island (Morgan‐Richards et al., 2010), where a coastal refugia for this species during the LGM has been inferred (Buckley et al., 2010). The pattern of geographic parthenogenesis seen in C. hookeri (Morgan‐Richards et al., 2010) is the expected result of asexual reproduction being under positive selection when populations are growing during range expansion (Law & Crespi, 2002).

Origin of Isles of Scilly (UK) population

To infer the reproductive mode of the ancestral population from which the UK C. hookeri were derived, we examined genetic variation within natural populations in New Zealand. Evidence from both mitochondrial haplotypes and nuclear genotypes revealed the UK specimens were genetically most similar to Taranaki stick insects, a result that is consistent with plant collecting in New Zealand for the Tresco Abbey Gardens, Isles of Scilly, UK (Godley, 1997), providing an accidental source of stick insect eggs. Three of the five populations of C. hookeri sampled in Taranaki had an even sex ratio, and high heterozygosity, evidence that sexual reproduction is the predominant mode of reproduction in that region. In contrast, all C. hookeri specimens sampled in the Isles of Scilly (UK) were female and our sample contained neither mtDNA nor nuclear genetic variation, as expected of an obligate parthenogenetic population. It is possible that this population resulted from a single female or single egg transferred with the soil of a plant seedling from New Zealand, as parthenogenetic reproduction would have provided reproductive assurance in colonization (Baker's Rule: Baker, 1955).

Rapid acquisition of resistance to fertilisation by a parthenogen

In sexual populations of C. hookeri, mated females produce only sexual offspring, although each female is capable of parthenogenetic reproduction if males are absent (Morgan‐Richards et al., 2010; Salmon, 1955). Our C. hookeri data suggest that the Isles of Scilly (UK) population was derived from a sexually reproducing population (Taranaki, NZ), therefore, individuals from the UK population could be expected to return to 100% sexual reproduction if mated. However, we showed that captive females derived from the UK population had a barrier to fertilisation although no resistance to mating; their eggs and the offspring that hatched were mostly the result of parthenogenetic reproduction (~97%). Their barrier to fertilisation was similar to the ~5% sexually‐produced offspring previously documented for mated females of the geographic parthenogenetic lineage in New Zealand (Morgan‐Richards et al., 2010). Although apparently derived recently (~100 generations) from a Taranaki population, the UK lineage has developed resistance to fertilisation without apparent behavioural change. If there are costs associated with sperm storage or sexual receptivity we would expect selection to drive decay of sexual traits (Schwander, Crespi, Gries, & Gries, 2013). Therefore, the UK C. hookeri may have lost the propensity for sexual reproduction through selection involving sexual conflict.

Transition from parthenogenetic to sexual reproduction

Populations that naturally experience competition between individuals reproducing in different ways offer opportunities to determine what factors contribute to the local selective advantage (or disadvantage) of sexual reproduction. We used mtDNA and nuclear markers to determine the likely origin of two sexual C. hookeri populations within the range of the NZ parthenogenetic lineage. We expected that if these sexual populations were derived from long‐distance dispersal of sexual individuals (scenario 1: swamping), mtDNA and nuclear markers would show concordance in the placement of the samples within the genetic diversity of sexual C. hookeri. In contrast, the introduction of allopatric males (scenario 2: introgression), would result in retention of local maternal mtDNA and addition of alleles from a sexual population. This introgression scenario is what we observed in Otaki; males and females had mtDNA haplotypes that were part of the typical parthenogenetic lineage in this region, but their nuclear genomes contained alleles from two different genotypic clusters (sexual Taranaki and local parthenogenetic), resulting in a signal of mixture expressed as low assignment probabilities to clusters (Figure 3). Population samples separated by 10 years, (2003 and 2013) indicated introgression and establishment of sexual reproduction occurred prior to 2003, but high heterozygosity was retained in 2013.

If males had arisen via spontaneous loss of an X chromosome in situ (scenario 3: male genesis) only local nuclear alleles and local mtDNA haplotypes would be expected in their descendants. At Wilton we observed just this pattern; the parthenogenetic mtDNA haplotype and all nuclear alleles were those expected from the region. While Wilton samples from 2003 when males were rare and 2016 when males were common shared the same alleles at all loci, the 2016 sample had higher heterozygosity (0.08 compared to 0.025). These data suggest the most likely origin of males at Wilton was through in situ loss of an X chromosome, although we cannot exclude introgression and coincidental loss of all novel alleles. Given our small 2003 sample it is possible that we would miss alleles from outside the region if invasion involved just a few males. Over 13 years the sex ratio at Wilton changed from highly female skewed to an equal number of males and females, suggesting a rapid increase in successful sexual reproduction. Together, these observations suggest that the sexual populations at Wilton and Otaki constitute two independent and recent transitions to sexual reproduction from asexual reproduction in C. hookeri.

Local advantages of sexual reproduction

The best opportunities to understand the advantages of sexual reproduction over asexual reproduction come from studying natural replicates of the competition that plays out between sexual and asexual lineages (Neiman & Schwander, 2011). A local advantage for a particular reproductive strategy can be inferred when one strategy switches to another. The all‐female population in the UK was apparently derived from a population in the Taranaki region of New Zealand and was therefore likely to have been sexual (or recently derived from a sexual population). Accordingly, we infer a recent change from sexual to asexual reproduction. Intriguingly, our experimental work suggests this recent switch to asexual reproduction has been accompanied by a rapid loss of sexual propensity in the form of a barrier to fertilisation. An alternative scenario to consider is that the C. hookeri population was historically asexual at the time of transfer to the UK. In this instance the UK population would be derived from an asexual Taranaki population. This is plausible considering we have not sampled the exact UK mtDNA haplotype in New Zealand, however, the UK haplotype does nest within the Taranaki mtDNA diversity, suggesting either a recent switch to parthenogenesis in the UK lineage, or a switch to sexual reproduction in some Taranaki populations. Given our evidence from Otaki and Wilton males, it is indeed possible that Taranaki populations have switched from asexual to sexual reproduction, but the nuclear data give no indication of recent invasion. Despite the reproductive history of the original Taranaki population, finding the UK C. hookeri produced >90% of their offspring asexually after mating with conspecifics is significant considering the promiscuous nature of Clitarchus species (Langton‐Myers et al., 2019). Conflict between the sexes over egg fertilisation might influence the outcome of competition between reproductive strategies (Burke & Bonduriansky, 2017; Gerber & Kokko, 2016; Kawatsu, 2013). Models suggest that parthenogenetic populations are unlikely to establish in the presence of males if sexual conflict results in males coercing facultative parthenogens into sexual reproduction (Kawatsu, 2013, 2015).

Two separate New Zealand populations of C. hookeri provide evidence that sexual reproduction can replace parthenogenesis despite the numerical reproductive advantage provided by the latter. Heterozygote deficit within our parthenogenetic samples suggests automictic (rather than apomictic) reproduction. A switch to sexual reproduction from automictic parthenogenesis provides the benefits associated with outcrossing; increased heterozygosity and increased allelic diversity. Evidence for the switch in our study suggests these benefits can locally outcompete the numerical advantage of parthenogenetic reproduction.

At Otaki, the success of sex for C. hookeri might be linked to introgression from allopatric males that resulted in novel genotypes, increased allelic variation, and higher heterozygosity. By contrast, at Wilton, putative‐local males with local alleles reveal that rapid and successful switching to sexual reproduction cannot be attributed solely to increased allelic diversity provided by outsourced males. In this study system, we have spatial and temporal sampling to infer the outcome of competition among reproductive strategies. Continued observation will reveal whether switches from asexual to sexual reproduction are stable (Innes & Ginn, 2014) and whether males continue to expand their range. The value of temporal sampling for documenting the rapid spread of sex and identifying short‐term local effects in wild populations are clear. Our study suggests the short‐term advantages of outcrossing can enable sex to successfully outcompete automictic parthenogenetic reproduction. Females with lower resistance to fertilisation might gain indirect reproductive success through their sons inheriting an ability to overcome fertilisation barriers (Kawatsu, 2015). This new study system with multiple independent transitions in reproductive strategy presents a rare opportunity to focus on location‐specific forces that provide sexual reproduction with a competitive advantage over parthenogenetic reproduction.

AUTHOR CONTRIBUTIONS

M.M.R., and S.A.T. conceived the project, collected the insects and raised their offspring in captivity. M.M.R. received funding and generated genotypes and haplotype sequences. M.M.R., S.A.T., and S.S.L.‐M. analysed data. All authors contributed to design and writing the manuscript.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.