Phylogeny and biogeography of the carnivorous plant family Droseraceae with representative Drosera species from Northeast India

Survey, collection and taxon sampling

Insectivorous plant species in the genus Drosera were collected from different regions of Meghalaya, according to their present availability. The collected plants included two species of Drosera viz. D. peltata and D. burmannii (Figure 1 and Figure 2). Drosera burmanii Vahl. was collected from Jarain, Jaintia Hills District, Meghalaya (N 25°36ʹ, E 92°15ʹ) and Drosera peltata Sm. was collected from Cherrapunjee, East Khasi Hills District, Meghalaya (N 25°07ʹ, E 91°28ʹ). Identification of these insectivorous plants was carried out at the Botanical Survey of India (BSI), Eastern circle, Shillong, Meghalaya. Herbariums were prepared and submitted in BSI and Department of Botany, North-Eastern Hill University (NEHU), Shillong. Specimen voucher numbers (NEHU) and accession numbers (BSI) of Drosera burmanii are 11924 and 86843; and Drosera peltata are 11962 and 86840, respectively. We amplified the whole ITS and rbcL regions from all the above-mentioned plants for the proposed work. In addition, we collected GenBank data that included these markers from representative species belonging to the genus Nepenthes, Drosera, Aldrovanda, Dionaea and Sarracenia, along with their geographical distribution information (Table 1).

Figure 1. Drosera peltata.

(a) Plant in natural habitat, (b) close up view of leaf, (c) trapped insect, and (d) whole plant.

Figure 2. Drosera burmanii. 

(a) Plant in natural habitat, (b) close up view of leaf, (c) trapped insects (blue arrows), and (d) whole plant.

DNA extraction, PCR amplification and sequencing

For Drosera sp., the leaves and the stems were taken, washed thoroughly with water removing all the dirt and debris of insects, kept in 70% alcohol for a few minutes, dried, and then wrapped in aluminum foil and stored in liquid nitrogen for further use. Total genomic isolation from Drosera was carried out using DNeasy Plant Mini Kit (Qiagen, USA), according to the manufacturer’s instructions with minor modifications (combination of a borate extraction buffer with the DNA extraction kit, and a proteinase K treatment during extraction). The chosen markers were subjected to PCR amplification with desired forward and reverse primer pairs, as listed in Table 2. Reactions were performed in PCR tubes with a final volume of 100 µl. Each reaction mixture contained 4 µl of genomic DNA (30 ng/µl), 6 µl of dNTP* (2mM), 6 µl of 10X taq buffer B*, 4 µl of MgCl₂ (25mM)*, 0.8µl of taq polymerase (3 units/µl)*, 8µl of each primer pair (10 pm) (Metabion, Germany), final volume was made by adding sterilized Millipore water (*Bangalore Genei, India). DNA amplification was performed in an Applied Biosystems Gene-Amp PCR System 2700 programmed for 35 cycles (4 min at 94°C, 30 sec at 94°C, 1 min at 56°C , 40 sec at 72°C and 10 min at 72°C). Sequencing was carried out at Macrogen, Inc, Korea.

Phylogenetic analysis

Maximum likelihood. All ITS and rbcL sequences were first aligned using MUSCLE²¹ and subsequently concatenated using MESQUITE V3.03²². Highly variable sequence regions were excluded from analyses of the extended data set. Because initial separate calculations using noncoding spacer regions and coding matK sequences, yielded congruent but incompletely resolved topologies, respectively, both partitions were combined in all subsequent analyses. Maximum likelihood (ML) analyses were carried out using MEGA 7²³. To find the best substitution model for our analyses, ML fits of 24 different nucleotide substitution models were performed. Models with the lowest BIC scores (Bayesian Information Criterion) are considered to describe the best substitution pattern. For each model, AICc value (Akaike Information Criterion, corrected), Maximum Likelihood value (lnL), and the number of parameters (including branch lengths) were also computed. The following models were verified for this study: General Time Reversible (GTR); Hasegawa-Kishino-Yano (HKY); Tamura-Nei (TN93); Tamura 3-parameter (T92); Kimura 2-parameter (K2); Jukes-Cantor (JC).

