Unravelling the evolutionary history of Eucalyptus cordata (Myrtaceae) using molecular markers

Additional keywords: Diversity Array Technology (DArT), Eucalyptus, gene flow, genetic drift, heteroblastic, homoblastic, microsatellite, spatial genetic structuring.

Eucalyptus cordata Labill. (heart-leaved silver gum) is a rare tree species of the subgenus Symphyomyrtus, endemic to the island of Tasmania. The species is notable compared with other Tasmanian endemics of the series Orbiculares (E. archeri, E. gunnii, E. morrisbyi, and E. urnigera), being homoblastic and retaining the striking, glaucous, opposite and sessile juvenile leaves into reproductive maturity in both wild and ornamental populations, even as a large tree (Brooker 2000; Nicolle 2006; Slee et al. 2006; Fig. 1). Several populations are known where adult foliage does develop on the tops of trees, but flowering still occurs in the distinctive juvenile leaf stage (Potts 1989). Such neotenic homoblasty is relatively rare in eucalypts, but occurs in eucalypts from diverse taxonomic groups (Potts and Wiltshire 1997) as well as in two other species of the series Orbiculares (E. perriniana and E. pulverulenta; Brooker 2000; Nicolle 2006; Slee et al. 2006).

Fig. 1.  (a) The tallest known native Eucalyptus cordata at Moogara (Population 18, Table 1) (indicated by arrow), where homoblastic trees bearing only the glaucous, juvenile-type foliage characteristic of the species may reach over 20 m tall and grow with E. regnans. (b) The poorly drained habitat of E. cordata at Brown Mountain (Population 1); arrow points to E. cordata growing with the sedge Gahnia gradnis. (c) The coastal habitat of E. cordata subspecies cordata at Cape Queen Elizabeth (Population 4), on the western side of Storm Bay; white arrows point to E. cordata growing on the southern edge of Cape Queen Elizabeth, and the black arrow shows the type locality of E. cordata at Penguin Island (Population 10) in the distance. (d) Eucalyptus cordata subspecies cordata on the summit of Perpendicular Mountain (Population 11), where it can grow in close proximity to E. globulus (foliage above the glaucous E. cordata); arrow shows the location of the nearest E. cordata population at Hellfire Bluff (Population 6) in the distance. (e) The mallee form of E. cordata at Square Mountain (Population 13) (white arrows) growing under the canopy of surrounding E. obliqua. Photos by PAH.

Eucalyptus cordata occurs in small, isolated populations across a remarkably diverse range of habitats in south-eastern Tasmania. It grows on dry cliff-tops or hill slopes at lower altitudes and on poorly drained marsh edges at higher altitudes, and has a corresponding diversity of habits from stunted, multi-stemmed shrub to tall, single-stemmed forest trees more than 20 m high (Potts 1989; Williams and Potts 1996). A morphological cline has been described across this range of habitats (Potts 1989), and has been summarised by the recognition of two subspecies (Nicolle et al. 2008). The main discriminating characters between these two subspecies are stem shape and leaf size, which are square (as opposed to round) and larger, respectively, in subspecies quadrangulosa, and there is a genetic basis to this as shown in a glasshouse seedling trial (Potts 1989). Subspecies quadrangulosa occurs at higher altitude and wetter sites in the south-west of the species range (Nicolle et al. 2008).

Eucalyptus cordata has been a key species in studies of chloroplast sharing and reticulate evolution in eucalypts (McKinnon et al. 2004a). Extensive sharing of chloroplast haplotypes of the ‘Southern’ type has been reported among all Symphyomyrtus species occurring in south-eastern Tasmania, and is believed to be the result of hybridisation during species range contractions or expansions in glacial refugia (McKinnon et al. 2001a; McKinnon et al. 2004a). Eucalyptus cordata is one of the only species that are fixed for this chloroplast clade (McKinnon et al. 2004b), and detailed population studies have provided evidence that the chloroplast and rare components of the nuclear genome have been locally captured by the more widespread, co-occurring E. globulus (McKinnon et al. 2004b, 2010).

The only endemic Symphyomyrtus species with contemporary natural distributions confined to this south-eastern forest refugium are E. cordata and the endangered E. morrisbyi. Traditionally, studies have shown that populations within glacial refugia tend to be more genetically diverse than populations that have expanded into more recently available habitat (Comps et al. 2001). However, this pattern may well be case-specific, and studies on eucalypts have suggested that the restricted and fragmented distribution in putative glacial refugia can result in genetically depauperate populations because of genetic erosion (Byrne and Hopper 2008). This is indeed the case with E. morrisbyi (Jones et al. 2005), and it suggests that the disjunct populations of E. cordata may also demonstrate low genetic diversity and high levels of differentiation due to their isolation.

The present study used Diversity Array Technology (DArT), nuclear microsatellites, and chloroplast markers to investigate the population genetic processes that have shaped the genetic architecture of E. cordata. The microsatellite and chloroplast markers were chosen so that the results of this study could be compared with previous works (McKinnon et al. 2004b; Jones et al. 2005, 2013; Rathbone et al. 2007). We used DArT markers, which have wide genome coverage (Sansaloni et al. 2011; Steane et al. 2011), so as to achieve high resolution for the phylogenetic component of the study. The goals of this study were to determine (1) the genetic affinities of E. cordata within series Orbiculares; (2) the genetic affinities of populations within E. cordata; and (3) the evolutionary origin of subspecies quadrangulosa. We present evidence that multiple expansion and contraction events through evolutionary time have shaped the genetic architecture of this species, and that the newly described subspecies quadrangulosa is a neo-endemic.

