Young fresh leaves from the insectivorous plants (Drosera and Nepenthes species) were washed under a stream of water to remove dead insects and plant debris. The surface of the leaves was sensitized with a cotton swab to activate the release of hydrolytic enzymes. Leaves were thoroughly rinsed under running water to wash off the mucilage. Subsequently the leaves were gently blotted into a paper and then 1g of leaves were ground in liquid nitrogen. Further it was incubated for 1 h at 60° C in a 10 mL extraction buffer (1.4M sodium chloride (NaCl); 100 mM Tris hydrochloride (HCl); 20mM ehtylene diamine tetraacetic acid (EDTA); 2% cetyltrimethyl ammonium bromide, CTAB), followed by 30 min and 2 h incubation steps at 37°C with RNase A (0.4 mg) and proteinase K (0.4 mg), respectively. DNA was extracted with a ratio of phenol: chlorophorm: isoamylic acid (25: 24: 1) and chlorophorm: isoamylic acid (24: 1) and precipitated with isopropanol. The pellet was washed with ethanol (70%), dried, dissolved in 50 μL of water and cleaned with a QIAquick PCR purification kit (QIAGEN, Clifton Hill, Victoria, Australia,) and stored at −20°C. The chosen markers were subjected to PCR for amplification with desired forward and reverse primer pairs as listed in Table 4. Amplification for the target regions were carried out using the thermal cycle Applied Biosystems Gene-Amp PCR System 2700. A QIAquick PCR purification kit (QIAGEN) was used for cleansing the PCR mixtures as per the manufacturer's instructions. Sequencing was conducted at Macrogen. Inc., Korea using the primers listed in Table 4.

Nucleotide sequences representing Nepenthaceae, Droseraceae, Drosophyllaceae, Dioncophyllaceae and the non-carnivorous plant family Ancistrocladaceae, were first aligned using MUSCLE (Edgar, 2004) and ClustalX (Larkin et al., 2007) and then manually refined by eye. The alignments for each of the markers (mtDNA, cpDNA, and rDNA) were concatenated using MESQUITE V3.03 (Maddison and Maddison, 2016). Additionally for taxon gap analysis we included species from the carnivorous plant families Sarraceniaceae and Roridulaceae. ML bootstrap percentages were estimated from 1000 rapid bootstrapping pseudo-replicates and Bayesian posterior probabilities were obtained from BEAST (Drummond and Rambaut, 2007). The combined dataset was partitioned by locus and analyzed using the General Time Reversible (GTR) model. Rate heterogeneity was modeled by assuming that some sites are invariable and that the rate of evolution at other sites approximates a discrete gamma distribution [GTR+I+Γ]). It was pre-determined to be the best fitting model based on a likelihood ratio test for the concatenated data, as well as for each of the individual partitions. ML and Bayes trees were inferred from concatenated dataset of ITS and matK regions only, as data from other markers for the three genomic regions (mtDNA, cpDNA, rDNA) were not sufficiently available in public domains. For ITS and matK dataset the evolutionary history was inferred by using the ML method based on the Tamura-Nei model (Tamura and Nei, 1993). 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. The bootstrap consensus tree inferred from 500 replicates (Felsenstein, 1985) is taken to represent the evolutionary history of the taxa analyzed. Branches corresponding to partitions reproduced in less than 55% bootstrap replicates were collapsed and the percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (500 replicates) are shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach, and then selecting the topology with superior log likelihood value. Evolutionary analyses were conducted in MEGA 7 (Kumar et al., 2016). The ML tree was further used in divergence time analysis.

The concatenated dataset from the previous analysis (in nexus format) was further analyzed with MrBayes v.3.1.2 (Ronquist and Huelsenbeck, 2003). Two parallel Bayesian analyses with four chains each and partitioned by DNA region were run for 50 million generations. Tree construction was carried out using general time reversible substitution model (GTR) with substitution rates estimated by MrBayes. Metropolis-Coupled Markov Chain Monte Carlo (MCMCMC) sampling was performed with two incrementally heated chains that were combinatorially run for 100,000 generations. Coalescence of substitution rate and rate model parameters were also examined. Average standard deviation of split frequencies was carried out and generations were added until the standard deviation value was below 0.01. Posterior Probabilities indicated clade support. A cladogram with the posterior probabilities for each split and a phylogram with mean branch lengths was generated and subsequently read by FigTree v1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/). The node labels were color coded with a consensus 50 major rule based on percent probabilities.

