The ‘IAPAR-123’ passion fruit (P. edulis) accession described in Carneiro et al. (2002) was used in the present study. The other complete cp genomes and cp gene sequences were downloaded from GenBank. A list of species and GenBank accession numbers are provided in the Supplementary Table S1.

The large-insert BAC library of P. edulis was constructed and maintained at the French Plant Genomic Resources Centre³. It was previously accessed using the BES approach and comparative genome mapping, and two clones (Pe69Q4G9 and Pe85Q4F4) were found to match the Arabidopsis thaliana chloroplast genome (Santos et al., 2014). In the present study, these clones were selected so that their entire inserts could be sequenced. The DNA was then isolated using the Nucleobond Xtra Midi Plus kit (Macherey-Nagel, Düren, Germany) following the manufacturer’s instructions, after growing the bacterial clones on LB medium (100 mL) containing chloramphenicol as the selective marker (12.5 μg/mL).

Around 1.5 μg of each individual cpDNA were pooled together with other P. edulis BAC inserts (12 in total) for the construction of an SMRT library using the standard Pacific Biosciences (San Francisco, CA, USA) preparation protocol for 10 kb libraries. The pool was then sequenced in one SMRT Cell using the P4 polymerase in combination with C2 chemistry, following the manufacturer’s standard operating procedures and using the Pacific Biosciences PacBio RS II platform. Sequencing was performed by GATC Biotech⁴.

The reads were assembled following a hierarchical genome assembly process HGAP workflow (Chin et al., 2013), and using the SMRT^® Analysis (v2.2.0) software suite for HGAP implementation. Reads were first aligned by the PacBio long read aligner or BLASR (Chaisson and Tesler, 2012) against the complete genome of Escherichia coli strain K12 substrain DH10B (GenBank: CP000948.1). The E. coli reads, as well as low quality reads (minimum read length of 500 bp and minimum read quality of 0.80), were removed from the data set. Filtered reads were then preassembled to yield long, highly accurate sequences. To perform this step, smallest and longest reads were separated from each other to correct read errors by mapping single-pass reads to longest reads (seed reads), which represent the longest portion of the read length distribution. Next, the sequences were filtered against vector (BAC) sequences, and the Celera assembler was used to assemble data and obtain draft assemblies. The last step of the HGAP workflow is performed in order to significantly reduce the remaining InDel and base substitution errors in the draft assembly. The Quiver algorithm was used for this purpose. It is a quality-aware consensus algorithm that uses rich quality scores (QV scores) embedded in Pacific Biosciences’ bas.h5 files. Once the polished assembly was obtained, each BAC sequence was individualized by matching its paired-end sequences to the assembled sequences using BLAST. Read coverage was assessed by aligning the raw reads on the assembled sequences with BLASR.

The sequences obtained from both inserts (Pe69Q4G9 and Pe85Q4F4) were aligned using ClustalX software (Larkin et al., 2007) to obtain a single contig. Specific primers were designed at the sequence ends of this contig in order to find out whether the circular chloroplast genome was complete. PCR reactions were then performed using a 9700 thermal cycler (Applied Biosystems, Foster City, CA, USA) in reaction mixtures containing 20 ng template DNA (P. edulis accession ‘IAPAR-123’), 1× buffer, 1 mM MgCl₂, 0.2 mM of each dNTP, 0.3 μM of the forward and reverse primers, 1.2 U Go Taq Flex DNA polymerase (Promega, Madison, WI, USA), and ultra-pure water to bring the final volume up to 20 μL. The thermal profile for amplification was: 95°C for 5 min, 35 cycles at 95°C for 40 s, 55°C for 40 s and 72°C for 1 min, followed by a final 8 min incubation at 72°C. The amplified fragments were checked on 1% (w/v) agarose gel with a 100 bp molecular size standard Invitrogen (Carlsbad, CA, USA). The PCR product was purified using the Wizard^® SV Gel kit and PCR Clean-Up System (Promega), and then used as a template for the sequencing reaction based on the Sanger method. It was subjected to capillary electrophoresis in the ABI Prism 3100 sequencer (Applied Biosystems). The sequence of the PCR product was aligned with the single contig to obtain the complete sequence of the cp genome.

