Population genetic analysis of a medicinally significant Australian rainforest tree, Fontainea picrosperma C.T. White (Euphorbiaceae): biogeographic patterns and implications for species domestication and plantation establishment

Genetic diversity

Measures of observed heterozygosity (H _O ) were relatively low across populations, ranging from 0.298–0.487 (mean H _O   = 0.397) and were more or less concordant with levels of expected heterozygosity (H _E ) (0.264 to 0.507; mean H _E  = 0.407) calculated under conditions of Hardy-Weinberg Equilibrium (HWE) (Table 2). Consequently, combined populations of the dioecious F. picrosperma displayed an overall negligible level of inbreeding (mean F = −0.003), however individual population values ranged between F = −0.139 to 0.149, indicating a low to moderate excess of either heterozygotes or homozygotes at particular sites (Table 2).

Although the level of genetic diversity resolved in the 218 samples tested was reasonably low, the statistical confidence for individual identification using the 11 loci employed in this study was quite high (P _ID  = 1.72 × 10⁻⁶) with only two individuals from East Barron found to share the same multilocus genotype. The other 216 samples had unique multilocus genotypes. Several microsatellite markers displaying minimal polymorphism (2–3 alleles; Table 1) were removed from the analysis to assess its sensitivity to a reduction (and by inference, increase) in loci; whilst there was minimal impact on fundamental genetic diversity outputs, a considerable proportion of the discriminatory power to accurately identify individuals was lost. The validation of the ability to discriminate individuals using the complete set of 11 markers identified for this study is therefore significant with regards to future selection and breeding programs.

Population structure and gene flow

Analysis of Molecular Variance (AMOVA) found most (75 %) of the species diversity to reside within populations, with the rest of the variation due to differences between populations (PhiPT = 0.248, p = 0.001) (Additional file 1: Table S2; supporting information). Wright’s F-statistics further subdivided population differentiation into a combination of differences among individuals (F _IS  = 0.096) and populations (F _ST  = 0.153, p = 0.001), translating to a low to moderate level of historical gene flow (mean N _m  = 1.382 individuals/generation), sufficient to prevent or slow the rate of genetic drift between sites (Table 3). Pairwise population F _ST values were all significantly different from zero (p <0.001) and ranged from a level of minimal differentiation (F _ST  = 0.035; N _m  = 6.880) between the relatively proximate populations at Topaz and Towalla to a maximum distance (F _ST  = 0.302; N _m  = 0.579) between the two northern, most isolated populations, Gadgarra and East Barron (Table 3; Fig. 1). In fact, apart from a low level of contact suggesting Boonjie as the possible source population (Boonjie-Gadgarra N _m  = 1.585; Boonjie-East Barron N _m  = 1.360), neither Gadgarra nor East Barron displayed sufficient gene flow (N _m  < 1.000) with any of the other populations to prevent genetic drift [52]. Conversely, both Evelyn Highlands and Boonjie displayed evidence of genetic exchange with most other populations, supporting the hypothesis that both of these populations may be long term refugia.

	Evelyn Highlands	Boonjie	East Barron	Malanda	Topaz	Gadgarra	Towalla
Evelyn Highlands	0.000	2.122	0.850	2.677	1.985	0.956	1.291
Boonjie	0.105	0.000	1.360	2.125	2.911	1.585	1.859
East Barron	0.227	0.155	0.000	0.979	0.786	0.579	0.812
Malanda	0.085	0.105	0.203	0.000	3.672	0.682	3.104
Topaz	0.112	0.079	0.241	0.064	0.000	0.834	6.880
Gadgarra	0.207	0.136	0.302	0.268	0.231	0.000	0.745
Towalla	0.162	0.119	0.235	0.075	0.035	0.251	0.000

Mean F _ST  = 0.153. Mean N _m  = 1.382. Effective levels of past gene flow among the seven populations of F. picrosperma assessed are indicated in bold type. Values based on 999 permutations

