Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Using phylogenetic diversity to identify ancient rain forest refugia and diversification zones in a biodiversity hotspot

Using phylogenetic diversity to identify ancient rain forest refugia and diversification zones in... Introduction The conservation of biodiversity to date has largely focused on maximizing gain of feature diversity of some kind whether it is species, functional traits, genetic diversity or phylogenetic diversity. Although much progress has been made in the development of these types of indices, the science of maximizing biodiversity gain has operated separately from the science of understanding how species and ecosystems have evolved and dispersed across the globe through time until relatively recently. Phylogenetic measures in particular have become increasingly relevant to applied conservation (Mace et al ., ; Rodrigues et al ., ; Cadotte & Davies, ). Here, we use a dramatic example to illustrate how important it is to bring these two disciplines closer together in applied conservation prioritization contexts. We use a case study from one of the world's oldest known continuous rain forest regions, the Wet Tropics World Heritage bioregion in north‐east Queensland, Australia. We develop a new application of the phylogenetic diversity (PD) index (Faith, ), which accounts for the biogeographic origin of the component taxa, to understand patterns of biome assembly and long‐term rain forest continuity. PD is a measure of biodiversity that incorporates information on evolutionary relationships (phylogeny) between taxa, calculated as the total branch length connecting a subset of taxa on a phylogenetic tree. The phylogenetic trees are usually inferred using nucleotide sequence data, where branch lengths are proportional to the number of genetic mutations estimated to have occurred. Instead of counting the number of species in a particular location, PD sums the total branch length represented by those taxa in the phylogeny. In this sense, two areas with the same number of species may not have the same PD value. It has been argued that the value of biodiversity measures that incorporate phylogeny (such as PD) lies in the expectation that these will better predict feature diversity of organisms. Thus, maximizing PD should maximize feature diversity and hence biodiversity option value (Faith, ). Although several studies have reported a direct correlation between taxonomic diversity and PD (Polasky et al ., ; Rodrigues & Gaston, ; Sechrest et al ., ; Davies et al ., ) raising questions as to its usefulness, the residuals of a regression of taxonomic richness and PD have been shown to identify important evolutionary patterns in the landscape that are not otherwise apparent using taxonomic diversity alone (Forest et al ., ). Here, we expand upon this use of PD by incorporating data on the historical biogeography and community assembly into a model to explain the extant distribution of genera in north‐east Queensland. The Wet Tropics bioregion of north‐east Queensland (Fig. ) is globally recognized for its so‐called relict Gondwanan flora (UNESCO, ). Patches of continuous rain forest in this location have served as arks for remnant lineages from Australia's mesic past. It is generally accepted that rain forest and mesic habitat were more widespread across Australia up until c . 40–35 mya (Byrne et al ., ). Australia's separation from Antarctica initiated continent‐wide aridification with rain forest retreating to a fragmented east coast distribution. This large‐scale contraction of rain forest habitat was exacerbated by periodic expansions and contractions of rain forest coinciding with the glacial fluctuations of the Quaternary period (Hilbert et al ., ) that would have acted as an extinction filter. Notwithstanding these conditions, Australia's persisting rain forest fragments have remained stable enough to allow the survival of many ancient lineages of plants and animals. The Wet Tropics in particular contain the highest concentration of basal angiosperm groups in the world (Endress, ) with high levels of endemism including one monotypic family and 62 monotypic genera (Metcalfe & Ford, ). Further to this exceptional evolutionary legacy, accumulating molecular evidence (Crisp et al ., ; Byrne et al ., ; Sniderman & Jordan, ) now suggests that the ‘Gondwanan’ elements of the contracted rain forest flora have had to compete with an influx of species well adapted to warm, moist conditions from Southeast Asia, as a result of the Miocene ( c . 25 mya) collision of Sahul (Australian plate) and Sunda (Southeast Asian plate), facilitating intercontinental dispersal of the biota (Crayn et al., ). Phylogenetic diversity ( PD ) (left) and residuals of PD (middle) indicated by colour‐coded scale bar at 0.0625° grid cell resolution. Shaded areas (right) have maintained stable rain forest habitat since the Last Glacial Maximum. The Wet Tropics bioregion of north‐east Queensland (inset, upper right). Shaded areas indicate extant rain forest and vine thickets in Australia. This complex and ancient biotic history would have undoubtedly left a signature on the extant distribution of plants. We aimed to determine whether the integration of the PD index and data on historical biogeography could be used to explain the extant distribution of genera. Where are the hotspots of evolutionary history? Do they coincide with refugial areas thereby highlighting the museum effect (Erwin, ), is phylogenetic diversity concentrated in areas where active speciation is occurring suggesting the importance of the ‘cradle’ (May, ), or are high PD areas associated with biogeographic convergence zones indicating immigration and admixture of biotas is a key process to identify for preserving evolutionary potential? We explored these questions by first reconstructing the genus‐level phylogeny of the north‐east Queensland Wet Tropics flora using nucleotide sequences of the plastid ribulose‐1,5‐bisphosphate carboxylase/oxygenase large subunit ( rbcL ) for all genera present in a plot network spanning the region. The high amplification rate of the rbcL gene across all angiosperm families and its designation as a universal barcode for land plants (CBOL Plant Working Group ) make it an optimal choice for measuring feature diversity at a bioregional scale (Forest et al ., ). After calculating the PD at multiple spatial scales from our phylogeny, we modelled PD as a function of GR (genus richness) then mapped residuals from this model to investigate the spatial variation in PD that cannot be attributed to the effects of GR. On examination of the spatial distribution of residuals, we undertook a further round of model fitting identical to above but including two further predictor variables, altitude and an estimate of biogeographic history of the assemblage. Estimate of biogeographic history was based on per cent Sahul plate derived (Sahul) versus Sunda plate derived (Sunda) species per site. In further identifying a correlation between the spatial patterns of our results and the degree to which rain forest extent has fluctuated in association with glacial cycles in north‐east Queensland, we demonstrate a general and broad scale approach that can be used to identify ancient rain forest refugia and more recent speciation zones in hyperdiverse regions of the world. Methods Distribution and DNA sequence data We utilized presence/absence of families, genera, and species, elevation and GPS location data from a plot network of 238 0.1 hectare plots spanning the Queensland Wet Tropics bioregion. Leaf tissue was sampled from fresh leaf material or from herbarium specimens ≤ 10 years of age for at least one species per genus present in the plot network (Appendix S1 in Supporting Information). Total genomic DNA was extracted using the Machery Nagel Plant II DNA Extraction Kit at the Australian Genome Research Facility (AGRF). Amplification of the rbcLa fragment was performed following the Consortium for the Barcode of Life (CBOL) Plant Working Group PCR protocol (CBOL Plant Working Group ) using the primers rbcLa forward (ATGTCACCACAAACAGAGACTAAAGC) and rbcL a reverse (GTAAAATCAAGTCCACCRCG). Sanger sequencing was performed by the Australian Genome Research Facility and the Canadian Centre for DNA Barcoding according to their protocols and submitted to the Barcode of Life Data Systems (BOLD) database (Ratnasingham & Hebert, ). In total, 585 new rbcL sequences were generated, representing 585 genera, 43 orders and 129 families of flowering plants. Sequences for 73 genera were obtained from GenBank ( http://www.ncbi.nlm.nih.gov/genbank/ ). Phylogenetic tree and PD Consensus sequences were assembled using ChromasPro v.1.32 and DNA Baser Sequence Assembler v.3, aligned with MAFFT online v. 6 ( http://mafft.cbrc.jp/alignment/server/ ), then checked manually. The final alignment of rbcLa (560 base pairs) for 658 genera was analysed using maximum likelihood in PhyML (Guindon et al ., ) using the HKY85 substitution model with estimated gamma shape parameters and optimized topology and branch lengths. The final tree was then imported into Biodiverse (Laffan et al ., ) to calculate PD, and the richness of species and genera were plotted against PD per taxonomic unit at four different spatial scales (grid cells of 0.1 hectares, 0.0625°, 0.125° and 0.25°). To verify the accuracy of PD estimated from our original tree, we then ran a second maximum‐likelihood analysis in raxmlGUI (Silvestro & Michalak, ) with 1000 bootstrap pseudo replicates using the GTRGAMMA substitution model. As neither of the maximum‐likelihood trees were completely congruent with the Angiosperm Phylogeny Group III (Bremer et. al, ) topology, we ran a third analysis that was constrained to the APG III tree (Bremer et al. ( )). We constrained the tree by defining monophyletic groups at the level of order and above following APG III classification in beast (Drummond & Rambaut, ). We then ran a Bayesian analysis using the HKY substitution model with a lognormal relaxed clock and Yule process tree prior for 100 million iterations. The maximum clade credibility tree from this analysis and the most likely tree from the raxmlGUI analysis were then both imported into Biodiverse to compute PD as a proportion of total tree length (PD_P). These results were correlated against the results obtained for the PhyML tree. The PD_P values from all three phylogenetic methods were strongly correlated to each other (PhyML on raxmlGUI, R = 0.9976, P < 0.0001; PhyML on beast , R = 0.9681, P < 0.0001; raxmlGUI on beast , R = 0.9683, P < 0.0001) providing confidence that the results are independent of the method used to generate the phylogenetic structure underlying the calculation of PD. The APG III constrained tree was selected for display due to its substantially higher support values below the order level compared with the two unconstrained trees. The maximum clade credibility tree with posterior support values for each node is displayed as a cladogram in Fig. S1, and a phylogram of a randomly selected tree from the results of the same analysis is displayed in Fig. S2. Linear regression To examine the ability of GR to predict PD, we tested the relationship between them via regression at four spatial scales (0.1 Ha, 0.0625°, 0.125° and 0.25°). Presence/absence data from the 0.1 Ha spatial resolution were amalgamated into grid cells for the three higher spatial scales in Biodiverse. To account for any components of PD not explained by GR, we plotted the residual values from previous analyses and examined the distribution of residuals in geographic space. We separated the residuals into six classes (see ), colour coded to indicate sign (+/−) and size proportional to their relative values. The spatial analyses were performed, and the results mapped using Biodiverse and DIVA‐GIS v.5.2 (Hijmans et al ., ). Spatial autocorrelation Spatial autocorrelation in PD and residual PD was first examined via Moran's I correlogram using spatial eigenvector mapping (SEVM), generated through Principal Coordinates Neighbour Matrices (PCNM) (Diniz‐Filho & Bini, ). Geographic location information was based on decimal latitude and longitude of the centre of each grid cell, and a truncation distance (calculated in SAM – Spatial Analysis in Macroecology; Rangel et al ., ) of 0.428 decimal units was used to create spatial filters. Three eigenvector filters were chosen based on their influence on PD being both statistically significant ( P < 0.05) and having sufficient explanatory power ( r 2 > 0.2) (Huang et al ., ). Competing regression models were generated based on all possible combinations of GR and the three spatial filters and compared for their ability to describe PD based on the sample size‐corrected Akaike information criterion (AICc). Although the model with the lowest AICc score was considered to be most informative (Burnham & Anderson, ), all models having ΔAICc values ≤ 2 showed substantial support as approximating models. To visualize the geographic distribution of the component of PD that could not be described by either GR or any spatial filter, we then mapped the residuals from the best approximating model. Incorporating historical biogeography A published dataset of floristic origins data (Richardson et al ., ) was utilized which assigned three categories of floristic origin to all genera present in north‐east Queensland; (1) Sunda – plants which have dispersed southeast from the Sunda plate, (2) Sahul – plants which originated on the Sahul plate and persisted in Australia or dispersed northwest (also popularly known as ‘Gondwanan’) and (3) unresolved – lineages of plants for which data are lacking to assign a floristic origin. Here, we use the term Sahul instead of ‘Gondwanan’ to distinguish Australasian lineages from lineages that have dispersed to or dispersed back to Australia from the Sunda plate. In this sense, some Sunda lineages may have a more ancient Gondwanan origin having rafted to the Sunda region on the Indian plate. Richness of Sahul and Sunda taxa was calculated in Biodiverse. Per cent of each floristic element was then calculated by dividing these numbers by the total number of species in each plot. We then included two further predictor variables [altitude (alt)] and the proportion of species within each local community that were determined to be from Sahul lineages [p(Sahul)] in another round of model fitting with spatial filters identical to that described above. Model selection was based on AICc. We also checked for the ability of inclusion of spatial filters to remove possible influence of spatial autocorrelation by plotting residuals from the best approximating model via Moran's I correlogram (Diniz‐Filho & Bini, ). We used Spearman's rank correlation test comparing our floristic origin data per site, p(Sahul), and the rain forest stability index modelled for the Wet Tropics bioregion during the Quaternary to one hectare resolution (Hilbert et al ., ) to test for correlation between the per cent Sahul species per site and the level of rain forest continuity. Results Hotspots of PD and genus richness PD was found to be concentrated in two primary areas in north‐east Queensland (see Fig. , blue). A large southern hotspot in the southern Cairns–Cardwell lowlands (hotspot 1) and a smaller northern hotspot concentrated in the Daintree region (hotspot 2). Moderately high PD values were also observed in smaller isolated areas in the Paluma range and the northern edge of the Black Mountain Corridor. Low PD is observed between the northern and southern hotspots, south of hotspot 1 and north of hotspot 2 (Fig. , yellow). Hotspots of genus richness (GR) are nearly identical (Fig. S3) and when GR is plotted against PD, GR is an accurate predictor of PD at four spatial scales (Fig. ). The relationship between phylogenetic diversity ( PD ) and genus richness ( GR ) at four different spatial resolutions; 0.1 hectares (black), 0.065° (red), 0.125° (green) and 0.25° (blue). Significant spatial autocorrelation was associated with PD among sites (Fig. S4). In our first round of analyses, GR was identified as a significant explanatory variable in determining PD (Table ), although in all cases, the best approximating model included (at least) one spatial filter (Table ) and GR. Coefficients for GR in models were 0.05 at 0.1 Ha scale, but slightly lower for scales amalgamating sites (0.033, 0.032 and 0.028 for 0.0625°, 0.125° and 0.25° grid cells, respectively) (Table ). In all cases, a high proportion of the variation in PD could be explained by GR and spatial filters alone [ R 2 values; 87.5% (0.1 Ha); 91.9% (0.065°); 97.3% (0.125°); 97.6% (0.25°)]. A total of 12.5% of observed variation in PD from our point data (0.1 Ha) were not explained by GR and spatial filters. When this 12.5% variation (PD residuals) was mapped, a spatial pattern was revealed. Summary of the best fit multiple regression model at four spatial scales Variable b t P Scale 0.1 ha Constant 1.943 18.788 < 0.001 GR 0.05 39.673 < 0.001 SF3 1.92 4.087 < 0.001 Scale 0.0625° Constant 2.675 15.99 < 0.001 GR 0.033 34.789 < 0.001 SF2 2.466 4.735 < 0.001 Scale 0.125° Constant 2.78 15.299 < 0.001 GR 0.032 42.243 < 0.001 SF1 −0.799 −1.832 0.071 SF2 1.848 4.691 < 0.001 SF3 −2.303 −5.161 < 0.001 Scale 0.250° Constant 3.6 13.707 < 0.001 GR 0.028 34.225 < 0.001 SF1 −1.528 −3.297 0.002 SF2 2.718 6.398 < 0.001 Standardized regression coefficients ( b ), t statistics, and associated P ‐values of the best fit multiple regression model to explain phylogenetic diversity (PD) based on a single predictor (GR = generic richness), and three spatial filtering variable conducted at four different spatial resolutions. The single model presented for each spatial resolution was that with the lowest AICc value. Mapping the residuals of PD Three spatial patterns were identified from the residuals of the regression models. The area encompassing hotspot 2 had high PD overall but lower PD than expected based on genus richness (negative residuals). The most northern and southern areas and the area between the two hotspots (Black Mountain Corridor) all had low PD but higher PD than expected (positive residuals). Hotspot 1 and the Paluma range areas notably had high PD, which was higher PD than expected (positive residuals). There also appeared to be a spatial pattern correlated to elevation. Most of the positive PD residuals appeared to occur in lowland areas. Floristic origin data When the per cent Sunda and per cent Sahul species present per site were plotted against elevation, the per cent Sunda species decreased with an increase in elevation (Fig. S5), although a correlation with per cent Sahul species (Fig. S6) is not evident. Sites that are rich in Sahul elements are present in both the uplands and lowlands; however, sites that are rich in Sunda elements tend to be restricted to the lowlands. These results supported the inclusion of floristic origin and elevation data in a more complex model to explain the variation of PD. The per cent Sahul species was also significantly correlated to the rain forest stability index estimated for each site ( R s = 0.15, Z 2 = 2.25, P = 0.024, n = 216). Modelling PD with historical biogeography data In the second round of model fitting conducted for 0.1 Ha resolution, which included the two new terms altitude (Log 10 Alt) and proportion of Sahul taxa [p(Sahul)], three spatial filters were again identified as potentially important. Of the 61 candidate models, four were considered to have substantial support as approximations for PD (Table ). As in the initial analysis, all four models included both GR and SF3. All four models also included one of the new predictor variables p(Sahul). The most informative model contained these three factors only, while the other three models included these plus an additional spatial filter or altitude, either SF2 (model 2), SF1 (model 3) or log 10 Alt (model 4) (Table ). Based on Akaike weights ( W i ), we selected model 1 for more detailed analysis. The final model showed significant positive association between PD and GR, and PD and SF3. However, the standardized coefficient for p(Sahul) was negative which supports the hypothesis that increased proportion of immigrant lineages within local species assemblages is associated with higher PD values (Table ). This model was a significant improvement over the original ( R adj 2 = 0.908). Further, the spatial correlogram based on Moran's I index suggested that the inclusion of the single spatial filter was sufficient to overcome the pattern of spatial autocorrelation in the residuals of the regression (Fig. S4). Summary of model selection for the top four models describing the distribution of PD Model Adjusted R 2 AICc ΔAICc K W i GR, p(Sahul), SF3 0.908 227.114 0 4 0.304 GR, p(Sahul), SF2, SF3 0.908 228.679 1.565 5 0.139 GR, p(Sahul), SF1, SF3 0.908 228.785 1.67 5 0.132 GR, p(Sahul), logAlt, SF3 0.908 229.065 1.951 5 0.115 ΔAICc values were compared with the best fitting model. W i is the Akaike weight, K is the number of variables, including intercept. Predictor variables are GR = generic richness, p(Sahul) = proportion of species derived from Gondwanan lineages, logAlt = log10 altitude of 1Ha study site and SF1–SF3 = spatial filtering variables. Summary statistics of the final multiple regression model Variable b t P Constant 3.084 19.913 < 0.01 GR 0.05 46.135 < 0.01 p(Sahul) −1.873 −9.012 < 0.01 SF3 1.032 2.473 < 0.01 Standardized regression coefficients ( b ), t statistics, and associated P ‐values of the best fit multiple regression model identified in Table . GR = generic richness, p(Sahul) = proportion of species in community sample known to be derived from Gondwanan lineages, SF3 = spatial filtering variable. While we found that the best model did not include both Log 10 Alt and p(Sahul), one model with substantial support did include both variables (model 4) (Table ). Spatial correlation between p(Sahul) and Log 10 Alt was significant and positive ( n = 236, Pearson's r = 0.452, F corrected = 58.607, d.f. corrected = 227.773, P < 0.001), suggesting that the proportion of species in a local assemblage that are from Sahul lineages increases with increasing altitude in these forests. Equivocal lineages A proportion of taxa could not be resolved to a particular geographic origin with currently available data. The per cent unresolved taxa that were calculated as p(Sahul) and p(Sunda) (proportion of Sunda taxa) per site were plotted against p(unresolved species) per site to assess whether unresolved taxa were more likely to fall within one or the other categories (Figs S7 and S8). There was no significant spatial correlation between p(Sunda) and p(unresolved) ( n = 236, Pearson's r = 0.018, F corrected = 0.066, d.f. corrected = 209.31, P < 0.798). There was, however, a significant spatial correlation between p(Sahul) and p(unresolved), and the correlation was negative ( n = 236, Pearson's r = −0.835, F corrected = 415.47, d.f. corrected = 180.89, P < 0.001). Thus, as the proportion of all taxa identified from Sahul lineages went up, the proportion of unresolved went down, while the proportion identified as Sunda lineages had little influence on p(unresolved). This suggests that unresolved species are more likely to turn out to be of non‐Sahul origin. The results thus present a potentially conservative picture of the influence of the immigrant element on the Queensland Wet Tropics flora. It also supports our selection of p(Sahul) as opposed to p(Sunda) as a more accurate influential factor in the model test to explain the distribution of PD. Discussion The two hotpots of phylogenetic diversity (PD) and genus richness (GR) identified are consistent with two centres of endemism previously identified for north‐east Queensland (Crisp et al ., ). The correlation between hotspots of PD and hotspots of taxonomic richness (GR) is statistically significant and initially suggested that both of these areas may be significant long‐term rain forest refugia. Plotting the linear regression residuals, however, showed that there were significant differences in the evolutionary history of these two hotspots. Hotspot 1 has a higher concentration of established Indomalayan or Sunda lineages, and the Daintree region (hotspot 2) has a much stronger ‘Gondwanan’ character with a higher proportion of Sahul lineages. In general, we found that areas that have higher PD than expected based on GR are positively correlated with areas that have a high percentage of Sunda lineages and to a lesser extent also correlated with low elevation. That the per cent Sahul taxa was a more accurate predictor of residual PD than elevation explains exceptions to the elevation trend, particularly in the lowlands of hotspot 2, which had negative PD residuals. The fossil record A wealth of macro and micro fossil data supports the dominance of a widespread cool‐adapted rain forest flora for much of Australia's history (Dettmann, ; Greenwood & Christophel, ). Forests went through warmer (more tropical) and cooler (more temperate) phases but persisted until Australia's final separation from Antarctica at c. 38 mya brought about the formation of the Antarctic circumpolar current. This resulted in an increase in the latitudinal temperature gradient in the Southern Hemisphere, and in Australia, cooler temperatures and aridification caused a massive contraction of rain forest to small refugia along the eastern coast and ranges (Greenwood & Christophel, ). From the largest of these refugia in north‐east Queensland's uplands, a prevalence of Nothofagus and Podocarpaceae in the fossil record up until the early Pliocene suggests a cooler, more temperate climate than today for much of the Miocene (Kershaw, ). The extant lowland flora of the region has also been considered the closest analogue to numerous macrofossil sites of south‐eastern Australia (Christophel, , ), thus the extant flora maintains a diverse assemblage of vegetation types, many of which bear much similarity to fossil floras of different eras during Australia's history (Kershaw et al ., ). Both upland sites rich in cooler temperate affinities and lowland sites in the bioregion have been inferred as long‐term refugia (Webb & Tracey, ; Hilbert et al ., ). This is consistent with our findings of the Sahul lineage‐rich sites occurring in both the uplands and lowlands. Hotspot 2, a lowland area in the Daintree rain forest, has long been regarded as an important long‐term refugium due to its concentration of rare, narrow range endemic genera. Its distinction from hotspot 1 in the south, in terms of richness of PD and GR, however, was only evident after plotting the residuals from our analysis. Affinities of the extant flora Although the ‘Gondwanan’ heritage of this region is well accepted, the significance of Sunda or ‘intrusive’ elements has remained contentious. Previous authors have noted fossil evidence for long distance dispersals from Southeast Asia to Australia during its isolation phase in the Oligocene, but argued this invasion was insignificant (Truswell et al ., ). Our data, however, support a strong presence of the Indomalayan element in the extant flora. Although few of these lineages have left a fossil record in Australia, nearly all of the lineages designated as Sahul derived in this study are characteristic of Australia's fossil floras (Hill, ). This is consistent with the premise that lineages that have recently dispersed and adapted to the Queensland Wet Tropics are less likely to have left a fossil record. Elevation and glacial cycle effects Our results indicate that Sunda plant lineages have more successfully colonized the lowland areas of Australia's Wet Tropics, although with notable exceptions in the northern Daintree area, whereas the upland areas have retained a higher percentage of Sahul lineages. Although the proportion of species of Sunda origin decreases with elevation, the Sahul element is abundant at higher elevations and also in specific localities in the lowlands. This is consistent with Australia's fossil record, which attests to the resemblance of both extant upland and lowland sites to fossil floras (see above). Previous authors have argued that invading Sunda lineages were unlikely to have been able to compete with an already established rain forest flora in Australia, one that was adapted to habitat in gradual decline as aridification swept across the continent (Adam, ). However, periodic glacial cycles forced rain forest to contract and re‐expand. The re‐expansion periods would have provided opportunities for invading lineages (Richardson et al ., ). We found direct support for this hypothesis in the correlation between the per cent Sahul species per site and the level of rain forest stability throughout the Quaternary period. Areas that were less stable through time had a higher proportion of immigrant lineages. Figure includes a summary map showing areas of continuous rain forest during the Quaternary based on a previously published study that modelled rain forest expansion and contraction in north‐east Queensland (Hilbert et al ., ). Our PD hotspot 1 clearly overlaps with large areas of low stability (white). These areas are currently rain forest but were not consistently covered by rain forest vegetation over the Last Interglacial cycle. Hotspot 1 contains high stability areas (black) but they are notably fragmented and smaller in comparison with the larger refugia areas directly to the north and in hotspot 2. Hotspot 2 is notably intact with a high proportion of high stability areas and lower concentration of introduced lineages. The Black Mountain Corridor, the Cooktown area north of hotspot 2, and the region from the Paluma range south to Townsville are all mostly low stability and high PD residual areas with a high concentration of introduced lineages. These patterns lend support to the view that glacial cycles influenced the establishment of the Sunda element in north‐east Queensland and can help explain the distribution of the extant species assemblages. In addition and perhaps most striking are the patterns visible in the high stability area (black) directly north of hotspot 1 (See Fig. , Atherton‐Bellenden Ker uplands). This area, mostly upland habitat stretching from the Atherton Tablelands to the Bellenden Ker and Lamb mountain ranges, has notably low PD but a strong Sahul or ancient ‘Gondwanan’ character. This case clearly illustrates the limitations of using diversity measures (species or phylogenetic) alone. Both of the two largest rain forest refugia or high stability zones in north‐east Queensland are mostly comprised of Sahul floristic elements, but one is a diversity hotspot rich in species, whereas the other is notably less diverse perhaps due to its more upland character. Of the two major diversity hotspots, hotspot 2 is clearly a long standing ancient rain forest refugium whilst the other, hotspot 1, is diverse due to the intermixing of refugial Sahul elements and an exogenous flora that likely became established there during periods when rain forest re‐expanded from refugia. Conservation implications Although PD has mostly been promoted for its use in guiding conservation policy by identifying priority areas, recent studies (Davies & Buckley, ; Fritz & Rahbek, ; Kissling et al ., ) have identified and interpreted links between PD and historical processes shaping extant species assemblages. Our study integrates historical biogeography data into a large‐scale PD analysis for an entire flora and demonstrates how important this is for interpreting biodiversity patterns in a species‐rich region. In this case, we identified distinct zones of both evolutionary relicts and fronts. But which are more important? ancient refugia of ‘Gondwanan’ heritage or zones of recent intercontinental intermixing and evolutionary potential? Although there is no clear answer to this question, it brings up important points for the field of conservation. Simply choosing a site or a series of sites that maximize PD in this scenario misses a fundamental story and the underlying processes that created the PD to begin with. In this case, maximizing PD without historical interpretation would deprioritize Australia's largest remaining tropical rain forest refugium (the Atherton‐Bellenden Ker uplands Fig. ). This area has low genus diversity but has remained stable enough through time to serve as an ark for the survival of relict lineages from Australia's ancient past. Being one of the last major fragments of continuous rain forest in Australia with a history stretching back over 40 million years, an approach which deprioritizes this area in favour of areas with high diversity due to the incursion of foreign lineages from Southeast Asia in comparatively more recent times may result in undesirable or unintended conservation outcomes. Alternatively, if the origin of lineages that dispersed into Australia comparatively recently was equally ancient or if the habitat of origin became extinct, these incursion zones (high PD residuals, Fig. ) could be viewed as a different type of refugium. Indeed, with the unabated habitat loss occurring in Southeast Asia and other tropical regions, this logic may prove to be a useful argument in other case studies around the world. Fortunately, the majority of extant rain forest in north‐east Queensland is well protected, thus we have explored these complex questions without the ‘agony of choice’. Many other species‐rich tropical regions of the world are not in the same position and may benefit from the methods proposed here to help make informed long‐term management decisions under a range of different conservation criteria. In Fig. , we outline a generalized summary of our approach so that it may be tailored and applied in other regions. Where generating new community phylogenies from sequence data is unfeasible, bioinformatics tools such as Phylomatic (Webb & Donoghue, ) can be used to generate trees suitable for estimating PD. As data on historical biogeography accumulates, the biogeographical origin of groups can be determined at a minimal cost with a species list and a dedicated literature review. This integrated approach may be used to verify or provide further support for other rain forest refugia inferred using traditional approaches, such as climate modelling and population genetics for individual species, and can rapidly advance knowledge on the evolutionary history of large species assemblages in other parts of the globe. Summary of our methodology combining phylogenetic diversity ( PD ) and historical biogeography data. Acknowledgements We thank the University of Adelaide, Australian Tropical Herbarium and Barcode of Life Data Systems molecular biology labs for helping generate the sequence data for this project. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Diversity and Distributions Wiley

Using phylogenetic diversity to identify ancient rain forest refugia and diversification zones in a biodiversity hotspot

Download PDF
11 pages

Loading next page...
 
/lp/wiley/using-phylogenetic-diversity-to-identify-ancient-rain-forest-refugia-sXergnpPxB

References (52)

Publisher
Wiley
Copyright
Copyright © 2015 John Wiley & Sons Ltd
ISSN
1366-9516
eISSN
1472-4642
DOI
10.1111/ddi.12266
Publisher site
See Article on Publisher Site

Abstract

Introduction The conservation of biodiversity to date has largely focused on maximizing gain of feature diversity of some kind whether it is species, functional traits, genetic diversity or phylogenetic diversity. Although much progress has been made in the development of these types of indices, the science of maximizing biodiversity gain has operated separately from the science of understanding how species and ecosystems have evolved and dispersed across the globe through time until relatively recently. Phylogenetic measures in particular have become increasingly relevant to applied conservation (Mace et al ., ; Rodrigues et al ., ; Cadotte & Davies, ). Here, we use a dramatic example to illustrate how important it is to bring these two disciplines closer together in applied conservation prioritization contexts. We use a case study from one of the world's oldest known continuous rain forest regions, the Wet Tropics World Heritage bioregion in north‐east Queensland, Australia. We develop a new application of the phylogenetic diversity (PD) index (Faith, ), which accounts for the biogeographic origin of the component taxa, to understand patterns of biome assembly and long‐term rain forest continuity. PD is a measure of biodiversity that incorporates information on evolutionary relationships (phylogeny) between taxa, calculated as the total branch length connecting a subset of taxa on a phylogenetic tree. The phylogenetic trees are usually inferred using nucleotide sequence data, where branch lengths are proportional to the number of genetic mutations estimated to have occurred. Instead of counting the number of species in a particular location, PD sums the total branch length represented by those taxa in the phylogeny. In this sense, two areas with the same number of species may not have the same PD value. It has been argued that the value of biodiversity measures that incorporate phylogeny (such as PD) lies in the expectation that these will better predict feature diversity of organisms. Thus, maximizing PD should maximize feature diversity and hence biodiversity option value (Faith, ). Although several studies have reported a direct correlation between taxonomic diversity and PD (Polasky et al ., ; Rodrigues & Gaston, ; Sechrest et al ., ; Davies et al ., ) raising questions as to its usefulness, the residuals of a regression of taxonomic richness and PD have been shown to identify important evolutionary patterns in the landscape that are not otherwise apparent using taxonomic diversity alone (Forest et al ., ). Here, we expand upon this use of PD by incorporating data on the historical biogeography and community assembly into a model to explain the extant distribution of genera in north‐east Queensland. The Wet Tropics bioregion of north‐east Queensland (Fig. ) is globally recognized for its so‐called relict Gondwanan flora (UNESCO, ). Patches of continuous rain forest in this location have served as arks for remnant lineages from Australia's mesic past. It is generally accepted that rain forest and mesic habitat were more widespread across Australia up until c . 40–35 mya (Byrne et al ., ). Australia's separation from Antarctica initiated continent‐wide aridification with rain forest retreating to a fragmented east coast distribution. This large‐scale contraction of rain forest habitat was exacerbated by periodic expansions and contractions of rain forest coinciding with the glacial fluctuations of the Quaternary period (Hilbert et al ., ) that would have acted as an extinction filter. Notwithstanding these conditions, Australia's persisting rain forest fragments have remained stable enough to allow the survival of many ancient lineages of plants and animals. The Wet Tropics in particular contain the highest concentration of basal angiosperm groups in the world (Endress, ) with high levels of endemism including one monotypic family and 62 monotypic genera (Metcalfe & Ford, ). Further to this exceptional evolutionary legacy, accumulating molecular evidence (Crisp et al ., ; Byrne et al ., ; Sniderman & Jordan, ) now suggests that the ‘Gondwanan’ elements of the contracted rain forest flora have had to compete with an influx of species well adapted to warm, moist conditions from Southeast Asia, as a result of the Miocene ( c . 25 mya) collision of Sahul (Australian plate) and Sunda (Southeast Asian plate), facilitating intercontinental dispersal of the biota (Crayn et al., ). Phylogenetic diversity ( PD ) (left) and residuals of PD (middle) indicated by colour‐coded scale bar at 0.0625° grid cell resolution. Shaded areas (right) have maintained stable rain forest habitat since the Last Glacial Maximum. The Wet Tropics bioregion of north‐east Queensland (inset, upper right). Shaded areas indicate extant rain forest and vine thickets in Australia. This complex and ancient biotic history would have undoubtedly left a signature on the extant distribution of plants. We aimed to determine whether the integration of the PD index and data on historical biogeography could be used to explain the extant distribution of genera. Where are the hotspots of evolutionary history? Do they coincide with refugial areas thereby highlighting the museum effect (Erwin, ), is phylogenetic diversity concentrated in areas where active speciation is occurring suggesting the importance of the ‘cradle’ (May, ), or are high PD areas associated with biogeographic convergence zones indicating immigration and admixture of biotas is a key process to identify for preserving evolutionary potential? We explored these questions by first reconstructing the genus‐level phylogeny of the north‐east Queensland Wet Tropics flora using nucleotide sequences of the plastid ribulose‐1,5‐bisphosphate carboxylase/oxygenase large subunit ( rbcL ) for all genera present in a plot network spanning the region. The high amplification rate of the rbcL gene across all angiosperm families and its designation as a universal barcode for land plants (CBOL Plant Working Group ) make it an optimal choice for measuring feature diversity at a bioregional scale (Forest et al ., ). After calculating the PD at multiple spatial scales from our phylogeny, we modelled PD as a function of GR (genus richness) then mapped residuals from this model to investigate the spatial variation in PD that cannot be attributed to the effects of GR. On examination of the spatial distribution of residuals, we undertook a further round of model fitting identical to above but including two further predictor variables, altitude and an estimate of biogeographic history of the assemblage. Estimate of biogeographic history was based on per cent Sahul plate derived (Sahul) versus Sunda plate derived (Sunda) species per site. In further identifying a correlation between the spatial patterns of our results and the degree to which rain forest extent has fluctuated in association with glacial cycles in north‐east Queensland, we demonstrate a general and broad scale approach that can be used to identify ancient rain forest refugia and more recent speciation zones in hyperdiverse regions of the world. Methods Distribution and DNA sequence data We utilized presence/absence of families, genera, and species, elevation and GPS location data from a plot network of 238 0.1 hectare plots spanning the Queensland Wet Tropics bioregion. Leaf tissue was sampled from fresh leaf material or from herbarium specimens ≤ 10 years of age for at least one species per genus present in the plot network (Appendix S1 in Supporting Information). Total genomic DNA was extracted using the Machery Nagel Plant II DNA Extraction Kit at the Australian Genome Research Facility (AGRF). Amplification of the rbcLa fragment was performed following the Consortium for the Barcode of Life (CBOL) Plant Working Group PCR protocol (CBOL Plant Working Group ) using the primers rbcLa forward (ATGTCACCACAAACAGAGACTAAAGC) and rbcL a reverse (GTAAAATCAAGTCCACCRCG). Sanger sequencing was performed by the Australian Genome Research Facility and the Canadian Centre for DNA Barcoding according to their protocols and submitted to the Barcode of Life Data Systems (BOLD) database (Ratnasingham & Hebert, ). In total, 585 new rbcL sequences were generated, representing 585 genera, 43 orders and 129 families of flowering plants. Sequences for 73 genera were obtained from GenBank ( http://www.ncbi.nlm.nih.gov/genbank/ ). Phylogenetic tree and PD Consensus sequences were assembled using ChromasPro v.1.32 and DNA Baser Sequence Assembler v.3, aligned with MAFFT online v. 6 ( http://mafft.cbrc.jp/alignment/server/ ), then checked manually. The final alignment of rbcLa (560 base pairs) for 658 genera was analysed using maximum likelihood in PhyML (Guindon et al ., ) using the HKY85 substitution model with estimated gamma shape parameters and optimized topology and branch lengths. The final tree was then imported into Biodiverse (Laffan et al ., ) to calculate PD, and the richness of species and genera were plotted against PD per taxonomic unit at four different spatial scales (grid cells of 0.1 hectares, 0.0625°, 0.125° and 0.25°). To verify the accuracy of PD estimated from our original tree, we then ran a second maximum‐likelihood analysis in raxmlGUI (Silvestro & Michalak, ) with 1000 bootstrap pseudo replicates using the GTRGAMMA substitution model. As neither of the maximum‐likelihood trees were completely congruent with the Angiosperm Phylogeny Group III (Bremer et. al, ) topology, we ran a third analysis that was constrained to the APG III tree (Bremer et al. ( )). We constrained the tree by defining monophyletic groups at the level of order and above following APG III classification in beast (Drummond & Rambaut, ). We then ran a Bayesian analysis using the HKY substitution model with a lognormal relaxed clock and Yule process tree prior for 100 million iterations. The maximum clade credibility tree from this analysis and the most likely tree from the raxmlGUI analysis were then both imported into Biodiverse to compute PD as a proportion of total tree length (PD_P). These results were correlated against the results obtained for the PhyML tree. The PD_P values from all three phylogenetic methods were strongly correlated to each other (PhyML on raxmlGUI, R = 0.9976, P < 0.0001; PhyML on beast , R = 0.9681, P < 0.0001; raxmlGUI on beast , R = 0.9683, P < 0.0001) providing confidence that the results are independent of the method used to generate the phylogenetic structure underlying the calculation of PD. The APG III constrained tree was selected for display due to its substantially higher support values below the order level compared with the two unconstrained trees. The maximum clade credibility tree with posterior support values for each node is displayed as a cladogram in Fig. S1, and a phylogram of a randomly selected tree from the results of the same analysis is displayed in Fig. S2. Linear regression To examine the ability of GR to predict PD, we tested the relationship between them via regression at four spatial scales (0.1 Ha, 0.0625°, 0.125° and 0.25°). Presence/absence data from the 0.1 Ha spatial resolution were amalgamated into grid cells for the three higher spatial scales in Biodiverse. To account for any components of PD not explained by GR, we plotted the residual values from previous analyses and examined the distribution of residuals in geographic space. We separated the residuals into six classes (see ), colour coded to indicate sign (+/−) and size proportional to their relative values. The spatial analyses were performed, and the results mapped using Biodiverse and DIVA‐GIS v.5.2 (Hijmans et al ., ). Spatial autocorrelation Spatial autocorrelation in PD and residual PD was first examined via Moran's I correlogram using spatial eigenvector mapping (SEVM), generated through Principal Coordinates Neighbour Matrices (PCNM) (Diniz‐Filho & Bini, ). Geographic location information was based on decimal latitude and longitude of the centre of each grid cell, and a truncation distance (calculated in SAM – Spatial Analysis in Macroecology; Rangel et al ., ) of 0.428 decimal units was used to create spatial filters. Three eigenvector filters were chosen based on their influence on PD being both statistically significant ( P < 0.05) and having sufficient explanatory power ( r 2 > 0.2) (Huang et al ., ). Competing regression models were generated based on all possible combinations of GR and the three spatial filters and compared for their ability to describe PD based on the sample size‐corrected Akaike information criterion (AICc). Although the model with the lowest AICc score was considered to be most informative (Burnham & Anderson, ), all models having ΔAICc values ≤ 2 showed substantial support as approximating models. To visualize the geographic distribution of the component of PD that could not be described by either GR or any spatial filter, we then mapped the residuals from the best approximating model. Incorporating historical biogeography A published dataset of floristic origins data (Richardson et al ., ) was utilized which assigned three categories of floristic origin to all genera present in north‐east Queensland; (1) Sunda – plants which have dispersed southeast from the Sunda plate, (2) Sahul – plants which originated on the Sahul plate and persisted in Australia or dispersed northwest (also popularly known as ‘Gondwanan’) and (3) unresolved – lineages of plants for which data are lacking to assign a floristic origin. Here, we use the term Sahul instead of ‘Gondwanan’ to distinguish Australasian lineages from lineages that have dispersed to or dispersed back to Australia from the Sunda plate. In this sense, some Sunda lineages may have a more ancient Gondwanan origin having rafted to the Sunda region on the Indian plate. Richness of Sahul and Sunda taxa was calculated in Biodiverse. Per cent of each floristic element was then calculated by dividing these numbers by the total number of species in each plot. We then included two further predictor variables [altitude (alt)] and the proportion of species within each local community that were determined to be from Sahul lineages [p(Sahul)] in another round of model fitting with spatial filters identical to that described above. Model selection was based on AICc. We also checked for the ability of inclusion of spatial filters to remove possible influence of spatial autocorrelation by plotting residuals from the best approximating model via Moran's I correlogram (Diniz‐Filho & Bini, ). We used Spearman's rank correlation test comparing our floristic origin data per site, p(Sahul), and the rain forest stability index modelled for the Wet Tropics bioregion during the Quaternary to one hectare resolution (Hilbert et al ., ) to test for correlation between the per cent Sahul species per site and the level of rain forest continuity. Results Hotspots of PD and genus richness PD was found to be concentrated in two primary areas in north‐east Queensland (see Fig. , blue). A large southern hotspot in the southern Cairns–Cardwell lowlands (hotspot 1) and a smaller northern hotspot concentrated in the Daintree region (hotspot 2). Moderately high PD values were also observed in smaller isolated areas in the Paluma range and the northern edge of the Black Mountain Corridor. Low PD is observed between the northern and southern hotspots, south of hotspot 1 and north of hotspot 2 (Fig. , yellow). Hotspots of genus richness (GR) are nearly identical (Fig. S3) and when GR is plotted against PD, GR is an accurate predictor of PD at four spatial scales (Fig. ). The relationship between phylogenetic diversity ( PD ) and genus richness ( GR ) at four different spatial resolutions; 0.1 hectares (black), 0.065° (red), 0.125° (green) and 0.25° (blue). Significant spatial autocorrelation was associated with PD among sites (Fig. S4). In our first round of analyses, GR was identified as a significant explanatory variable in determining PD (Table ), although in all cases, the best approximating model included (at least) one spatial filter (Table ) and GR. Coefficients for GR in models were 0.05 at 0.1 Ha scale, but slightly lower for scales amalgamating sites (0.033, 0.032 and 0.028 for 0.0625°, 0.125° and 0.25° grid cells, respectively) (Table ). In all cases, a high proportion of the variation in PD could be explained by GR and spatial filters alone [ R 2 values; 87.5% (0.1 Ha); 91.9% (0.065°); 97.3% (0.125°); 97.6% (0.25°)]. A total of 12.5% of observed variation in PD from our point data (0.1 Ha) were not explained by GR and spatial filters. When this 12.5% variation (PD residuals) was mapped, a spatial pattern was revealed. Summary of the best fit multiple regression model at four spatial scales Variable b t P Scale 0.1 ha Constant 1.943 18.788 < 0.001 GR 0.05 39.673 < 0.001 SF3 1.92 4.087 < 0.001 Scale 0.0625° Constant 2.675 15.99 < 0.001 GR 0.033 34.789 < 0.001 SF2 2.466 4.735 < 0.001 Scale 0.125° Constant 2.78 15.299 < 0.001 GR 0.032 42.243 < 0.001 SF1 −0.799 −1.832 0.071 SF2 1.848 4.691 < 0.001 SF3 −2.303 −5.161 < 0.001 Scale 0.250° Constant 3.6 13.707 < 0.001 GR 0.028 34.225 < 0.001 SF1 −1.528 −3.297 0.002 SF2 2.718 6.398 < 0.001 Standardized regression coefficients ( b ), t statistics, and associated P ‐values of the best fit multiple regression model to explain phylogenetic diversity (PD) based on a single predictor (GR = generic richness), and three spatial filtering variable conducted at four different spatial resolutions. The single model presented for each spatial resolution was that with the lowest AICc value. Mapping the residuals of PD Three spatial patterns were identified from the residuals of the regression models. The area encompassing hotspot 2 had high PD overall but lower PD than expected based on genus richness (negative residuals). The most northern and southern areas and the area between the two hotspots (Black Mountain Corridor) all had low PD but higher PD than expected (positive residuals). Hotspot 1 and the Paluma range areas notably had high PD, which was higher PD than expected (positive residuals). There also appeared to be a spatial pattern correlated to elevation. Most of the positive PD residuals appeared to occur in lowland areas. Floristic origin data When the per cent Sunda and per cent Sahul species present per site were plotted against elevation, the per cent Sunda species decreased with an increase in elevation (Fig. S5), although a correlation with per cent Sahul species (Fig. S6) is not evident. Sites that are rich in Sahul elements are present in both the uplands and lowlands; however, sites that are rich in Sunda elements tend to be restricted to the lowlands. These results supported the inclusion of floristic origin and elevation data in a more complex model to explain the variation of PD. The per cent Sahul species was also significantly correlated to the rain forest stability index estimated for each site ( R s = 0.15, Z 2 = 2.25, P = 0.024, n = 216). Modelling PD with historical biogeography data In the second round of model fitting conducted for 0.1 Ha resolution, which included the two new terms altitude (Log 10 Alt) and proportion of Sahul taxa [p(Sahul)], three spatial filters were again identified as potentially important. Of the 61 candidate models, four were considered to have substantial support as approximations for PD (Table ). As in the initial analysis, all four models included both GR and SF3. All four models also included one of the new predictor variables p(Sahul). The most informative model contained these three factors only, while the other three models included these plus an additional spatial filter or altitude, either SF2 (model 2), SF1 (model 3) or log 10 Alt (model 4) (Table ). Based on Akaike weights ( W i ), we selected model 1 for more detailed analysis. The final model showed significant positive association between PD and GR, and PD and SF3. However, the standardized coefficient for p(Sahul) was negative which supports the hypothesis that increased proportion of immigrant lineages within local species assemblages is associated with higher PD values (Table ). This model was a significant improvement over the original ( R adj 2 = 0.908). Further, the spatial correlogram based on Moran's I index suggested that the inclusion of the single spatial filter was sufficient to overcome the pattern of spatial autocorrelation in the residuals of the regression (Fig. S4). Summary of model selection for the top four models describing the distribution of PD Model Adjusted R 2 AICc ΔAICc K W i GR, p(Sahul), SF3 0.908 227.114 0 4 0.304 GR, p(Sahul), SF2, SF3 0.908 228.679 1.565 5 0.139 GR, p(Sahul), SF1, SF3 0.908 228.785 1.67 5 0.132 GR, p(Sahul), logAlt, SF3 0.908 229.065 1.951 5 0.115 ΔAICc values were compared with the best fitting model. W i is the Akaike weight, K is the number of variables, including intercept. Predictor variables are GR = generic richness, p(Sahul) = proportion of species derived from Gondwanan lineages, logAlt = log10 altitude of 1Ha study site and SF1–SF3 = spatial filtering variables. Summary statistics of the final multiple regression model Variable b t P Constant 3.084 19.913 < 0.01 GR 0.05 46.135 < 0.01 p(Sahul) −1.873 −9.012 < 0.01 SF3 1.032 2.473 < 0.01 Standardized regression coefficients ( b ), t statistics, and associated P ‐values of the best fit multiple regression model identified in Table . GR = generic richness, p(Sahul) = proportion of species in community sample known to be derived from Gondwanan lineages, SF3 = spatial filtering variable. While we found that the best model did not include both Log 10 Alt and p(Sahul), one model with substantial support did include both variables (model 4) (Table ). Spatial correlation between p(Sahul) and Log 10 Alt was significant and positive ( n = 236, Pearson's r = 0.452, F corrected = 58.607, d.f. corrected = 227.773, P < 0.001), suggesting that the proportion of species in a local assemblage that are from Sahul lineages increases with increasing altitude in these forests. Equivocal lineages A proportion of taxa could not be resolved to a particular geographic origin with currently available data. The per cent unresolved taxa that were calculated as p(Sahul) and p(Sunda) (proportion of Sunda taxa) per site were plotted against p(unresolved species) per site to assess whether unresolved taxa were more likely to fall within one or the other categories (Figs S7 and S8). There was no significant spatial correlation between p(Sunda) and p(unresolved) ( n = 236, Pearson's r = 0.018, F corrected = 0.066, d.f. corrected = 209.31, P < 0.798). There was, however, a significant spatial correlation between p(Sahul) and p(unresolved), and the correlation was negative ( n = 236, Pearson's r = −0.835, F corrected = 415.47, d.f. corrected = 180.89, P < 0.001). Thus, as the proportion of all taxa identified from Sahul lineages went up, the proportion of unresolved went down, while the proportion identified as Sunda lineages had little influence on p(unresolved). This suggests that unresolved species are more likely to turn out to be of non‐Sahul origin. The results thus present a potentially conservative picture of the influence of the immigrant element on the Queensland Wet Tropics flora. It also supports our selection of p(Sahul) as opposed to p(Sunda) as a more accurate influential factor in the model test to explain the distribution of PD. Discussion The two hotpots of phylogenetic diversity (PD) and genus richness (GR) identified are consistent with two centres of endemism previously identified for north‐east Queensland (Crisp et al ., ). The correlation between hotspots of PD and hotspots of taxonomic richness (GR) is statistically significant and initially suggested that both of these areas may be significant long‐term rain forest refugia. Plotting the linear regression residuals, however, showed that there were significant differences in the evolutionary history of these two hotspots. Hotspot 1 has a higher concentration of established Indomalayan or Sunda lineages, and the Daintree region (hotspot 2) has a much stronger ‘Gondwanan’ character with a higher proportion of Sahul lineages. In general, we found that areas that have higher PD than expected based on GR are positively correlated with areas that have a high percentage of Sunda lineages and to a lesser extent also correlated with low elevation. That the per cent Sahul taxa was a more accurate predictor of residual PD than elevation explains exceptions to the elevation trend, particularly in the lowlands of hotspot 2, which had negative PD residuals. The fossil record A wealth of macro and micro fossil data supports the dominance of a widespread cool‐adapted rain forest flora for much of Australia's history (Dettmann, ; Greenwood & Christophel, ). Forests went through warmer (more tropical) and cooler (more temperate) phases but persisted until Australia's final separation from Antarctica at c. 38 mya brought about the formation of the Antarctic circumpolar current. This resulted in an increase in the latitudinal temperature gradient in the Southern Hemisphere, and in Australia, cooler temperatures and aridification caused a massive contraction of rain forest to small refugia along the eastern coast and ranges (Greenwood & Christophel, ). From the largest of these refugia in north‐east Queensland's uplands, a prevalence of Nothofagus and Podocarpaceae in the fossil record up until the early Pliocene suggests a cooler, more temperate climate than today for much of the Miocene (Kershaw, ). The extant lowland flora of the region has also been considered the closest analogue to numerous macrofossil sites of south‐eastern Australia (Christophel, , ), thus the extant flora maintains a diverse assemblage of vegetation types, many of which bear much similarity to fossil floras of different eras during Australia's history (Kershaw et al ., ). Both upland sites rich in cooler temperate affinities and lowland sites in the bioregion have been inferred as long‐term refugia (Webb & Tracey, ; Hilbert et al ., ). This is consistent with our findings of the Sahul lineage‐rich sites occurring in both the uplands and lowlands. Hotspot 2, a lowland area in the Daintree rain forest, has long been regarded as an important long‐term refugium due to its concentration of rare, narrow range endemic genera. Its distinction from hotspot 1 in the south, in terms of richness of PD and GR, however, was only evident after plotting the residuals from our analysis. Affinities of the extant flora Although the ‘Gondwanan’ heritage of this region is well accepted, the significance of Sunda or ‘intrusive’ elements has remained contentious. Previous authors have noted fossil evidence for long distance dispersals from Southeast Asia to Australia during its isolation phase in the Oligocene, but argued this invasion was insignificant (Truswell et al ., ). Our data, however, support a strong presence of the Indomalayan element in the extant flora. Although few of these lineages have left a fossil record in Australia, nearly all of the lineages designated as Sahul derived in this study are characteristic of Australia's fossil floras (Hill, ). This is consistent with the premise that lineages that have recently dispersed and adapted to the Queensland Wet Tropics are less likely to have left a fossil record. Elevation and glacial cycle effects Our results indicate that Sunda plant lineages have more successfully colonized the lowland areas of Australia's Wet Tropics, although with notable exceptions in the northern Daintree area, whereas the upland areas have retained a higher percentage of Sahul lineages. Although the proportion of species of Sunda origin decreases with elevation, the Sahul element is abundant at higher elevations and also in specific localities in the lowlands. This is consistent with Australia's fossil record, which attests to the resemblance of both extant upland and lowland sites to fossil floras (see above). Previous authors have argued that invading Sunda lineages were unlikely to have been able to compete with an already established rain forest flora in Australia, one that was adapted to habitat in gradual decline as aridification swept across the continent (Adam, ). However, periodic glacial cycles forced rain forest to contract and re‐expand. The re‐expansion periods would have provided opportunities for invading lineages (Richardson et al ., ). We found direct support for this hypothesis in the correlation between the per cent Sahul species per site and the level of rain forest stability throughout the Quaternary period. Areas that were less stable through time had a higher proportion of immigrant lineages. Figure includes a summary map showing areas of continuous rain forest during the Quaternary based on a previously published study that modelled rain forest expansion and contraction in north‐east Queensland (Hilbert et al ., ). Our PD hotspot 1 clearly overlaps with large areas of low stability (white). These areas are currently rain forest but were not consistently covered by rain forest vegetation over the Last Interglacial cycle. Hotspot 1 contains high stability areas (black) but they are notably fragmented and smaller in comparison with the larger refugia areas directly to the north and in hotspot 2. Hotspot 2 is notably intact with a high proportion of high stability areas and lower concentration of introduced lineages. The Black Mountain Corridor, the Cooktown area north of hotspot 2, and the region from the Paluma range south to Townsville are all mostly low stability and high PD residual areas with a high concentration of introduced lineages. These patterns lend support to the view that glacial cycles influenced the establishment of the Sunda element in north‐east Queensland and can help explain the distribution of the extant species assemblages. In addition and perhaps most striking are the patterns visible in the high stability area (black) directly north of hotspot 1 (See Fig. , Atherton‐Bellenden Ker uplands). This area, mostly upland habitat stretching from the Atherton Tablelands to the Bellenden Ker and Lamb mountain ranges, has notably low PD but a strong Sahul or ancient ‘Gondwanan’ character. This case clearly illustrates the limitations of using diversity measures (species or phylogenetic) alone. Both of the two largest rain forest refugia or high stability zones in north‐east Queensland are mostly comprised of Sahul floristic elements, but one is a diversity hotspot rich in species, whereas the other is notably less diverse perhaps due to its more upland character. Of the two major diversity hotspots, hotspot 2 is clearly a long standing ancient rain forest refugium whilst the other, hotspot 1, is diverse due to the intermixing of refugial Sahul elements and an exogenous flora that likely became established there during periods when rain forest re‐expanded from refugia. Conservation implications Although PD has mostly been promoted for its use in guiding conservation policy by identifying priority areas, recent studies (Davies & Buckley, ; Fritz & Rahbek, ; Kissling et al ., ) have identified and interpreted links between PD and historical processes shaping extant species assemblages. Our study integrates historical biogeography data into a large‐scale PD analysis for an entire flora and demonstrates how important this is for interpreting biodiversity patterns in a species‐rich region. In this case, we identified distinct zones of both evolutionary relicts and fronts. But which are more important? ancient refugia of ‘Gondwanan’ heritage or zones of recent intercontinental intermixing and evolutionary potential? Although there is no clear answer to this question, it brings up important points for the field of conservation. Simply choosing a site or a series of sites that maximize PD in this scenario misses a fundamental story and the underlying processes that created the PD to begin with. In this case, maximizing PD without historical interpretation would deprioritize Australia's largest remaining tropical rain forest refugium (the Atherton‐Bellenden Ker uplands Fig. ). This area has low genus diversity but has remained stable enough through time to serve as an ark for the survival of relict lineages from Australia's ancient past. Being one of the last major fragments of continuous rain forest in Australia with a history stretching back over 40 million years, an approach which deprioritizes this area in favour of areas with high diversity due to the incursion of foreign lineages from Southeast Asia in comparatively more recent times may result in undesirable or unintended conservation outcomes. Alternatively, if the origin of lineages that dispersed into Australia comparatively recently was equally ancient or if the habitat of origin became extinct, these incursion zones (high PD residuals, Fig. ) could be viewed as a different type of refugium. Indeed, with the unabated habitat loss occurring in Southeast Asia and other tropical regions, this logic may prove to be a useful argument in other case studies around the world. Fortunately, the majority of extant rain forest in north‐east Queensland is well protected, thus we have explored these complex questions without the ‘agony of choice’. Many other species‐rich tropical regions of the world are not in the same position and may benefit from the methods proposed here to help make informed long‐term management decisions under a range of different conservation criteria. In Fig. , we outline a generalized summary of our approach so that it may be tailored and applied in other regions. Where generating new community phylogenies from sequence data is unfeasible, bioinformatics tools such as Phylomatic (Webb & Donoghue, ) can be used to generate trees suitable for estimating PD. As data on historical biogeography accumulates, the biogeographical origin of groups can be determined at a minimal cost with a species list and a dedicated literature review. This integrated approach may be used to verify or provide further support for other rain forest refugia inferred using traditional approaches, such as climate modelling and population genetics for individual species, and can rapidly advance knowledge on the evolutionary history of large species assemblages in other parts of the globe. Summary of our methodology combining phylogenetic diversity ( PD ) and historical biogeography data. Acknowledgements We thank the University of Adelaide, Australian Tropical Herbarium and Barcode of Life Data Systems molecular biology labs for helping generate the sequence data for this project.

Journal

Diversity and Distributions – Wiley

Published: Mar 1, 2015

There are no references for this article.