A visualization analysis tool, TaxonGap 2.4.1 (Slabbinck et al., 2008), was used to illustrate the sequence divergences within and between species of the candidate markers from the cp, mt and rDNA regions representing the four carnivorous plant families (Nepenthaceae, Droseraceae, Sarraceniaceae, and Roridulaceae).

A novel way of representing secondary structures for ease of manual inspection was carried out via “rnaplot (RNA2ndStruct)” function in Matlab R2012a with specified format values, “Tree Diagram” for ITS2 sequences representing the four carnivorous families mentioned in previous sections.

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 (Yu et al., 2015).

A Timetree was generated using the RelTime method (Tamura et al., 2012). Divergence times for all branching points in the user-supplied topology for estimating phylogenetic history and divergence times of Nepenthaceae were calculated using the Maximum Likelihood me based on the Tamura-Nei model (Tamura and Nei, 1993). We combined the rDNA (ITS) and cp (matK) datasets. We implemented a relaxed molecular clock (uncorrelated lognormal) (Drummond et al., 2006). A fossil for Sarracenia is previously reported (Li, 2005). Earlier studies on 21 fossil constraints (series of secondary age constraints) from an angiosperm-wide analysis (Bell et al., 2010; Ellison et al., 2012) have also thrown light on evolutionary age of Caryophyllales. Based on these studies and fossil data information on Droseraceae and Aldrovanda, 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 Aldrovanda species) was set with a minimum and maximum divergence time to 38 and 55 respectively; and the MRCA of Drosera species was set to 22-5 MYA. Drosophyllum lusitanicum species were taken as outgroup for the ML timetree.

Our aligned matK and ITS datasets included 1,564 and 979 characters, respectively. Within Caryophyllales the four different carnivorous families—Nepenthaceae, Dioncophyllaceae, Droseraceae, Drosophyllaceae and one non-carnivorous family Ancistrocladaceae are clearly rooted with monophyletic groupings with very high support (100 percent bootstrap support values; 1.0 Bayesian posterior probability) (Figure 1). A common ancestor, presumably one with flypaper traps, links the families Nepenthaceae, Ancistrocladaceae, and Dioncophyllaceae in a large clade to the family Droseraceae. The tree topology was analyzed by using the ML method based on the Tamura-Nei model. The log likelihood of the tree is (−7261.35). The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The analysis involved 41 nucleotide sequences. All positions containing gaps and missing data were eliminated. There were a total of 1,389 positions in the final dataset. The same tree was incorporated into the RelTime tree analysis (Figure 2). The ML and Bayes phylogenies were largely congruent while N. khasiana (the only native species identified in the Indian sub-continent) and N. pervillei came out to be the primitive one, as they were placed basal to the overall topology of Nepenthes species. The Bayesian tree was drawn and color-coded with percentage probability of node values (Figure 1). Corresponding to the tree topology there were basal polytomies viz. the N. distillatoria, endemic to Sri Lanka, N. pervillei from Seychelles, N. truncata from Philippines and N. ampullaria from Malayan Archipelago. A separate distinct branch is indicated by N. khasiana, endemic to North-east India (Meghalaya). Sister group to N. khasiana with all other Nepenthes taxa that are predominantly distributed on the islands of the Malay Archipelago in a range from Malaysia to New Guinea extending as Far East as northern Australia and New Caledonia can be seen (Figure 1). These internal clades represent three main evolutionary lineages, each well supported by high probability values.

Suitability of a DNA barcode marker is judged by the intra- and inter-specific sequence divergence amongst the markers chosen for study. Comparisons of the six marker candidates (ITS, matK, rbcL, atpB, rps4, and coxI) for species representing the families Nepenthaceae, Droseraceae, Sarraceniaceae, and Roridulaceae were performed. Sarracenia and Roridula species were taken as outgroups. The results are summarized in Figure 3. For each species, sequence similarity of each gene within the same species was high; therefore, the relevant intra-specific variation (shown as gray bars) was low. For the matK and ITS regions, the inter-specific variations were very similar, and are apparently higher than those for other marker regions. There was a clear gap between the maximum intra-specific variations and the minimum inter-specific variations in Nepenthes and Drosera species. From the TaxonGap analysis it emerges that ITS, matK, rbcL, and atpB are reliable markers that can be further investigated for barcode analysis with potential species discriminatory power in carnivorous group of plants. CoxI emerges to be a poor candidate when it comes to barcoding and species identification in plants.