We performed two phylogenomic studies. The first was restricted to the Malpighiales based on the available entire chloroplast genomes of members of the four families that compose this order (22 species in total, Supplementary Table S1). In the second, a set of chloroplast genes was used to infer the relationships within the Fabids (42 species in total, Supplementary Table S1).

Entire cloroplastidial genomes of Malpighiales representing the families Passifloraceae, Salicaceae, Euphorbiaceae, and Chrysobalanaceae were used as the ingroup in phylogenomic comparisons, and the cp genome of V. vinifera (Vitaceae, Vitales) was used as the outgroup in order to root the phylogenetic tree. The data set consisting of 23 taxa was aligned at nucleotide level in server-based program MAFFT version 7.221, using the FFT-NS-2 algorithm with default settings. To generate the alignment, inverted sequences detected in the cp genomes of P. edulis and H. brasiliensis were reversed and the respective positions adjusted. The raw alignment was manually corrected in BioEdit (Hall, 1999) and further processed in GBLOCKS 0.91b (Castresana, 2000) in order to remove regions containing gap positions, with a minimum block length of five, and maximum number of contiguous non-conserved positions of eight. The resulting alignment was analyzed in jModelTest software version 2.1.7 (Darriba et al., 2012) to determine the optimal model of molecular evolution and gamma rate heterogeneity using the AIC.

Maximum likelihood analysis was performed using PAUP version 4.0a146 (Swofford, 2002), based on the transversional (TVM) substitution model, gamma distribution of rate heterogeneity with five discrete categories (+G). To estimate the level of support for the ML topology, bootstrap analysis was performed on 1,000 replicates. BI was run in MrBayes, version 3.1.2 (Ronquist and Huelsenbeck, 2003) but based on the GTR +G model, the second model chosen by jModelTest taking the AIC values into account. The Markov chain Monte Carlo (MCMC) algorithm was run for 5,000,000 generations, sampling one tree every 100 generations. The first 25% of trees were discarded as burn-in to estimate the values of posterior probabilities. Convergence diagnostics were monitored on the basis of an average standard deviation of split frequencies below 0.01, potential scale reduction factor (PSRF) values close to 1.0, and ESS values above 200. Trees were visualized using FigTree version 1.4.0.

A set of 43 nucleotide sequences of homologous protein-coding chloroplast genes from 42 species representing all orders of the Fabids clade (Rosales, Fagales, Cucurbitales, Fabales, Malpighiales, Celastrales, Oxalidales, Zygophillales) were used as the ingroup in the phylogenomic comparisons. V. vinifera (Vitaceae, Vitales) was chosen to serve as outgroup to produce a rooted tree. A list of the chloroplast gene sequence sources is provided in Supplementary Table S1.

First, each protein-coding gene sequence was aligned using ClustalW with default settings and the raw alignments manually corrected in program BioEdit (Hall, 1999) and further processed in GBLOCKS (Castresana, 2000), excluding gap positions from the data set, with a minimum block length of five, and a maximum number of contiguous non-conserved positions of eight. Next, all individual filtered alignments were concatenated into a single alignment matrix. Both conserved and variable nucleotide positions in the alignment matrix were analyzed in MEGA6 (Tamura et al., 2013).