The UPGMA cluster analysis (Fig. 2) further confirmed Gadgarra and East Barron as more divergent, with approximately 88 % and 90 % similarity, respectively, to the remaining populations of the species (Fig. 2). The most western and eastern peripheral populations of Evelyn Highlands (on and around Mt. Hypipamee, 1125 m asl) and Boonjie (on the western slopes of Mt Bartle Frere, 1622 m asl) displayed a moderate gene flow (N _m  = 2.122) strongly suggesting that similarity (F _ST  = 0.105) may be linked to their putative status as long-term interglacial refugia (Figs. 2, 3 and 4; Table 3), rather than recent gene flow per se.

Principal coordinates analysis (PCoA) detected a close genetic relationship between individuals within populations due to low levels of diversity (Fig. 3). The first three principal components were the main axes of variation as indicated by the scree plot (Additional file 2: Figure S1) and broken stick analysis according to Jackson [27], explaining a cumulative variation between individuals of only 37.17 %. The scree plot indicated a gradual decay in eigenvalues rather than a steep decline, further highlighting the low levels of diversity and genetic structure observed in this species. However, despite failing to clearly separate populations into discrete clusters, the PCoA analysis mostly concurred with the UPGMA cluster analysis and supports our hypothesis of the existence of three main groups; a western group (Evelyn Highlands), an eastern group (Boonjie) and a central group (Topaz, Towalla, Malanda and Gadgarra). For example, individuals from Boonjie and Evelyn Highlands form two separate but genetically overlapping groups that combined overlap the majority of individuals from Topaz, Towalla, Malanda and Gadgarra, which themselves cluster tightly together. Although Gadgarra clustered with the central populations, it seems to be somewhat inbred and genetically divergent from this group, containing a depauperate subset of the genetic variation found within the central populations (Table 2). In contrast, East Barron’s genetic distinctiveness was likely due to a relatively high proportion of private alleles (Table 2) and random founder effects that took place during its establishment (Figs. 2, 3 and 4). Results of the STRUCTURE analysis indicated that ln likelihoods of the data plateaued quickly from K = 3 to K = 4 (Additional file 3: Table S1, supporting information). Hence, K = 3 was selected as the best estimate of the number of genetic clusters following implementation of the Evanno et al. [17] method in STRUCTURE HARVESTER. However, additional genetic structure of biological relevance at different levels of K is also apparent (Fig. 4). While each of the three genetic clusters was present at the seven sites assessed, proportions differed substantially among the populations. Calculation of the average proportionality of each genetic cluster for each population support the UPGMA and PCoA analyses and the presence of three main groups of F. picrosperma (Fig. 4). For example, the average proportionality of each genetic cluster for the western group (Evelyn Highlands) was approximately 34 % K₁ (pink), 58 % K₂ (orange) and 8 % K₃ (blue), compared to 31 % K₁, 13 % K₂ and 56 % K₃ for the eastern group (Boonjie) and 4 % K₁, 32 % K₂ and 64 % K₃ for the central group (Topaz, Towalla, Malanda and Gadgarra). Evelyn Highlands and Boonjie therefore have a more uniform but differing spread of the three genetic clusters as compared to the central group, while one of the genetic clusters that is strongly represented in both Evelyn Highlands and Boonjie (K₁) is only minimally represented in the plateau group, together providing support to their assignment as putative refugial populations. The Mantel test found that the geographic structuring of F. picrosperma’s genetic variation did not follow a predictable pattern, and no relationship was detected between genetic and geographic distance matrices among populations (R _xy  = 0.282; r ²  = 0.0795; p > 0.05).

Results of the Bottleneck analysis did not detect any signs of recent bottlenecks in five of the seven populations assessed (p > 0.05). However, a significant (p = 0.004) heterozygosity excess at ten of the eleven loci was found in both the Malanda and Gadgarra populations. This data suggests that individuals at these sites are showing effects of disruption to ‘continuous’ populations and are no longer in mutation-drift equilibrium. These effects likely reflect their long-term isolation from populations in the two putative refugial areas (Boonjie and Evelyn Highlands) for this species and may have been further exacerbated by anthropogenic activities such as aboriginal burning since the Last Glacial Maximum and large-scale rainforest clearing in more recent times.