To investigate the genetic affinities between populations of E. cordata, total genomic DNA was extracted from fresh leaf tissue collected from 533 individuals from 21 natural populations representing the two subspecies (subspecies quadrangulosa and subspecies cordata) and their intermediates (Table 1, Fig. 2). Additionally, to investigate the genetic affinities between populations of E. cordata and E. morrisbyi, DNA was extracted from fresh leaves collected from 24 individuals from the Risdon Hill population of E. morrisbyi (Fig. 2), and DNA samples that were available from a previous study on the Calverts Hill E. morrisbyi population (Jones et al. 2005) were re-genotyped. Trees were sampled at least two canopy heights apart to avoid sampling clonal and closely related individuals (Jones et al. 2005). All trees were photographed and tagged, with their geographic position recorded. Contemporary population sizes were estimated based on the number of individuals (either clumps of stems or individual stems) observed during sampling. Phenotypic affinities of E. cordata populations were classified according to Nicolle et al. (2008).

Table 1.  Population code, altitude, population sizes, and genetic diversity parameters for the 23 analysed Eucalyptus cordata and E. morrisbyi populations
Code, population identification code; N, estimated population size; n, number of individuals sampled; H_E and H_O, expected and observed heterozygosity; F, Wright’s fixation index; SSR A_R and CP A_R, microsatellite and chloroplast allelic richness. Allelic richness was rarefied to 17 individuals for the microsatellite data and four individuals for the chloroplast data

Fig. 2.  The distribution of E. cordata (closed symbols) and E. morrisbyi (+) in south-eastern Tasmania. Populations studied are numbered, with the colour of the number (blue, Eucalyptus cordata subspecies cordata; purple, E. cordata intermediate; black, E. cordata subspecies quadrangulosa; red, E. morrisbyi) and shape of symbol corresponding to morphological affinities of the population (see Table 1 for population details).

Total genomic DNA was extracted from 10 mg of ground tissue using a CTAB protocol (Doyle and Doyle 1990) with the modifications of McKinnon et al. (2004b), grinding fresh tissue in liquid nitrogen with either a mortar and pestle or a tungsten carbide bead in a 1.2-mL sample tube placed in a Mixer Mill MM400 (RETSCH, Haan, Germany) for two 1-min periods at 30 Hz. DNA quality and quantity were visually estimated by comparing the Lambda HindIII molecular weight standard (New England BioLabs, Ipswich, MA, USA) against the electrophoretically separated DNA.