For ITS and rbcL dataset, the evolutionary history was inferred by using the ML method based on the Tamura-Nei model²⁴. Sequence information for the aligned dataset pertaining to total number of sites (excluding sites with gaps/missing data), sites with alignment gaps or missing data, invariable (monomorphic) sites, G+C content, parsimony informative sites, number of haplotypes (h), haplotype gene diversity (Hd), Nucleotide diversity per site (Pi), average number of nucleotide differences (k) were computed and are shown in Table 3. The evolutionary history of the taxa analyzed are represented from the bootstrap consensus tree from 500 replicates of the original dataset²⁵. Branches with less than 55% bootstrap replicates were collapsed. Initial tree(s) for the heuristic search were obtained by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach. Evolutionary analyses were conducted in MEGA 7. The ML tree was further used in divergence time analysis.

Bayesian inference. The concatenated dataset from the previous analysis (in nexus format) was further analyzed with MrBayes v.3.1.2²⁶. The model of best fit was (GTR + Γ + I), determined with Modeltest. Posterior probabilities were estimated by sampling trees from the posterior probability distribution using the Metropolis-coupled Markov chain-Monte Carlo approach (MCMCMC) implemented in MrBayes²⁶, using default priors. The temperature of the heated chain was set to 0.2. Three times, four chains were run for 1 million generations. Consensus trees, clade posterior probabilities, and mean branch lengths were computed relying on the trees sampled every 10 generations after the burn-in and the tree was viewed in FigTree v1.3.1. A split decomposition median network graph tree was generated in SplitsTree4²⁷, with the variable positions in the same dataset.

TaxonGap analysis: Testing of markers for barcode suitability

A visualization analysis tool, TaxonGap 2.4.1²⁸, was used to illustrate the sequence divergences within and between species of the candidate markers from the ITS, rbcL and matK regions representing the family Droseraceae.

Secondary structure prediction and analysis

Consensus structures of ITS2 regions for the three genera in the Droseraceae family were predicted using LocARNA²⁹ from Freiburg RNA tools server, which outputs a multiple alignment together with a consensus structure. For the folding, a very realistic energy model for RNAs was used that features RIBOSUM-like similarity scoring and realistic gap cost. The high performance of LocARNA was mainly achieved by employing base pair probabilities during the alignment procedure.

Fossil data calibration, age estimation and ancestral area reconstruction

A Time-tree was generated using the RelTime method in MEGA 7^23,30. Divergence times for all branching points in the user-supplied topology for estimating phylogenetic history and divergence times of Droseraceae were calculated using the ML method based on the Tamura-Nei model²⁴. The estimated log likelihood value of the topology shown is -26876.8062. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 2.0910)). The rate variation model allowed for some sites to be evolutionarily invariable ([+I], 0.0000% sites). The tree is drawn to scale, with branch lengths measured in the relative number of substitutions per site. There were a total of 530 positions in the final dataset.

Fossils of Droseraceae pollens from the Eocene (55-38 MYA)³¹, present day Aldrovanda vesiculosa and its ancestral species (pollen fossil records) dating to Miocene in South Ural, Eocene in Kazakhstan and Belgium and Paleocene in East Germany were considered as internal calibration points while drawing the Time-tree phylogeny. The first real Drosera pollens have appeared in sediments from Miocene (25-5 MYA)³¹. Based on these studies and fossil data information on Droseraceae³¹, the following constraints were applied with a normal prior distribution that spanned the full range of nodal age estimates: the most recent common ancestor (MRCA) of Droseraceae (divergence between Drosera and Aldrovando species) was set with a minimum and maximum divergence time to 38 and 55, respectively; the MRCA of Drosera species was set to 22-5 MYA.