Genetic diversity

This is the first study to utilise microsatellites to examine genetic structure in the genus Fontainea. We investigated the levels and partitioning of genetic variation across the known range of F. picrosperma and found that the seven populations surveyed were genetically distinct despite having uniformly low levels of genetic diversity. This finding was not unexpected as many Australian plant species are characterised by low levels of genetic diversity, often as an adaptation to harsh environmental conditions [29, 51, 55], but also as a result of belonging to ancient lineages [45]. For instance, contrary to the accepted anthropomorphic view that a high level of genetic diversity bestows optimal evolutionary capability under conditions of environmental stress, James [29] found low levels of diversity in many successful species of Australia’s southwestern flora due to the purging of recombinational impedimenta (genetic load), allowing them to operate in harsh conditions at a highly adapted level. This counter-intuitive finding may also explain low genetic diversity in many of the ancient lineages in the Australian rainforest flora [22], including the results of this study for F. picrosperma. In essence, these rainforest taxa are highly adapted over long time periods to specific niches provided by the rainforest environment. As a consequence of this specialisation and niche differentiation in an essentially stable local environment, they experience only modest selection pressure during periods of climatic stability and when environmental conditions change, they retreat into the remaining environmental habitat to which they are so well adapted.

Inbreeding

In general, the results indicate extremely low levels of inbreeding (F = −0.003), which despite local populations having been isolated through glacial events, and more recently by anthropogenic habitat fragmentation, would be expected in an obligate outcrossing, dioecious species like F. picrosperma. Even though proximate trees are likely to be siblings or half-sibs, due to the limited dispersal capabilities of F. picrosperma’s relatively large drupaceous fruit, this suggests that deleterious mutations may have been purged over time, as most of the diversity resolved was between individuals within populations, not among populations.

The slight excess of heterozygosity detected in some populations suggests that recent bottlenecks with subsequent founder effects due to the expansion/contraction dynamics of small populations located outside of the main refugia may be responsible for a minor degree of genetic drift causing the random fixation of alleles. However, several populations were found to exhibit an equally slight excess of homozygosity, either as a result of the lack of overall genetic variation in the species or because of consanguineous matings. Although allelic diversity was found to be low, the fact that only two individuals shared the same multilocus genotype indicates that ‘selfing’ among proximate sibs or half-sibs was of limited occurrence; in fact these two individuals may be clones. Numerous studies have found pollen travel in continuous rainforest vegetation may be within the order of several kilometres [3]; more detailed, parent-progeny research to investigate fine-scale patterns of gene flow within wild populations, aimed at maintaining optimal among production seed crops of F. picrosperma, is required.

Population structure and gene flow

Stands of F. picrosperma occur in the upland and highland rainforests of the Atherton Tableland within a 15–20 km radius of Malanda. As such, the seven populations selected for population genetic analysis in this study likely represent a considerable proportion of the available genetic diversity within the species. It is entirely plausible that the low levels of genetic diversity and weak population structure that we have observed within F. picrosperma could merely be reflective of a random distribution of the diversity between individuals and populations. However, we believe that our observations reflect the existence of three distinct races or forms, including two long-term refugial races where suitable habitat is known to have persisted during less favourable times [25].