The genetic affinities of E. cordata to other species from the series Orbiculares was investigated using genome-wide DArT markers (Sansaloni et al. 2011). This included a subset of E. cordata samples from the current study (including both subspecies and their intermediates), along with representative samples from Tasmanian endemic species (E. gunnii, E. urnigera, E. morrisbyi, E. archeri) and species from mainland Australia (E. saxatilis, E. glaucescens, E. perriniana, E. pulverulenta) of the series Orbiculares (Appendix 1). DNA was prepared following the protocol of Steane et al. (2011). Genome complexity reduction was conducted following Sansaloni et al. (2011), with polymorphisms detected by DArT Pty Ltd (http://www.diversityarrays.com/, verified 3 April 2014) based on the presence or absence of the DArT marker in the sample. The dissimilarity between samples was calculated using an additive Dollo distance (ADD), a distance measure appropriate for the DArT data (Woodhams et al. 2013). The matrix of ADD distances was summarised using a principal coordinate analysis undertaken with GenAlEx version 6.501 (Peakall and Smouse 2012).

All E. cordata and E. morrisbyi individuals were genotyped using 13 microsatellite loci. Primer pairs were designed in E. grandis (prefix EMBRA; Brondani et al. 1998; Brondani et al. 2006) and E. globulus (prefix EMCRC; Steane et al. 2001) and were separated into three multiplex mixes (Appendix 2). Microsatellite loci were amplified by polymerase chain reaction (PCR) using a QIAGEN Multiplex PCR kit (QIAGEN, Melbourne, Vic., Australia). Each reaction contained 1 × QIAGEN Multiplex PCR Master Mix, 1 × Q-solution, 5–15 ng DNA template, primers (Appendix 2), and RNase-free water to a final volume of 5 μL. Amplification of microsatellite fragments was performed on a Bio-Rad C1000 Thermal Cycler (Bio-Rad Laboratories, Hercules, CA, USA) using two PCR programs. Multiplex mixes A and C used the following program: 95°C for 15 min; 30 cycles of 94°C for 30 s, annealing temperature of 57°C for 90 s, and 72°C for 60 s; a final extension at 60°C for 10 min; finishing with a final hold temperature of 15°C. Amplification of Multiplex mix B used an annealing temperature of 55°C for 90 s. Final PCR product concentration and quality was estimated by comparing a 100-bp molecular weight standard (New England BioLabs) with the electrophoresis separation of the PCR product on a 1% agarose gel pre-stained with Goldview (Guangzhou Geneshun Biotech, Guangzhou, China). PCR products were diluted at a ratio of 1 μL of PCR product to 10 μL of Milli-Q water. One μL of diluted PCR product was dried and sent to the Australian Genome Research Facility (AGRF), South Australia. There, it was resuspended and sized using an AB3730 DNA Analyser (Applied Biosystems, Foster City, CA, USA). Allele sizes were assessed in GeneMapper version 3.7 (Applied Biosystems) using GS500 k(–250)LIZ as a size standard.

To estimate the frequency of null alleles, we used the program INEst (Chybicki and Burczyk 2009). This detected null alleles in EMCRC2 and EMCRC7. However, the removal of these loci did not affect the genetic diversity parameters, so all loci were included. The microsatellite loci were powerful enough to distinguish related and clonal genotypes in E. morrisbyi (Jones et al. 2005) and E. cordata (data not presented). Where clonal genotypes were detected, one representative genotype was retained from each ‘clonal patch’ for further analysis. This resulted in the removal of a sample from Chimney Pot Hill that was 3.6 m from its adjacent clonal sample. It also allowed validation of a large clone ~2 m across on Penguin Island where multiple stems had been specifically sampled from a large lignotuber remnant. Only individuals with allele scores for at least eight loci were retained, resulting in a sample size of 567 individuals from E. cordata and E. morrisbyi. Repeatability of microsatellite sizing was assessed by regenotyping 10% of the E. cordata samples, which showed that the genotyping error rate was only 1.8%.

Genetic affinities between individuals from populations of E. cordata and E. morrisbyi were investigated using the Bayesian methods implemented in STRUCTURE version 2.3.3 (Pritchard et al. 2000; Falush et al. 2003). An admixture model with no a priori population grouping was run, assuming that: (1) there are K populations characterised by allele frequencies that are correlated among populations; (2) individuals have inherited a proportion of their ancestry from each population; (3) there is linkage equilibrium; and (4) populations are in Hardy–Weinberg equilibrium. All STRUCTURE analyses used a burn-in period of 250 000 iterations followed by 100 000 Markov chain Monte Carlo iterations. Twenty independent realisations at each value of K (K = 2 to K = 23) were performed. STRUCTURE Harvester (Earl and vonHoldt 2012) was used to detect the optimum value of K using the post hoc methods of Evanno et al. (2005). Solutions for K were then collapsed into a single representative solution using the Greedy algorithm implemented in CLUMPP version 1.1.2 (Jakobsson and Rosenberg 2007), weighted by the number of individuals in each population.

Population genetic diversity parameters were calculated and averaged across populations and loci by using Genetic Data Analysis version 1.1 (Lewis and Zaykin 2001). This included the observed number of alleles (A), observed (H_O) and expected (H_E) heterozygosity, and Wright’s fixation index (F). Allelic richness (A_R), rarefied to a minimum sample size of 17 individuals per population, and Nei’s (1987) estimate of total heterozygosity (H_T) were calculated in FSTAT version 2.9.3.2 (Goudet 2002). Correlations among the genetic diversity measures (H_E, H_O, A_R, and F), altitude, and population size were tested using a Pearson’s product–moment test in the R statistical package (R Core Team 2012). To determine the relationships among populations, Nei’s (1972) genetic distance matrix among all populations was calculated in Genetic Data Analysis. The F-statistics F_IS (inbreeding within individuals relative to their population), F_IT (inbreeding within individuals relative to all other populations), and F_ST (inbreeding in populations relative to all other populations) were estimated in FSTAT with 99% bootstrap confidence intervals derived after 1000 replications across loci. A pair-wise F_ST matrix was also produced using FSTAT, with statistically differentiated populations indicated at the α = 0.05 level after 4200 permutations. The level of differentiation among subspecies cordata, subspecies quadrangulosa and intermediate populations of E. cordata was investigated by pooling populations based on their taxonomic classification and calculating the F_ST between these three groups in FSTAT. The Mantel test in GenAlEx (version 6.501) was used to determine whether the level of genetic differentiation, based on the matrix of F_ST/(1 – F_ST) values (Rousset 1997), was related to the matrix of geographic distances (natural logarithm transformed) among populations.

The chloroplast genome in Eucalyptus has been shown to be maternally inherited, meaning it can only be dispersed by seed (Byrne et al. 1993; McKinnon et al. 2001b), which, in Eucalyptus, is usually gravity dispersed (Cremer 1977). The hypervariable J_LA+ (~630 bp) region of the chloroplast genome was sequenced in 4–7 geographically well-spaced trees for each of the E. cordata and E. morrisbyi populations studied (n = 127); the forward primer rpl2 and reverse primer eucpsbA were used (Freeman et al. 2001). Chloroplast sequences from previous studies of E. cordata (McKinnon et al. 2004b) and E. morrisbyi (Jones et al. 2005) were used where possible, with additional samples sequenced to increase the total sample size to a minimum of five per population. For each DNA sample, ~5–15 ng DNA was used in a final volume of 25 μL PCR reactant mix as in Freeman et al. (2001), except that 2.5 mM MgCl₂, 0.15 μM of each primer, 1 unit MangoTaq polymerase (Bioline Australia, Alexandria, NSW), and deionised Milli-Q water to 25 μL was used. Amplification was performed on a Bio-Rad C1000 Thermal Cycler (Bio-Rad Laboratories) using the following PCR cycling program: 95°C for 5 min; 30 cycles consisting of 94°C for 60 s, 61°C for 60 s, and 72°C for 60 s; 72°C for 5 min; finishing with a final hold temperature of 15°C. PCR product concentration and quality was estimated by electrophoresis separation on a 2% agarose gel pre-stained with Goldview (Guangzhou Geneshun Biotech), using a 100-bp molecular weight standard (New England BioLabs). Purification and sequencing of PCR products were performed by Macrogen (Seoul, South Korea) using an AB3730 DNA Analyser (Applied Biosystems).

Forward and reverse sequences were analysed in Sequencher version 4.7 (Gene Codes, Ann Arbor, MI, USA) and CLC Workbench version 4 (CLC Bio, Aarhus, Denmark). The chloroplast haplotypes were then compared with other cpDNA sequences of Symphyomyrtus species held in the School of Biological Sciences DNA database at the University of Tasmania. New chloroplast haplotype sequences were classified and named following the nomenclature of Freeman et al. (2001). Chloroplast haplotype diversity (allelic richness rarefied to four individuals per population, CP A_R) and F_ST were estimated in FSTAT. The Pearson’s correlation of CP A_R with the nuclear genetic parameters, altitude, and population size was calculated as described above.

To determine the relative historic rates of pollen and seed migration between populations of E. cordata and E. morrisbyi, the ratio of pollen to seed gene flow was estimated using eqn 5a from Ennos (1994). This equation uses the level of differentiation among populations (F_ST) in the nuclear and maternally inherited chloroplast markers and takes into account the level of inbreeding by using the F_IS derived from the nuclear microsatellite data.

The genetic relationship of E. cordata to other species from series Orbiculares was determined using a principal coordinate analysis based on the ADD dissimilarity matrix derived from 3172 non-redundant polymorphic DArT markers assayed in 49 representative samples of 10 species (Fig. 3). The first axis separated the mainland samples of the series Orbiculares from the Tasmanian endemics, with one E. gunnii sample (Sample 31) intermediate between these two groups (Fig. 3). Eucalyptus cordata (Samples 5–21, Fig. 3) showed closer genetic affinities with the Tasmanian endemic species from the series Orbiculares than with the homoblastic species from the mainland (E. pulverulenta, Sample 38, and E. perriniana, Sample 39 in Fig. 3). Axis 2 in the principal coordinate analysis showed a cline in genetic affinities within the Tasmanian endemic group (Fig. 3). At one end of this molecular cline were the majority of E. cordata samples, which were well differentiated from the other Tasmanian endemics. The other Tasmanian endemics were at the other end of this cline and were poorly differentiated from each other in this two-dimensional space, and remained poorly differentiated in the third dimension (data not shown). The E. cordata samples from Brown Mountain (Samples 8 and 18 in Fig. 3), a population that transitions into petiolate adult leaves while still flowering in the juvenile phase, grouped with the other Tasmanian endemics, as did the intermediate samples at Knights Creek and Moogara (Samples 7 and 11 in Fig. 3). However, samples from the heteroblastic E. morrisbyi populations at Calverts Hill (Samples 35 and 36 in Fig. 3) and Risdon Hill (Samples 33 and 34 in Fig. 3) had close affinities to E. cordata and were embedded within the principal coordinate space occupied by the E. cordata samples, which suggested that E. morrisbyi should be included in further analyses of E. cordata.

Fig. 3.  Principal coordinates analysis of the additive Dollo distance (ADD ) matrix derived from 3172 genome-wide Diversity Array Technology (DArT) markers assayed in 49 individuals representing nine species from the series Orbiculares. The number colour corresponds to each of the species included in this study: blue, Eucalyptus cordata; red, E. morrisbyi; green, E. gunnii; orange, E. archeri; purple, E. urnigera; black, mainland species (see Appendix 1 for population details).

The 13 microsatellite loci used to genotype individuals from E. cordata and E. morrisbyi were highly polymorphic, ranging from 18 to 31 alleles detected per locus (mean A across loci = 24.4, data not shown). Total genetic diversity (H_T) averaged across the 13 loci was slightly higher in E. cordata (H_T = 0.86) than in E. morrisbyi (H_T = 0.80), with many more alleles per loci, on average, detected in E. cordata (A = 24) than in E. morrisbyi (A = 10). The H_E and H_O within the E. cordata populations averaged 0.79 and 0.74, respectively (Table 1). There was higher genetic diversity in populations of subspecies quadrangulosa than populations of subspecies cordata, as measured by the H_E and H_O as well as allelic richness (A_R) (Table 1). The small E. cordata population on Penguin Island (Population 10 in Fig. 2) and both E. morrisbyi populations (Populations 22, 23) had the lowest genetic diversity (Table 1). Spatial plotting of population genetic diversity revealed two regions of elevated allelic richness, one in the north-east and the other in the south-west of the range of E. cordata (Fig. 4a). Similar patterns were evident in the H_E and H_O of populations (Table 1). Populations peripheral to these regions of elevated diversity had lower genetic diversity (A_R < 7.6), but there were exceptions such as Brown Mountain (Population 1, A_R = 8.7, Fig. 4a). Population size could not explain this spatial pattern, with the Pearson’s product–moment correlation test showing no linear relationship between the number of individuals in each population and any of the genetic diversity parameters (Table 2). However, there was a tendency for H_O to be significantly lower and the inbreeding levels (F) significantly higher in lower altitude populations of E. cordata (Table 2).

Fig. 4.  (a) Map of nuclear microsatellite allelic richness (A_R) and (b) chloroplast A_R for populations used in this study. Increasing size of the triangles corresponds to increasing levels of population A_R above the species mean A_R, and an increase in the size of the circles correspond to a decrease in population A_R below the species mean A_R. Population codes are detailed in Table 1.

Table 2.  Estimated Pearson’s product–moment correlation matrix among the population genetic parameters, altitude, and the logarithmic number of individuals for the 21 Eucalyptus cordata populations
Significant correlations (P < 0.05) are shown in bold. log(N), logarithm of population size (Table 1); H_E and H_O, microsatellite expected and observed heterozygosity; F, Wrights fixation index; SSR A_R and CP A_R, microsatellite and chloroplast allelic richness. Significant (P < 0.05) correlations among variables when the two E. morrisbyi populations were included in the analysis are shown in bold

The populations of E. cordata were well differentiated, with an average F_ST of 0.09 (range 0.02–0.20; Appendix 3). This compares with an F_ST of 0.14 between the two E. morrisbyi populations and an average F_ST between the E. cordata and E. morrisbyi populations of 0.13. Populations of subspecies quadrangulosa were less differentiated (F_ST = 0.04, n = 3) than the E. cordata intermediates (F_ST = 0.09, n = 4) and subspecies cordata (F_ST = 0.10, n = 14). Populations of these subspecies were equally differentiated, on average, from the E. morrisbyi populations, but the smallest pair-wise values of F_ST involved the subspecies cordata populations from Corbetts Hill (Population 5) and Square Mountain (Population 13). All estimates of pair-wise F_ST between population pairs were statistically significant, even the lowest pair-wise F_ST (0.02), which was between the geographically adjacent populations of subspecies quadrangulosa at Snug Plains (21) and Sherwood Hill (20). Most of the higher pair-wise F_ST estimates involved the small E. cordata population at Penguin Island (10) (average pair-wise F_ST from other E. cordata populations, 0.16; Appendix 3). Penguin Island is the type locality of E. cordata (Potts 1988) but it is clearly a molecular outlier. An isolation-by-distance model could not explain the pattern of genetic differentiation between populations, with a Mantel test showing no relationship between geographic and genetic distance (R_xy = 0.18, P = 0.08, n = 21).

The genetic affinities of the E. cordata and E. morrisbyi individuals were analysed using STRUCTURE with no a priori population assignment. The post hoc ΔK method of Evanno et al. (2005) showed that the assignment of individuals into five genetic clusters (K = 5) had the highest ΔK value, suggesting this as the most optimal summary of their genetic affinities (Appendix 4). Although other K values had reasonably high ΔK estimates, additional genetic clusters tended only to separate individuals from specific populations and were not as meaningful as K = 5. Most individuals were assigned to two geographically separated clusters and comprised individuals from populations in the north-east (13, 5, 6, 11, 12, and 2) and south-west (17, 19, 16, 20, 21, 15, 9, and 7) of the E. cordata range (Fig. 5). The south-western cluster was more phenotypically diverse than the pure subspecies cordata cluster in the north-east and included individuals from subspecies cordata, subspecies quadrangulosa and their intermediates. The south-western cluster also shared partial genetic affinities with a smaller, third cluster of individuals from subspecies cordata populations at Taranna (14), Chimney Pot Hill (3), and Moogara (18) (Fig. 5). Two smaller clusters reflected affinities between E. cordata and E. morrisbyi, with individuals from the E. cordata populations at Penguin Island (10) and Brown Mountain (1) clustering with the E. morrisbyi population at Calverts Hill (23), whereas individuals from the E. cordata population at Meehan Range (8) clustered with the nearby E. morrisbyi population at Risdon Hill (22) (Fig. 5).

Fig. 5.  The average proportion of membership for sampled individuals from each of the Eucalyptus cordata and E. morrisbyi populations into the K = 5 STRUCTURE groups. The colour of the population number corresponds to each species: E. cordata subspecies cordata (blue), E. cordata intermediate (purple), E. cordata subspecies quadrangulosa (black), and E. morrisbyi (red). Population codes are detailed in Table 1.

In total, 27 chloroplast haplotypes were identified following sequencing of the J_LA+ chloroplast region in 127 samples from 23 populations of E. cordata and E. morrisbyi (Appendix 5), and these all showed the Southern haplotype (suffix JS; Freeman et al. 2001; McKinnon et al. 2004b). Twenty-four of these chloroplast haplotypes were private to a single population. Only JS81 was shared between E. cordata and E. morrisbyi populations and this sharing involved only O’Briens Hill (9) and the E. morrisbyi population at Risdon Hill (22), which are 30 km apart. Eleven of the 21 E. cordata populations and both E. morrisbyi populations were fixed for a single chloroplast haplotype (Appendix 5), and in four cases, this involved fixation for the most common haplotype (JS43). Nevertheless, higher diversity of chloroplast haplotype was found in E. cordata (H_T = 0.89) than in E. morrisbyi (H_T = 0.50). The chloroplast haplotype (allelic) richness calculated from this E. cordata dataset was not significantly correlated with the microsatellite allelic richness parameters, altitude or population size (Table 2). The majority of populations had low diversity of the chloroplast genome, except for populations at Square Mountain (13), Sherwood Hill (20) and Snug Plains (21) (Table 1, Fig. 4b). The E. cordata populations were highly differentiated in their chloroplast haplotypes, as indicated by an F_ST of 0.70. This high level of chloroplast differentiation among populations compared with the differentiation in the nuclear SSRs resulted in an estimation of 22 times more pollen- than seed-mediated gene flow in E. cordata.

Our genome-wide genotyping of multiple species from the series Orbiculares showed that E. cordata is more aligned to the Tasmanian endemic species of series Orbiculares than to mainland species, and suggests E. cordata most likely evolved from a common ancestor shared with the Tasmanian endemics. This result confirms findings based on amplified fragment length polymorphism (AFLP) markers that three samples of E. cordata formed a monophyletic group within the Tasmanian Orbiculares species (McKinnon et al. 2008). Similarly, the phylogenetic study of Steane et al. (2011) found that the homoblastic E. cordata was more aligned to the heteroblastic Tasmanian Orbiculares species than to other homoblastic species of the Orbiculares, such as E. perriniana and the mainland E. pulverulenta. Indeed, these molecular studies suggest that these two homoblastic species from the mainland are more aligned to series Viminales than to Orbiculares (McKinnon et al. 2008; Steane et al. 2011). This also implies that neotenic homoblasty in E. cordata has evolved independently from a heteroblastic ancestor, and is consistent with the hypothesis that the retention of the juvenile foliage into reproductive maturity has evolved independently on multiple occasions within Eucalyptus (Potts and Wiltshire 1997). Although most E. cordata populations are homoblastic, trees in the Brown Mountain (1) and Penguin Island (10) populations, still flowering in the juvenile leaf stage, have been observed with adult leaves in the wild (Potts 1988). The heteroblastic characteristic of these populations, coupled with their genetic affinities to individuals from the larger, Calverts Hill (23) population of E. morrisbyi, suggests that these peripheral E. cordata populations may be relicts of the ancestral lineage from which the contemporary E. cordata evolved. This would explain the genetic affinities we detected between the southern and northern extremities of the E. cordata distribution (Link 1, Fig. 6), despite the small population at Penguin Island (10) being a genetic outlier. Penguin Island is the type locality for E. cordata (Potts 1989), and is clearly not typical of the species from a genetic (present study) and phenotypic (Potts 1988) perspective.

Fig. 6.  Hypothesised glacial expansion routes that have structured the genetic diversity in Eucalyptus cordata and E. morrisbyi populations. Line colour refers to the molecular affinities of populations based on the K = 5 STRUCTURE groups defined in Fig. 5. The lines are numbered according to a hypothesised chronology of evolutionary events based on the extant molecular variation (see Discussion). Line 1 reflects the ancestral links between populations of E. cordata and E. morrisbyi; line 2 reflects the ancient gene flow between E. cordata and E. morrisbyi; line 3 reflects the east–west expansion across Storm Bay and the River Derwent; line 4 reflects the eastern expansion of E. cordata populations in the north-eastern lineage (dashed line shows extent of north-eastern core); and line 5 reflects the expansion of the south-west lineage (dashed line shows extent of south-western core). Black dashed line shows the hypothesised extent of landmass during the Last Glacial Maximum, based on a ~120-m depression of sea levels (Jackson 2005).

Although E. cordata is locally distributed as a series of isolated populations, the nuclear microsatellite diversity within populations as measured by H_E is similar to that reported for widespread eucalypt species with more continuous distributions (average H_E of three species 0.79; Byrne 2008). For example, the H_E (0.79) and H_O (0.74) in E. cordata populations are similar to those observed in Tasmanian populations of E. obliqua (H_E = 0.80, H_O = 0.79; Bloomfield et al. 2011), E. pauciflora (H_E = 0.80, H_O = 0.75; A. Gauli, unpubl. data), and E. globulus (H_E = 0.81, H_O = 0.70; Jones et al. 2013). By contrast, the two other eucalypt species studied that also occur in Tasmania as small, isolated populations have lower H_E and H_O than the average observed for the E. cordata populations: E. perriniana (H_E = 0.72, H_O = 0.62; Rathbone et al. 2007) and E. morrisbyi (H_E = 0.69, H_O = 0.64; Jones et al. 2005). These observations suggest that, on average, most of the E. cordata populations have not been isolated and/or have persisted as small populations for a long time from an evolutionary perspective. This conclusion is also supported by the moderate level of population differentiation (F_ST = 0.09), which is below the average reported for regionally (average F_ST of two species, 0.25) and locally (average F_ST of two species, 0.15) distributed eucalypt species, and more similar to the mean for widespread species (average F_ST of two species, 0.06) (Byrne 2008). Loss of heterozygosity for neutral molecular markers requires isolation and persistence of small population sizes for multiple generations (Ellstrand and Elam 1993), and none of these Tasmanian species comprising isolated populations, including E. cordata, have the low population genetic diversity and high population differentiation evident in species such as E. curtisii (F_ST = 0.30; Smith et al. 2003) and E. caesia (F_ST = 0.53; Byrne and Hopper 2008), which have persisted as isolated populations in ancient landscapes (Byrne 2008). These species have smaller population sizes than E. cordata, and they have no doubt experienced prolonged bottlenecks through their longer evolutionary history.

A positive relationship between genetic diversity parameters (A_R, H_E and H_O) and population size may be expected due to genetic drift and inbreeding in small populations (Ellstrand and Elam 1993; Leimu et al. 2006). However, the absence of such a relationship is common in eucalypts (reviewed by Butcher et al. 2005), and this was the case for E. cordata. Allelic richness is a highly sensitive genetic diversity parameter and is influenced by the loss of rare alleles following bottlenecks in population size (Maruyama and Fuerst 1985). There is strong evidence that the small E. cordata population at Penguin Island has experienced a persistent bottleneck, based on its low allelic richness (A_R = 5.4, compared with the species mean of 8). However, little support was found for major bottlenecks affecting genetic diversity in other small populations, suggesting that these populations have only recently become small and isolated, at least in terms of tree generations.

Eucalypts can have long generational times because of regeneration from underground lignotubers, which can fragment with time and result in disconnected, clonal patches (Tyson et al. 1998; Byrne 2008). In the case of E. morrisbyi, Jones et al. (2005) estimated that clonal patches in the Risdon Hill (22) population were at least 670–1523 years old. Many of the E. cordata populations (particularly of subspecies cordata) similarly comprise individuals that have regenerated vegetatively from large basal lignotubers (Potts 1989). Applying the rate of radial growth of lignotubers reported by Tyson et al. (1998) to the isolated clonal stems detected at Chimney Pot Hill (3), which were separated by 3.6 m, suggests that a single genotype can be at least 529–857 years old. This evidence argues that individuals within many populations of E. cordata could be very long lived, and represent genotypes remaining from previously larger populations that were not affected by small population processes.

Introgressive hybridisation is another possible mechanism for the maintenance of high levels of genetic diversity in the small populations of E. cordata. Hybridisation between co-occurring species has been suggested to increase genetic diversity through the introgression of alleles (Rieseberg et al. 2007) and has been observed in animals (Canestrelli et al. 2010) and in plants (Klier et al. 1991), including eucalypts (Butcher et al. 2002). In several eucalypt studies, the presence of F₁ hybrids has been reported to increase with decreasing population size (Field et al. 2009; Larcombe et al. 2014; Potts and Wiltshire 1997), and contemporary and historic hybridisation of E. cordata has been reported (McKinnon et al. 2004b; P. Harrison pers. obs.). However, in the well-studied case of E. cordata–E. globulus hybridisation, McKinnon et al. (2010) provided evidence of cryptic introgression of nuclear markers from E. cordata into the surrounding E. globulus, but there was no evidence of cryptic introgression from E. globulus into pure E. cordata phenotypes.

Although multiple processes may be operating to affect population genetic diversity, we suggest that, on a broad scale, the patterns of molecular genetic diversity and affinities within E. cordata can be largely explained by differences in the time since populations have been isolated. The complex spatial patterns observed reflect the co-existence of relictual and more recently derived populations that have had various degrees of isolation following multiple phases of population extinction or contraction, and expansion. The spatial patterns of genetic diversity suggest that the populations of E. cordata comprise two differentiated regions. For differentiation to occur between such sizable groups of populations, a long separation with a lack of gene flow would be required. This suggests that in the past (probably before the last glacial epoch) the distribution of the two regions in the north-east and south-west were well separated (Links 4 and 5, Fig. 6). Surrounding the genetically diverse central core of each region are peripheral, often more divergent, populations of lower genetic diversity. Populations near the centre of a species geographical range generally exhibit higher genetic diversity and lower divergence than populations on the periphery (central–peripheral hypothesis; Brussard 1984; Eckert et al. 2008), which appears to hold for both the eastern and western cores of E. cordata. Indeed, Vucetich and Waite (2003) have shown that peripheral populations can experience up to 30 times greater genetic drift than central-core populations. Although peripheral populations have traditionally been studied in relation to a continuous central core (Eckert et al. 2008; Jones et al. 2013), in the case of E. cordata, the cores are not continuous but have become recently fragmented. This hypothesis is supported by the low differentiation between populations within the ‘phantom core’, which is suggestive of a more continuous distribution in the cores in the recent past.

Populations of low microsatellite diversity tend to separate the two central core areas and form a band through the central part of the E. cordata distribution, which includes E. morrisbyi (Fig. 4a). These populations are centred on the River Derwent valley and fringes of Storm Bay (even near sea level, e.g. Populations 10 and 23), and are likely the remnants of a more widespread distribution in this area during the various Quaternary glaciations. Lowered sea levels during the Quaternary glacials would have resulted in large expanses of low-altitude habitat suitable for eucalypt forest (Fig. 6; Kirkpatrick and Fowler 1998; McKinnon et al. 2004a; Jackson 2005). The affinities of the population at Taranna (14) on the Tasman Peninsula with western populations at Moogara (14) and Chimney Pot Hill (3) (Link 3, Fig. 6) rather than to geographically closer populations to the north provides support for this hypothesised past, more-or-less continuous distribution of E. cordata in the now-flooded Storm Bay. On the other hand, the genetic affinities between Penguin Island (10), Brown Mountain (1) and the Calverts Hill E. morrisbyi population (23) (Link 1, Fig. 6) provides evidence for an even older, more ancestral distribution in this area, which probably corresponds to an older glacial cycle. The antiquity of this link is suggested by the generally high differentiation of these populations from conspecific populations and the ancestral heteroblastic features in the two E. cordata populations. The genetic distinctiveness of the two isolated E. morrisbyi populations suggests long occupancy and isolation in this region, and the genetic affinities of the Meehan Range E. cordata population (8) and Risdon Hill E. morrisbyi (22) may reflect historical nuclear gene exchange between these two proximal populations early in the speciation process.

Chloroplast genome diversity was relatively low in many E. cordata populations, with most being fixed for a specific haplotype. This pattern was previously noted by McKinnon et al. (2004b) in a smaller scale study. However, the E. cordata population at Square Mountain (13) was an exception to this and had the highest chloroplast diversity. The unique chloroplast haplotypes found in Square Mountain have not been reported in any other eucalypt (McKinnon et al. 2001b; University of Tasmania DNA database 2014, unpubl. data, comprising 1780 samples from 44 species) and were not shared with other E. cordata populations. Chloroplast sharing with co-occurring eucalypts could explain the high haplotype diversity at Square Mountain (i.e. McKinnon et al. 2004b); however, the uniqueness of these chloroplast haplotypes to this E. cordata population argues against this diversity originating from hybridisation. The high haplotype diversity at Square Mountain contrasts with surrounding populations in the north-east, which tended to have high microsatellite diversity, yet are fixed for different chloroplast haplotypes.

The variation in the chloroplast genome in the north-eastern core populations likely reflects historical population processes associated with past population isolation or bottlenecks. Newly founded or expanded populations are expected to experience severe bottlenecks, which erode the genetic diversity (Excoffier et al. 2009). Because the chloroplast genome is maternally inherited in eucalypts (Byrne et al. 1993; McKinnon et al. 2001b), it is more susceptible to genetic drift than the nuclear genome (Petit et al. 1993). This pattern of low diversity of chloroplast haplotype in the north-eastern lineage could arise if the core had been established by founder events along an eastward expansion axis or re-established from multiple small remnants from a previous era during which chloroplast diversity had been eroded. Under this scenario, the greater pollen- than seed-mediated gene flow would allow nuclear genetic diversity to be rapidly re-established as populations expand, as demonstrated in several studies of isolated tree populations (Buschbom et al. 2011; Hampe et al. 2013). The high chloroplast diversity in the more western population at Square Mountain (13) suggests that it has not experienced a major bottleneck in the past and could be the source of the eastern expansion of E. cordata subspecies cordata in this north-eastern region (Link 4, Fig. 6).

In contrast to the north-eastern lineage, genetic diversity within the south-west lineage appears to have been shaped to some extent by recent, seed-mediated migration of populations across a mid-altitude plateau, south of the Mount Wellington Range (Fig. 6). Recent colonisation and maintenance of large populations is supported by the close microsatellite affinities and high diversity, and the sharing of a common chloroplast haplotype (JS43) among most populations in this core region. The low microsatellite diversity in two geographic outliers to this core, Cape Queen Elizabeth (4) and Mount Grosse (7), suggests that these populations have lost gene flow connectivity with this core. Although having greatest microsatellite affinities to the south-western lineage, the affinities of the low-altitude populations at Cape Queen Elizabeth (4) and Electrona (15) to the populations transgressing the Derwent Valley (Cluster 3, Fig. 4; Link 3, Fig. 6) may reflect their more ancestral, intermediate status.

The marked morphological differentiation between populations of the core of the south-western molecular lineage without major, putatively neutral microsatellite differentiation implies that selection has played a role in the evolution and maintenance of the subspecies quadrangulosa. This south-western core includes populations that are fixed for quadrangular branchlets and classify as subspecies quadrangulosa (e.g. Snug Plains, 21; Sherwood Hill, 20; Combes Hill, 19), populations that have round branchlets (e.g. O’Briens Hill, 9; Mount Grosse, 7) and are of subspecies cordata, and intermediates (Electrona, 15; Herringback (top), 16). This finding raises the question of the origin of the quadrangular branchlet phenotype. A possible scenario for the evolutionary origin of the quadrangular branchlets could be ancient hybridisation coupled with selection. Quadrangular branchlets have evolved several times within Maidenaria and they are most notable in juvenile E. globulus. Recent studies have found evidence of genome-wide introgression of adaptive traits in Helianthus (Yatabe et al. 2007; Vekemans 2010; Whitney et al. 2010) and Iris (reviewed in Arnold and Martin 2009), and some genomic regions of E. globulus could have been introgressed into an E. cordata-like hybrid. However, although introgression in the reverse direction (e.g. E. cordata markers into E. globulus) has been demonstrated, no evidence of nuclear introgression into E. cordata was evident (McKinnon et al. 2010). An alternative and perhaps more likely explanation for the quadrangular branchlets of subspecies quadrangulosa is that the trait has independently evolved within E. cordata as a direct or indirect response to directional selection for increased vegetative and reproductive organ size in the cooler, wetter and more productive environments occupied by this subspecies (Potts 1989) and may be of structural significance.

In conclusion, our genetic study of E. cordata showed that most populations maintain high levels of genetic diversity relative to other localised eucalypt species. This indicates that the species appears to be ‘genetically healthy’. The genetic clustering of individuals into five distinct groups provides a genetic framework for the management of the species. The preservation of populations in each of these groups as well as all of the relict populations will be important. Although most populations are well represented in reserves, the low number of individuals in some populations suggests that their ability to respond to stochastic environmental change may be limited, in particular, the small, relict populations at Penguin Island and Brown Mountain, which have preserved signals of the early speciation stage of E. cordata. The absence of marked genetic differentiation between many morphologically differentiated populations recognised as different E. cordata subspecies argues for recent evolution of the localised subspecies quadrangulosa through directional selection. The current distribution of subspecies quadrangulosa is more restricted than subspecies cordata, and although this subspecies is currently well reserved, the status and health of unreserved populations, including the intermediates, should be monitored to ensure long-term maintenance of the genetic diversity of this endemic species.