jModelTest software version 2.1.7 (Darriba et al., 2012) was used to determine the optimal model of molecular evolution and gamma rate heterogeneity based on AIC. Both ML and BI were used to infer the phylogenomic relationships within the Fabids clade. ML analysis was performed using RAxML version 8.2.4 (Stamatakis, 2014). The GTR model of nucleotide substitution was selected for ML analysis, taking into account the gamma distribution of rate heterogeneity with five discrete categories (+G). To estimate the support of the ML topology, a bootstrap analysis was performed on 1,000 replicates. BI was run on MrBayes software, version 3.1.2 (Ronquist and Huelsenbeck, 2003) with the GTR +G model. The MCMC algorithm ran for 5,000,000 generations, sampling one tree every 100 generations. The first 25% of trees were discarded as burn-in to estimate the values of posterior probabilities. Convergence diagnostics were monitored on the basis an average standard deviation of split frequencies below 0.01, PSRF values close to 1.0, and ESS values above 200. Trees were visualized using FigTree version 1.4.0.

Following extraction of reads containing only chloroplast genome sequence data and subsequent error correction, 2,340 PacBio RS reads from the clone insert Pe69Q4G9 were recovered, ranging from 500 to 22,458 bp, and containing a total of 94,052 bp assembled into a contig. The average depth of coverage of the consensus sequence was 96× and the average GC content was 37%. The final quality of the assembly corresponded to a nominal QV of 48.48 (approximately 99.999% accuracy). Similarly, 4,972 reads from the clone insert Pe85Q4F4 were recovered, ranging from 500 to 26,035 bp, and containing a total of 91,155 bp assembled into a contig. The average depth of coverage of the consensus sequence was 172× and the average GC content was 38%. The final quality of the assembly corresponded to a QV of 48.57.

It is worth noting the usefulness of the HGAP workflow as a solution for successfully resolving long repeat regions, as already pointed out by Chin et al. (2013). The high quality sequence data generated by the PacBio RS II Platform, added to its capability to assemble long reads, allowed us to obtain a single contig for each clone insert and then the complete cp genome of P. edulis. Both contigs were aligned and merged into a single long contig of 151,362 bp with an overlapping region of 33,848 bp. The amplification reaction using the primer pair designed to anneal at the extremities of the large contig and genomic DNA from P. edulis accession ‘IAPAR-123’ as a template produced an amplicon of 325 bp, which was subsequently sequenced and aligned with the long contig sequence. Thus, a 44-nucleotide sequence was obtained and added for closing the gap to produce the complete sequence of the circular molecule of 151,406 bp. The chloroplast genome sequence of P. edulis has been submitted to the NCBI GenBank and has received the accession number KX290855.

The chloroplast genomes of Potentilla micrantha (Ferrarini et al., 2013) and Aconitum barbatum var. puberulum (Chen et al., 2015) were entirely sequenced using the PacBio RS II Platform, and both studies described the advantages of using this method. For instance, Ferrarini et al. (2013) highlight the generation of a single contig covering the entire cp genome of P. micrantha, whereas Illumina HiSeq2000 sequencing data have lower genome coverage (some regions have zero or very low coverage) and the resulting assembly consisted of seven contigs. Similarly, Chen et al. (2015) emphasize the low level of errors (~0.0027%) in PacBio reads, since after applying the necessary filter, the result is an error-corrected consensus read with a higher intra-molecular accuracy. Interestingly, even in the absence of any other biological information on the target species, the authors stated that it took less than half an hour to finish the genome assembly step, confirming that long read lengths are superior, especially for de novo assemblies.

Chloroplast genome annotation resulted in the identification of 105 unique genes, including 30 tRNAs, 4 rRNAs and 71 protein coding genes (Table 1). The molecule has a typical quadripartite structure characterized by the presence of two copies of an IRA and IRB each of 26,154 bp (34.6% in total) separating a SSC region of 13,378 bp (8.8%) and a LSC region of 85,720 bp (56.6%). Genes arising from duplication events created perfect IRs, each containing the same 16 genes, totalling 120 genes in the complete chloroplast DNA molecule (Figure 1).