The population genetic structure of F. picrosperma is likely heavily influenced by the species’ life-history attributes and the effects of a long history of rainforest attrition followed by successive cycles of glacial-induced expansion and contraction upon the distribution of remaining populations. The Quaternary glacial cycles of recent geological times are known to have played a significant role in the current distributions and genetic signatures of many species [24] and based on our results this would seem to apply to F. picrosperma. Episodes of range expansion and contraction can have considerable genetic consequences [42] and the dynamics of the Wet Tropics rainforests corresponding to the glacial cycles of the Plio-Pleistocene are well documented [23, 33, 57]. Hence, the present-day configuration of F. picrosperma’s population genetic structure is likely a direct product of re-colonisation of dry sclerophyllous vegetation by tropical rainforest from refugial pockets of suitable habitat, following amelioration of the cool, dry conditions associated with past glacial cycles [8, 15, 23, 25, 33, 34, 37, 38, 51, 54, 56, 57]. It is likely that during this period several of the central populations assessed here have undergone at least some degree of geographic and genetic isolation.

The fruits of F. picrosperma disperse primarily by gravity with secondary long-distance dispersal facilitated either by hydrochory along drainage lines or zoochorous vectors [11, 14]. Populations therefore do not spread as a continuous wave of advance but rather are found as small and often isolated clumps or clusters, which may help to explain patterns in the geographical distribution of alleles. Nonetheless, the population genetic structure of F. picrosperma and the degree of historical gene flow between populations has been sufficient to maintain species’ integrity, suggesting populations were likely more continuously distributed in the past. The fact that the genus, originally described as containing a single taxon, F. pancheri, is composed of several highly similar taxa [21, 30], suggests vicariance due to habitat contraction occasioning genetic drift and the eventual loss of species cohesion may have been responsible for species divergences in the past.

The UPGMA cluster, principal coordinates and STRUCTURE analyses all provide a clear indication about the genetic distribution of this species. When combined with the genetic diversity analysis, the data show that the geographically distant (~28 km), peripheral populations of Boonjie and Evelyn Highlands, are genetically most diverse in comparison to the other populations whilst having elements of similarity, and form two genetically similar groups. Four of the remaining populations (Topaz, Towalla, Malanda and Gadgarra) form another genetic group, whereas the population at East Barron is genetically more divergent. We speculate that the populations at Evelyn Highlands and Boonjie represent two, genetically similar races or forms representing the two main refugial areas, where F. picrosperma persisted during times of sclerophyll expansion, before re-radiating out across the landscape under more favourable climatic conditions. In contrast, we suggest that the central populations of Topaz, Towalla, Malanda and Gadgarra represent a ‘plateau’ race or form that have likely expanded from small refugia during less severe climatic cycles, forming a genetically divergent race or form of F. picrosperma. East Barron appears to be derived from the elevated population at Evelyn Highlands (~1100 m asl), but is a genetically more divergent population, probably due to random founder effects. Gadgarra on the other hand, is genetically distinct, most likely as it contains no unique alleles and is somewhat inbred; essentially Gadgarra is a genetically depauperate variation of the plateau form. Despite the fact that the data suggests the presence of these three groups, it is important to highlight that the genetic diversity within F. picrosperma is low and the genetic structure between these three groups is proportionately low. In fact, the pairwise F _ST values between Evelyn Highlands and Boonjie, Evelyn Highlands and the plateau group, and Boonjie and the plateau group range from only 0.039 to 0.060. However, each value was significantly different from zero (p <0.001) and within the context of the low levels of genetic variation within this species, this is suggestive of the presence of relevant genetic structure.

It is likely that the genetic relationship between the populations can be explained not so much by linear geographic distance but by their distribution within major river catchments radiating from the putative refugial sites of Evelyn Highlands and Boonjie. However, we also recognise the possibility that our analysis could merely be reflective of a random distribution of the observed genetic diversity. Therefore, future research to test our hypothesis that two refugial races or forms and a plateau race or form of F. picrosperma exist will necessarily involve chloroplast DNA analysis and the sampling of additional individuals sourced along potential gene flow corridors, such as major river systems originating from the putative refugia at Evelyn Highlands and Boonjie.

Selection, breeding, and plantation management