For biogeographic inference, Bayesian Binary MCMC (BBM) and Statistical Dispersal-Vicariance Analysis (S-DIVA) methods were employed in which biogeographic reconstructions were averaged over a sample of highly probable Bayesian trees³². In S-DIVA the occurrence of an ancestral range at a node was computed using all alternative reconstruction frequencies generated by the DIVA algorithm for each tree in the data set. To account for both phylogenetic and ancestral states uncertainty S-DIVA was utilized for an entire posterior distribution of trees. Different geographic areas of endemism for the carnivorous plants meant for this study consistent with the present distribution with both outgroup and in-group sampling are outlined in Table 1.

Phylogenetic analysis

The rbcL gene and ITS regions were separately aligned and ML and parsimonious trees were used for the phylogenetic analyses. The nucleotide sequences were aligned without any insertions or deletions. A total of 52 accessions from the Droseraceae family, including N. khasiana (Nepenthaceae) and S. flava (Sarraceniaceae) (used as out-groups), were considered for revisiting the Drosera phylogeny (Table 1). A separate concatenated dataset of ITS and rbcL were taken for Bayesian phylogeny reconstruction. The consensus core secondary structures of ITS regions for Drosera species as per their geographical distribution were drawn and shown in Figure 3. The ML tree was further used for time divergence studies.

Figure 3. Maximum likelihood (ITS + rbcL) tree of family Droseraceae.

Nepenthes khasiana and Sarracenia flava were taken as out groups. Substitutions per site between taxon groups according to the model of best fit (GTR + Γ + I) determined with Modeltest in MEGA 7. Consensus ITS2 core secondary structures drawn in LOCARNA for representative Drosera species based on geographical distribution are shown next to their respective phylogenetic groupings.

Major clades within Drosera determined via phylogenetic analysis were subjected to relative rate tests using ML estimates of substitutions per site between taxon groups according to the model of best fit (GTR + Γ + I), determined with Modeltest in MEGA 7²³. Relative rates were intended to provide evidence of whether certain longer branches within Drosera were the result of rate acceleration of individual species. Models with the lowest BIC scores were considered suitable for the analysis. For each model, AICc value, ML value and the number of parameters (including branch lengths) were also presented (Table 4). Evolutionary rates among sites and their non-uniformity were modeled by using a discrete Gamma distribution (+G) with 5 rate categories. Estimates of gamma shape parameter and/or the estimated fraction of invariant sites are shown (Table 4). For each model assumed or estimated values of transition/transversion bias (R) are shown and are followed by nucleotide frequencies (f) and rates of base substitutions (r) for each nucleotide pair. Sum of r-values is made equal to 1 for each model.

For estimating ML values, a user-specified topology was used. The analysis involved 52 nucleotide sequences. Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated. There were a total of 530 positions in the final dataset. All analyses supported the monophyly of Droseraceae. Aldrovanda and Dionaea species cladded well, and emerged as a sister branch to the genus Drosera. N. khasiana and S. flava emerged as outgroups. The Bayesian tree and the split tree also congrued with the results of ML tree, though there were slight changes in the overall topology of Drosera species, with very high support (100 percent bootstrap support [BS] values; 1.0 Bayesian posterior probability). Drosera species (D. burmannii and D. peltata) from Meghalaya grouped separately and emerged to be evolutionarily primitive (Figure 4 and Figure 5).

Figure 4. Bayesian phylogenetic tree of family Droseraceae based on concatenated whole ITS and rbcL markers.

Colour codes represent phylogenetic nodes with percent probability scores. Bayesian phylogeny reconstruction obtained by posterior probabilities for the nodes in the ML tree. GTR evolutionary model was implemented in MrBayes 3.2 with four chains during 50,000 generations and trees were sampled every 100 generations.