Protein coding genes constitute 37.3% of the chloroplast genome (75 genes, totalling 56,428 bp), in addition to tRNA and rRNA coding genes that constitute 1.8% (37 genes totalling 2,801 bp) and 6.0% (eight genes totalling 9,048 bp) respectively of the chloroplast genome. Introns constitute 10.7% (totalling 16,155 bp), while intergenic regions and pseudogenes form the remaining 44.2% (totalling 66,974 bp). The amount of GC was estimated at 37%.

tRNA genes are spread throughout the genome and encode the 20 amino acids incorporated into proteins. One is present in the SSC region, 21 in the LSC, and 7 in the IRs, which also have 5 protein-coding genes and 4 rRNA genes. A total of 18,799 codons represent the coding capacity of the 71 protein-coding genes (Table 2). The most frequent amino acid was found to be leucine (2,026 codons corresponding to 10.8% of total DNA) and the least frequent cysteine (218 codons corresponding to 1.2% of total DNA).

Fifteen unique genes (nine protein- and six tRNA-coding genes) have introns, and two introns were found in only one gene, ycf3. The largest intron occurs in the trnK-UUU gene (2,524 bp) in which the matK gene (1,506 bp) is inserted, and the smallest occurs in the rps12 gene (537 bp). The number of introns identified in the cp genome of P. edulis is similar to that of other species. For instance, in S. purpurea (Wu, 2015) there are 17 intron-containing genes, three of them containing two introns, including ycf3. Similar figures have been found for R. communis (Rivarola et al., 2011), whereas in H. brasiliensis (Tangphatsornruang et al., 2011), 23 introns are inserted in 22 cp genes.

The matK gene was found in the intronic region of the trnK-UUU gene in P. trichocarpa (Tuskan et al., 2006), H. brasiliensis (Tangphatsornruang et al., 2011), and M. esculenta (Daniell et al., 2008). This gene encodes the unique maturase found in the plastid genomes of land plant species, but it does not contain the reverse transcriptase domain and therefore is not able to promote intron mobility. The clpP gene has no intron in P. edulis, although it does have two introns in the following conordinal-related species: P. trichocarpa (Tuskan et al., 2006), S. purpurea (Wu, 2015), H. brasiliensis (Tangphatsornruang et al., 2011), M. esculenta (Daniell et al., 2008), J. curcas (Asif et al., 2010), and R. communis (Rivarola et al., 2011). Similar intron losses have already been reported in legumes, almost exclusively in the clade known as the IR-lacking clade (IRLC) (Jansen et al., 2008). Coincidently, the clpP gene is located at the beginning of first the inversion found in the P. edulis cp genome (Figure 1), which may lead to intron losses.

Group I and II introns have been derived from mobile genetic elements, which explains why they are lost or gained in the evolution of chloroplast genomes (Barkan, 2004). Both cyanobacteria and algae, as well as land plant species, share a single group I intron in the trnL-UAA gene, which is therefore considered the more ancestral (Simon et al., 2003). Usually, the atpF gene shows a conserved group II intron; however, it is absent in the atpF gene of P. edulis. This loss was also reported in other species of Passiflora (Jansen et al., 2007). Daniell et al. (2008) have suggested an association between C-to-T substitutions (at nt position 92) and the loss of the intron in M. esculenta and in other atpF gene sequences of Malphigiales, implying that recombination between an edited mRNA and the atpF gene may be a possible mechanism for this intron loss.

Two protein-coding genes show alternative initiation codons. The ‘GTG’ triplet was found in the rps19 and ndhD genes. This triplet was also found in the rps19 gene of J. curcas (Asif et al., 2010), H. brasiliensis (Tangphatsornruang et al., 2011), and Cynara cardunculus (Curci et al., 2015). The rps12 gene is a trans-spliced gene consisting of three exons: the first exon (5′-rps12) is located in the LSC region, far from the other two located in the IRs. This kind of organization was observed in P. trichocarpa, S. purpurea, H. brasiliensis, M. esculenta, J. curcas, and R. communis. Interestingly the infA gene, which codes for translation initiation factor 1, and the rps16 gene, which codes for a S16 ribosomal protein are absent in the P. edulis cp genome, as previously reported in other species of Passiflora (Jansen et al., 2007). These genes were also lost or are non-functional in related species M. esculenta (Daniell et al., 2008), J. curcas (Asif et al., 2010) and in seven species of Salicaceae (Wu, 2015). It is worth noting that Passifloraceae are known to have chloroplast gene and intron losses (Hansen et al., 2006; Jansen et al., 2007; Daniell et al., 2008).