Tree diagrams represented a rooted tree, where leaf nodes correspond to unpaired residues and internal nodes (except the root) correspond to base pairs. They exhibited the branching pattern with blue colored dots joined with each other as paired regions and red colored dots as leaves depicting unpaired regions. These diagrams gave an explicit picture of ITS2 secondary structures for a visual comparison of representative species belonging to the four carnivorous plant families (Figures 4A–C).

Our Timetree age divergence estimates (Figures 1, 2) indicate that the family Nepenthaceae diversified from its sister families Ancistrocladaceae-Dioncophyllaceae in the Early Eocene (~38.42) Million years ago (Mya). The three families Nepenthaceae, Ancistrocladaceae and Dioncophyllaceae diversified from Droseraceae in the Order Nepenthales suggesting a monophyletic origin of botanical carnivory in the Late Cretaceous (~84.07) Mya. The present day Nepenthes species diversified in the Miocene and early Miocene (15-5) Mya. Species rich Malesian clades evolved during the early Miocene (6-4) Mya (Figure 2).

It could be inferred from the ancestral area reconstructions that Nepenthes khasiana from Meghalaya, N. distillatoria from Srilanka and N. pervillei from Seychelles are isolated species but geographically serve as intermediate link between the western species and the Indo-Malayan taxa (Figure 5). The transitional position of N. khasiana from Northeast India with respect to its morphological features shows closer affinity to species in the Indo-Malayan taxa well supported by the sister group relationships with the subclades of Asian stem-group. Nepenthaceae probably originated in Asia and Africa with species in the crown-group widespread in the Malayan Archipelago. The most recent common ancestor of N. khasiana (Meghalaya) and N. distillatoria (Sri Lanka) was likely present in Africa. Our molecular data studies indicate colonization of Nepenthes species in South East Asia from the Indian subcontinent. The relict N. khasiana, endemic to Northeast India marks this migration to the Malay Peninsula. Subsequent processes of colonization started in the Indo-China and the Indo-Malayan Archipelago where three major lineages can be outlined. Initially the Malay Peninsula and the main land of Indo-China were probably occupied. Subsequent colonization in Sumatra from the Malay Peninsula marked a process of pompom radiation of Nepenthes species that we see today. If we assume this as the center of diversity, taxa were dispersed along the Sumatran track to Borneo where a similar process of divergence took place. Islands of the Malay Archipelago arose gradually during the Tertiary, making it much easier for the Nepenthes diversity to expand in the mainland than the islands accounting for the wide range distribution of some taxa (N. ampullaria, N. gracilis, N. mirabilis). Studies indicate spread of N. mirabilis to the Philippines, Sulawesi, Borneo, New Guinea and northern Australia especially at times of low sea level in the Quaternary (Figures 2, 5).

Evolution of Nepenthes and its spectacular pitchers is no less than a fairytale that has undergone a remarkable journey in the geological time scale to reach the present day diversity. Danser (1928) was probably the first one to study Nepenthes species throwing light on interrelationships of the numerous described species (summarized in Jebb and Cheek, 1997), which merits further extensive research. Molecular marker studies have confirmed a basal position of the western species (from Madagascar to Sri Lanka and Northeast India Khasi Hills), monophyly of a group of New Guinean species and of a clade containing the moderately derived species (Regiae, Borneo, Philippines, Sulawesi, New Guinea). Nepenthes colonized Southeast Asia from an ancient Indian stock and subsequently moved to Indo-China and the Malay Archipelago (Meimberg et al., 2001).

Our phylogenetic analysis of cp and nu genes is an extension of the work initiated by Meimberg et al. (2001, 2006) which lacked computational work on ancestral area reconstruction and age estimation of the monotypic family Nepenthaceae though they have discussed the biogeographic features of the Nepenthaceae family in detail. Our work corroborates most of their findings with a fully resolved phylogeny that supports the monophyletic nature of carnivorous plant families Nepenthaceae and Droseraceae in the order Caryophyllales. The biogeographic analyses reveal that the evolution of botanical carnivory in the order Nepenthales would have occurred in the Northern Tethys which is supported by fossil pollen records from the European Eocene and a Gondwanaland origin at a time when the Indian plate was alienated from Madagascar. Our estimates of divergence times within and among clades (Figure 5) also provide support for the vicariance hypothesis in previous studies explaining the biogeographic history of the family.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.