Knowledge of the genetic structure of source populations, mating system and patterns of gene flow are vital to the efficient establishment and management of seed orchard plantations and the production of improved open-pollinated seed [6, 7, 39]. Although the level of microsatellite variation detected in F. picrosperma was comparatively low, high exclusion probabilities (P _ID ) confirm that these markers will be useful in future paternity analyses and breeding programs; the former to determine patterns of gene flow in natural populations that will guide plantation design of this dioecious species, and the latter to ensure maximal genetic diversity is maintained during breeding of commercially significant agronomic traits, both of which are critical aspects of developing F. picrosperma as a niche tree crop for the supply of EBC-46.

Significant variation has been observed among F. picrosperma individuals with regard to several commercially significant agronomic traits such as growth, fruit production and EBC-46 content, suggesting that the species will be ideally suited for genetic improvement to optimise production. However, even in dioecious species, the genetic diversity of seed orchards can be eroded by a number of factors including a high proportion of ‘selfs’ arising from consanguineous matings between sibs or half-sibs, and departures from random mating due to unequal contributions of individuals to seed crops [39]. In fact, obligate outcrossing species such as F. picrosperma may be particularly vulnerable to losses in reproductive fitness stemming from elevated rates of inbreeding, leading to reductions in both vigour and yield [7, 35]. Therefore, to implement suitable plantation design and management options, it is necessary to have an intimate knowledge of a species’ mating system, reproductive biology, outcrossing rate and gene flow patterns in order to maximise breeding progress whilst preserving genetic diversity [5, 6, 39].

Traditionally, the most cost-effective manner of limiting inbreeding in ex situ populations was to position individuals in such a way that the possibility of close relatives mating would be small and hope for the best, however new techniques based on the minimisation of the global probability of consanguinity by considering the genetic relationships among trees within the entire planting have been developed [20]. Microsatellites are powerful tools for tracing pollen flow using parent/progeny arrays and work is continuing in both wild and artificial populations of F. picrosperma to establish which seed source and orchard variables are most likely to govern the efficiency of production plantations.

Study site and sample collection

Fontainea picrosperma occurs on soils derived from basaltic parent materials at altitudes of 700–1200 m above sea level (asl) and is restricted to an area of approximately 30 × 30 km on the southern Atherton Tableland, Queensland, Australia. Whilst it is geographically restricted and its distribution is fragmented, the species is relatively common at a local scale where suitable habitat exists.

We sampled 218 individuals from seven F. picrosperma populations selected to cover the geographical range of the species (Fig. 1). Leaf tissue was collected from between 17 and 68 mature plants per population, dependent upon site area and the numbers of individuals present (Table 1). The location of each individual was mapped using a handheld GPS and voucher specimens from each population have been lodged at the Queensland Herbarium (BRI). Total genomic DNA was extracted from silica-dried leaf tissue using a DNeasy™ Plant Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions.

Microsatellite analysis

A detailed description of marker development using GS-FLX Titanium chemistry (Roche Applied Science; Mannheim, Germany) is given in Agostini et al. [1]. Eleven polymorphic microsatellite loci (Table 1) with consistent PCR amplification, clear allelic variation, and clarity of electrophoretic signatures were selected to assess population genetic variation. The forward primer of each locus was direct-labelled with a fluorescent dye (VIC, PET, FAM, NED). Three multiplex PCR pools (Pool 1: FP39, FP40, FP62, FP64; Pool 2: FP21, FP44, FP56; Pool 3: FP32, FP47, FP49, FP59) were amplified using Multiplex PCR Plus Kits (Qiagen). Forward and reverse primers for each multiplex pool were combined in a 10× primer mix using 1–3 μM of each primer, dependent upon PCR product fluorescence intensities. Reactions, with volumes adjusted to 10 μL, each contained 1 μL of 10× primer premix, 3.0 μL of Qiagen Multiplex Buffer (2x), 3.5 μL of ddH₂O, and 2.5 μL of template gDNA (10 ng/μL). Amplification was performed using an Eppendorf Mastercycler (Hamburg, Germany) with cycling conditions as follows: initial denaturation at 95 °C for 5 min, followed by 35 cycles of 94 °C for 30 s, 57 °C for 90 s, and 72 °C for 30 s; with a final extension at 68 °C for 10 min. PCR products were separated by capillary electrophoresis on an AB 3500 Genetic Analyser (Applied Biosystems). Fragment sizes were determined relative to an internal lane standard (GS-600 LIZ; Applied Biosystems) using GENEMARKER v. 2.4.0 (SoftGenetics LLC, PA, USA) and double-checked manually. Individuals with low or missing peaks were amplified and genotyped a second time.