Eight pseudogenes were identified, five in the IRs and three in the LSC region. The nucleotide sequences of the rpl22, ycf1, and ycf2 pseudogenes are smaller in length compared to those of the functional genes identified in other genomes. We found a truncated portion of the rpl22 pseudogene that has also been identified in Passiflora biflora, P. quadrangularis, and P. cirrhiflora (Jansen et al., 2011), but no studies have been performed to demonstrate whether this partial copy is functional or whether there is an rpl22 functional copy in the nucleus. There is evidence that this gene was exported to the nucleus in some Rosids (Castanea, Prunus, and Theobroma) and Fagaceae, and it is assumed that the transfer resulted from two independent events. An additional event may have occurred in Passiflora (Jansen et al., 2011).

Only one of the two copies of ycf1 is a pseudogene in the H. brasiliensis, M. esculenta, J. curcas, and C. cardunculus cp genomes, as most of their sequences have been lost. Interestingly, in P. edulis both copies lost parts of their sequences, representing non-functional ycf1 genes (Figure 1). On the other hand, the rpl20, rps7, ycf15, and ycf68 pseudogenes have 5, 1, 3, and 3 internal stop codons, while the accD gene was found to have repetitive elements at the 5′ end.

A comparison of the P. edulis cp genome to those of the 22 species representing the four families of the Malpighiales order showed that each of the families (Passifloraceae, Salicaceae, Euphorbiaceae, and Chrysobalanaceae) constitutes a monophyletic group (Supplementary Table S1). The multiple-plastome alignment of the complete cp genome sequence was 208,695 bp in length. Removal of alignment columns containing gaps reduced the alignment length to 118,724 bp (nucleotide positions).

Both maximum likelihod and BI-based methods resulted in a similar tree (Figure 5), and all species were clustered in Malpighiales, in conformity with the APG II system. The nodes resulting from ML analysis were consistently supported in the tree. The bootstrap values at most of the nodes (15 of 20 nodes) were ≥85%. In the BI the value of the average standard deviation of split frequencies was 0.0001 after 5,000,000 generations. The PSRF values were close to 1.0 (ranging between 0.999 and 1.000), and the ESS values were above 200 (ranging between 4,909 and 29,555) for all the parameters. All these values indicated that the analysis has reached convergence. The nodes resulting from BI analysis were also highly supported (PP = 1 for 17 nodes). P. edulis has been placed near the Salicaceae (Salix and Populus species) with strong support (BS = 100%; PP = 1), but distant from members of the Chrysobalanaceae, that in turn are closer to the Euphorbiaceae members studied herein.

In Malpighiales, a phylogeny based on plastid, mitochondrial and nuclear gene regions (Wurdack and Davis, 2009) and a recent phylogenomic approach based on a set of 82 chloroplast genes (Xi et al., 2012) indicated the relationship between Passifloraceae and Salicaceae. According to Wurdack and Davis (2009), for instance, most of its members share parietal placentation. Our results confirm this association and the potential of complete chloroplast genome sequences to infer evolutionary relationships. Malpighiales is an order that underwent rapid basal radiation. Therefore the use of large molecular data sets to identify several phylogenetically informative sites is an important step in deducing its course of evolution.