Figure 5. Split tree graph using Neigbourjoin NET in SplitTree for the family Droseraceae based on concatenated whole ITS and rbcL markers.

TaxonGap analysis

A DNA barcode marker is judged by its resolving power to discriminate species at generic and infrageneric levels. The intra- and inter-specific sequence divergence amongst the candidate markers chosen for the present study showed a comparative pictorial barcode gap in form of taxon-plots for the marker candidates (ITS, matK, rbcL) for species representing the family Droseraceae. The results are summarized in Figure 6. For each species, sequence similarity of the same gene within the same species was high; therefore, the relevant intra-specific variation (shown as dark grey bars) was low. TaxonGap plots have the discriminatory power to gauge better marker barcodes when phylogenetic trees for multiple genes need to be compared. Moreover, it uses the same scaling for depicting distance values based on individual biomarkers, thus making it straightforward to evaluate multiple genes rather than the need for comparing separate gene trees drawn for each of the taxonomic units. In the present study, it emerged that a combination of different markers (rbcL+ITS+matK) would render them better discriminatory power for identifying species in the carnivorous plant diversity.

Figure 6. Comparisons of intra and inter-specific variations among ITS, matK and rbcL genes of the carnivorous plant family Droseraceae.

The grey and black bars represent the intra- and inter-specific variations, respectively. The thin, black lines denote the smallest inter-specific variation. Names appearing next to the dark bars denote the closest species to that listed on the left.

Molecular divergence time estimates

Several genera of Droseraceae have been blessed with fossil pollens. A single record called Fischeripollis from European Mid Miocene has been assigned to Dionaea³³. Fossil seed information and even leaves have contributed to the understanding of Drosera origin on the geological timescale. Drosera pollens have been recorded since Lower Miocene from New Zealand³⁴. Several findings of tertiary pollen in the Mid Miocene from Europe have been assigned to either Drosera (Droserapollis) or Nepenthes (Droseridites)³⁵. The molecular data calibration with cue from previous studies is in congruence with the fossil record information of Droseraceae pollen, thus testifying a wide distribution of the progenitors of Aldrovanda in the Droseraceae family since Late Cretaceous (Figure 7).

Figure 7. Rel TimeTree chronogram for the combined dataset of rbcL and ITS showing divergence time estimation of Droseraceae.

Divergence times for all branching points in the user-supplied calibrated topology were calculated using the ML method based on the Tamura-Nei model. A discrete Gamma distribution was used to model evolutionary rate differences among sites [5 categories (+G, parameter = 2.0910)]. Numbers at nodes are median ages in million of years (Ma) with two internal calibration points. Evolutionary analyses were conducted in MEGA7.

Ancestral area reconstruction

The RASP tree (Figure 8) indicated the phylogenetic roots of Dionaea and Aldrovanda to originate in the northern hemisphere, while Drosera species would have most probably had an Australasian origin. Apparently all palaeoendemics (D. meristocaulis, D. burmannii, D. arcturi) are scattered throughout the southern hemisphere and also in tropical America. In this respect, the extant fossil record, i.e. European Miocene fossils, are somewhat noteworthy. A very old age (Cretaceous) can therefore be hypothesized for the whole family, dating back to stages of tectonic development when South America, Africa, and Australia were in closer proximity compared to present day geographical barriers. From the recent studies, it emerges that Australia is perhaps the secondary center of diversity of the genus Drosera, and most of the Drosera descendants can be assumed to have originated here.

Figure 8. Ancestral area reconstruction and world map distribution of species belonging to the family Droseraceae using RASP (S-DIVA and BBM).

Alternative ancestral ranges of nodes (with frequency of occurrence) are shown in pie chart form. Bootstrap support values/Bayesian credibility values posterior probabilities (50% and higher) are indicated near the pie chart in the Bayesian tree. Colour key to possible ancestral ranges at different nodes represent other ancestral ranges. Nepenthes khasiana and Sarracenia flava were taken as out-groups.