Genetic diversity

Allelic frequencies for each population were generated in GenAlEx v. 6.5 [46] and used to determine population genetic parameters including: the mean number of alleles per locus (A), observed heterozygosity (H _O ), unbiased genetic diversity (H _E ), and the fixation index (F) as a measure of past inbreeding [58]. Allelic richness (A _R ) and private allelic richness (PA _R ) for each population were obtained via rarefaction using the program HP-RARE [31] to compute the mean number of alleles per locus and the frequency of private alleles within populations, based on a minimum sample size of 17 (Malanda). Polymorphic information content (PIC) and probability of identity (P _ID ), i.e., the chance of individuals sharing the same multilocus genotype, was calculated in CERVUS v. 3.0.3 [32].

Population structure and gene flow

We used a number of methods to analyse the population structure across F. picrosperma’s distribution. The average pair-wise level of genetic differentiation (F _ST ; [58]) between populations was determined using multi-locus comparisons in GenAlEx v. 6.5 [46] based on 999 permutations. As the F _ST statistic is an indirect measure of gene flow, inversely related to the effective migration rate, it was used in the following formula N _m  = 0.25 (1- F _ST )/F _ST [59] to estimate the number of migrants per generation between populations. Nei’s unbiased genetic distance (D; [41]) was calculated to examine patterns of genetic differentiation among populations. A hierarchical cluster analysis (UPGMA - unweighted pair group method with arithmetic averaging), using pairwise F _ST was performed employing 999 permutations using POPTREE2 [53]. Estimates of genetic similarity between populations were calculated from the cluster analysis.

To look for genetic relationships within and among populations, the genetic distance matrix [41] was also used in a principal coordinates analysis (PCoA; [44]). An analysis of molecular variance (AMOVA; [19]) was then applied to quantify the partitioning of genetic variation within and among populations. Both the PCoA and AMOVA were conducted using GenAlEx v. 6.5 [46]. Mantel tests [36] were used to examine the influence of between-population geographic distances on the observed patterns of genetic differentiation (Isolation by Distance; IBD) by regressing individual pairwise genetic distances against a matrix of geographic distances calculated from GPS coordinates. Levels of significance were derived from 999 random permutations using GenAlEx version 6.5 [46].

We investigated the presence of genetically differentiated groups of populations using the Bayesian genetic clustering algorithm implemented in STRUCTURE v. 2.3.4 [48]. An admixture model was applied with correlated allele frequencies and ten independent runs for each value of K (number of clusters) between 2 and 7 were performed, employing a burn-in of 100 000 followed by 500 000 Markov Chain Monte Carlo (MCMC) steps for each run. The geographic location of samples was not used in the clustering analysis. Results across each run were summarised to infer the optimal value of K using the method of Evanno et al. [17], as implemented in STRUCTURE HARVESTER web v. 0.6.93 [16], processed using CLUMPP v. 1.1.2 [28] to determine the optimal alignment of each of the ten iterations, and visualised with DISTRUCT v. 1.1 [49].

To ascertain the likelihood of recent bottlenecks due to severe range contraction and subsequent founder effects, the program BOTTLENECK v. 1.2.02 [47] was employed to investigate whether the genetic status of current small fragmented populations was outside mutation-drift equilibrium and if the excess of heterozygosity expected after a bottleneck was significant. We tested for significant departure from equilibrium of each population using Sign and Wilcoxon’s signed rank tests conducted under the assumptions of the intermediate two-phased model (TPM), due to its suitability for microsatellite data [47].