The Salicaceae species were clustered into two clades, both highly supported (BS = 100%; PP = 1). The most related were S. purpurea and S. suchowensis (BS = 100%; PP = 1), but positioned at some distance from those species is S. interior (BS = 100%, PP = 1). Populus species were split into two groups, the first incorporating Populus alba, P. tremula, and P. euphratica (BS = 78%; PP = 1), and the second P. balsamifera, P. trichocarpa, and P. fremontii (BS = 100%; PP = 1). The placement of Salix and Populus species is similar to that recently shown in the phylogenetic tree built from the cp DNA sequences of seven Salicaceae species (Wu, 2015).

With reference to Euphorbiaceae, M. esculenta was placed near to H. brasiliensis (BS = 100%; PP = 1), similar to the findings published by Xi et al. (2012) after examining a set of 82 chloroplast genes. Our findings placed J. curcas as the species most distant from the other Euphorbiaceae (BS = 100%; PP = 1), occupying an onward position relative to R. communis (BS = 82%; PP = 1). The placement of these species has weak support in two previous topologies based on 62 (Su et al., 2014) and 60 chloroplast genes (Kong and Yang, 2016). Remarkably, our study resolved the node positions in the Euphorbiaceae with strong support in both analyses (ML and BI). These findings also demonstrate the significance of using complete cp genomes in phylogeny reconstructions, since not only the genes but also non-coding sequences are examined. These sequences are informative, particularly at low taxonomic levels, due to their rapid evolution. In contrast, protein-coding genes evolve at a relatively slow rate (Folk et al., 2015).

Regarding the Chrysobalanaceae, the placement of Licania heteromorpha near to the species of Hirtella with strong support (BS = 100%; PP = 1) makes this genus paraphyletic. Parinari campestris was the most basal species in this family and this position is consistent with the findings of Bardon et al. (2013), based on six chloroplast markers and one nuclear marker. Moreover, the placement of the eight Chrysobalanaceae species with strong support is consistent with the phylogenetic hypothesis proposed by Malé et al. (2014) analyzing the complete plastid DNAs of eight species of this tropical family.

It is important to emphasize that our phylogenomic study is the first to take into account all the complete chloroplast genomes of the Malpighiales taxa available in databanks. Malpighiales includes 40 families, which are poorly represented in sequence databanks. Consequently, a limited taxonomic sampling of the order was used in the analysis and this may have contributed to the high support of the nodes in the phylogenetic trees obtained. Therefore, it is not possible to reach any definitive conclusions concerning the monophyly or position of the families within the order without examining the chloroplast genome sequences of members of the other families. However, our study lends support not only to the utility of complete chloroplast genomes to infer phylogenetic relationships, but as well as the enormous utility of long sequence reads from PacBio sequencing and assembly for generating very highly accurate biological information.

Subsequently, the complete nucleotide sequences of 43 chloroplast protein-coding genes of 42 species (including P. edulis) that belong to the eight orders of the Fabids clade were compared to reconstruct the phylogenomic relationships based on the GTR +G model (Supplementary Table S1), with V. vinifera (Vitales: Rosidae) as the outgroup. After alignment filtering and concatenating the nucleotide sequences into a single matrix, a total of 31,193 nucleotide positions were analyzed. Conservation analysis of positions was then performed and it is estimated that 61.65% correspond to conserved sites and 38.35% to variable sites. To summarize, these rates were estimated in Mega, a program that identifies a site as constant only if at least two sequences contain unambiguous nucleotides. In contrast, a variable site contains at least two types of nucleotide.

The method of analysis (ML or BI) had no substantial effect on the resulting trees. Their topologies were found to be highly similar (Figure 6). The accuracy of the inferred species’ phylogeny is strongly supported by the stability of the main clades obtained using different phylogenetic methods (ML and BI), and the best-scoring phylogenomic tree is shown in Figure 6. The nodes resulting from ML analysis were consistently supported in the tree (41 nodes in total). The bootstrap values at most of the nodes were higher than 85% (34 of 41), reaching 100% in some cases (28 nodes). In the BI the value of the average standard deviation of split frequencies was 0.0003 after 5,000,000 generations. The PSRF values were close to 1.0 (range between 0.999 and 1.000), and the ESS values were above 200 (range between 2,564 and 23,705) for all the parameters. All these values indicate that the analysis has reached convergence. The nodes resulting from BI analysis were also well-supported (PP = 1 for 31 nodes).

The Fabids consist of eight orders, which are all represented in this study and have been found to be monophyletic. For instance, Bulnesia arborea was the most distant species of the ingroup, thus placing the order Zygophyllales basal to the Fabids clade. This position of Zygophyllales has been previously reported (Wang et al., 2009).

Moreover, we were able to recognize two major monophyletic subgroups with strong support (BS = 100%; PP = 1). The first subgroup includes the orders Rosales, Fagales, Cucurbitales, and Fabales, known as the nitrogen-fixing clade, and the second consists of the orders Malpighiales, Celastrales, Oxalidales known as the COM clade. There is a trend in the literature toward classifying these subgroups as monophyletic (Wang et al., 2009; Su et al., 2014; Kong and Yang, 2016). Although there is morphological divergence among the species comprising the nitrogen-fixing clade, this capacity is shared only by angiosperms that belong to the orders Rosales, Fagales, Cucurbitales, and Fabales (Svistoonoff et al., 2013).

The members of the order Oxalidales (Averrhoa carambola and Oxalis latifolia) were placed as sister to a grade comprising the members of the order Celastrales (Elaeodendron orientale and Euonymus americanus) with strong support (BS = 100%; PP = 1). Thus, our studies indicate that Malpighiales shares more homologies with Celastrales than with Oxalidales, based on both methods used to infer species relationships within the Fabids. Previously, Hilu et al. (2003) reported a feasible relationship among the orders of the COM clade, based on sequence alignments of matK, a cp gene, but with weak support using maximum parsimony (BS = 60%). According to Sun et al. (2015), positioning the orders within the COM clade remains a great challenge. However, our results confirm other findings, for instance placing Celastrales nearer to Malpighiales than Oxalidales (Zhang and Simmons, 2006; Bell et al., 2010). The same scenario was described by Matthews and Endress (2005) observing floral structures and describing their implications for systematics in Celastrales. The use of a large set of genes is important for defining genomic relationships; otherwise the results could be interpreted as relating to the evolutionary course of particular genes.

It is worth noting that the Passifloraceae species comprise a distinct monophyletic group. In morphological terms, the monophyly of this family is supported by the occurrence of a flower corona (Judd et al., 2008). Additionally, our data suggest that P. edulis and P. ciliata are in close proximity. Both belong to the subgenus Passiflora. In taxonomic terms P. edulis is included in the supersection Passiflora, section Passiflora, and series Passiflora; P. ciliate is in the supersection Stipulata, section Dysosmia. However, P. biflora belongs to the subgenus Decaloba in supersection Decaloba, section Decaloba (Ulmer and MacDougal, 2004; Hansen et al., 2006). Passifloraceae appeared as a sister group to Salicaceae (BS = 100%; PP = 1) with a distant relationship to Chrysobalanaceae (BS = 52%; PP = 0.72) compared to Euphorbiaceae, confirming other findings based on four chloroplast, six mitochondrial and three nuclear gene sequences (Wurdack and Davis, 2009) and 82 chloroplast genes shared by 58 Malpighiales species (Xi et al., 2012).

However, we were unable to resolve the relationship between J. curcas and R. communis (BS ≤ 50%; PP = ≤ 0.50) in our work, confirming previous studies based on 60 chloroplast genes and the maximum likelihood method for estimating phylogeny inferences (Kong and Yang, 2016). Although the current phylogenetic reconstruction of Malpighiales is much improved compared to earlier versions, it is still incomplete, because of the limited taxonomic sampling available, and further phylogenetic and morphological studies are needed, focused especially on relations within Euphorbiaceae, by conducting, for instance, a phylogenomic analysis based on entire chloroplast genome sequences.