Orchidaceae is considered the largest family of flowering plants, with over 25,000 species (Dressler, 1990; Christenhusz & Byng, 2016). The family probably dates back to the Late Cretaceous, as indicated by fossil-calibrated molecular phylogenies (Gustafsson, Verola & Antonelli, 2010; Ramírez et al., 2007, 2011). Polystachya Hook. is an orchid genus containing 240 species, with most species found in Africa (Dressler, 1993). A total of 13 species are reported from the Neotropical region (Mytnik-Ejsmont, 2011), but this number may increase when considering the endemic species from Brazil that were not accounted for in Mytnik-Ejsmont (2011) or were considered synonymous (Barros et al., 2010).
Recent studies have shown a number of peculiar cytological, evolutionary and biogeographic aspects of Polystachya. The genus has diploid and polyploid species; the latter recently formed in the Neotropics and Madagascar (Rupp et al., 2010; Russell et al., 2010b). Unlike most genera of Orchidaceae, Polystachya has a wide geographical distribution range (Pridgeon et al., 2005; Fig. 1), having species that are Pantropical or have a transatlantic distribution. On the other hand, the Neotropics presents a high level of endemism. Brazil, as an example, has 12 species of which 10 are endemic (Barros et al., 2010). In addition, there is evidence of reticulate evolution in the genus and hybridization with independent origins (Russell et al., 2010a).
The monophyly of the genus has been reported in the latest studies (Russell et al., 2010b; Mytnik-Ejsmont, 2011), which contrasts starkly with the low level of monophyly observed in the taxonomic sections described within the genus. Those 15 sections (Kraenzlin, 1926; Summerhayes, 1942, 1947 apud Russell et al., 2011; Brenan, 1954; Cribb, 1978) are based on morphological characters and have been useful for field identification and inventories, but do not find support as natural groupings in the molecular studies currently available (Russell et al., 2010b; Mytnik-Ejsmont, 2011). According to those molecular studies, all sections are polyphyletic or paraphyletic, except sect. Isochiloides (Russell et al., 2010b).
Section Polystachya has been described as comprising 32 species worldwide and is the only section with species of pantropical distribution (Mytnik-Ejsmont, 2011). However, according to molecular analyses, some species of this section appear to be more related to species of other sections (Russell et al., 2010b; Mytnik-Ejsmont, 2011). These studies highlight the need for new infrageneric divisions based on robust molecular evidence. Russell et al. (2010b), using chloroplast markers, defined five different clades that could be used as the basis for a revised classification of new sections within the genus. Clade III (sensu Russell et al., 2010b) includes species from five different sections (sects. Polystachya, Eurychilae, Caulescentes, Superpositae, Polychaete) and is divided into distinct subclades of morphologically diverse plants. These species are Pantropical (such as P. concreta), Neotropical (such as P. foliosa), Malagasy endemics (such as P. henrici) or African (such as P. odorata). The relationships among Clade III species remain unresolved because specimens of P. concreta, P. foliosa, P. henrici and P. modesta form a large polytomy, due to low levels of divergence between sequences. In addition, several distinct subclades are unresolved at the base. The molecular phylogenetic studies produced so far have included about 35% of the recognized species within the genus, and used only a small number of nuclear (PgiC between exons 11 and 15, PhyC exon 1, Rpb2 intron 23 and ITS in Russell et al., 2010a; only ITS in Mytnik-Ejsmont, 2011) and chloroplast markers (rps16 intron, rps16 exon 2, rps16-trnK spacer, trnK intron excluding matK, matK and psbD-trnT spacer in Russell et al., 2010b; rps16-trnK, rps16 and rpl32-trnL in Mytnik-Ejsmont, 2011).
New methods of DNA sequencing as well as the development of more powerful algorithms are propelling the replacement of trees generated from one or a few genes to those constructed from hundreds of them (Edwards, 2009). The improvement of massively parallel sequencing techniques—or next generation sequencing (NGS)—has increased the amount of data available for biological research, whether the fully annotated reference genomes of species under study have been sequenced or not (Bräutigam & Gowik, 2010). However, despite its obvious potential, NGS technology is underused in most studies of plant systematics (Cronn et al., 2012; Carstens et al., 2013; Eaton & Ree, 2013), probably as a result of a prevailing focus on non-model organisms (which require de novo genomic sequencing and its inherent challenges), the need to sample many individuals per species and the absence of well-established protocols (McCormack et al., 2013). One method that increases the efficiency of NGS for non-model species compared to other genomic partitioning strategies is sequence capture (or hybridization-based enrichment), which is based on the prior selection of loci of interest (Lemmon & Lemmon, 2013). The main benefit of this technique is that the number of specific sequences obtained can be very high, which makes it an advantageous method compared to PCR-based approaches if the objective is to sequence several individuals and multiple loci. Furthermore, sequence capture when combined with NGS platforms, such as Illumina, also reduces the costs of the process (Lemmon, Emme & Lemmon, 2012). Sequence capture methods have successfully been used in other plant genera to generate large amounts of useful data for phylogenetic inference (Kamneva et al., 2017; Sousa et al., 2014; Stephens et al., 2015; Weitemier et al., 2014).
The necessity of a molecular phylogenetic framework for (and a morphological taxonomic revision of) Polystachya is clear. It requires a well-resolved phylogenetic hypothesis in order to clarify the relationships between species and also to redefine new infrageneric sections. In this paper, we explore the use of nearly complete plastid genomes (note that we use chloroplast and plastid interchangeably when referring to these genomes), obtained by sequence capture and massively parallel sequencing, to solve the phylogenetic inconsistencies found within Clade III of Polystachya (sensu Russell et al., 2010b). We also explore whether sequencing the entire chloroplast genome using NGS was worthwhile, compared to PCR and Sanger sequencing of a few fast-evolving loci. We hope that the results generated here can be extended to the rest of the genus and thus result in new interpretations of the evolutionary and biogeographic history of the group.
Materials and Methods
Sampling and DNA extraction
We sampled 19 species and 48 individuals (Table 1; Fig. 1), of which 15 were collected in different locations in continental Brazil and three were collected on Trindade Island in the South Atlantic. The DNA of the 15 Brazilian samples was extracted from 10 mg of tissue dried with silica gel and using the DNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA). DNA samples for the remaining 33 individuals were provided by the University of Vienna and the DNA bank of the Royal Botanic Gardens, Kew. To this sample of individuals, selected because they cluster in Clade III (as defined in Russell et al., 2010b), we added Polystachya tessellata Lindl., supposedly synonymous of P. concreta. We also included multiple samples of P. concreta, because previous studies have reported a lack of monophyly for this species and several synonymous species. Permits to collect were provided by the Ministério do Meio Ambiente (MMA), Instituto Chico Mendes de Conservação da Biodiversidade (ICMBio) and Sistema de Autorização e Informação em Biodiversidade (SISBIO), with registration number 29478-1.
|Species||Location||Voucher||ENA reference||Before genomic library assembly||After genomic library assembly|
|Concentration (ng/uL)||Purity 260/280||Volume||Concentration (ng/uL)||Purity 260/280|
|P. bicolor (=P. rosea)||–||Kew 25884||ERS2203588||20.1||1.51||50||31.7||1.93|
|P. concreta (1)||Brazil, Distrito Federal||N.L.Abreu 254||ERS2203559||39.2||1.92||75||17.2||1.92|
|P. concreta (2)||Brazil, Distrito Federal||N.L.Abreu 254||ERS2203560||33.6||1.96||75||31.9||1.9|
|P. concreta (3)||Brazil, Bahia State||N.L.Abreu 251||ERS2203564||18||1.78||70||29.9||1.97|
|P. concreta (4)||Brazil, Bahia State||N.L.Abreu 251||ERS2203565||8.5||1.52||130||38.3||1.9|
|P. concreta (5)||Cameroon||A. Russell 40 (YA)||ERS2203570||35||1.82||8||15||1.75|
|P. concreta (6)||Brazil||HBV ORCH 066004||ERS2203572||21.9||1.53||6||9.8||1.9|
|P. concreta (7)||Mauritius||HBV ORCH 07278||ERS2203575||26||1.78||13||32.4||1.85|
|P. concreta (8)||Réunion||HBV “Chase & Samuel 1”||ERS2203576||34||1.82||13||30||1.85|
|P. concreta (9)||Comoros||HBV ORCH 07417||ERS2203577||23.7||1.83||13||30.4||1.86|
|*P. concreta (10)||Madagascar||Fischer&Sieder FS3210 (WU)||–||21||1.81||13||9.4||1.63|
|*P. concreta (11)||Madagascar||Kew 17854||–||70||1.65||67||9.6||1.6|
|*P. cornigera||Madagascar||Fischer&Sieder FS3208 (WU)||–||22.2||1.76||11||7.7||1.58|
|P. dolichophylla||Cameroon||Kew 25886||ERS2203589||108.9||1.05||55||30.3||1.89|
|P. estrellensis (1)||Brazil, Minas Gerais State||N.L.Abreu 255||ERS2203551||29.8||1.81||65||33.7||1.9|
|P. estrellensis (2)||Brazil, Minas Gerais State||N.L.Abreu 255||ERS2203552||17.8||2.06||90||26||1.54|
|P. estrellensis (3)||Brazil, Minas Gerais State||N.L.Abreu 255||ERS2203553||11.9||1.82||75||35||1.86|
|P. estrellensis (4)||Brazil, Espírito Santo State||N.L.Abreu 253||ERS2203554||12.8||1.77||65||18.7||1.84|
|P. estrellensis (5)||Brazil, Espírito Santo State||N.L.Abreu 253||ERS2203555||22.4||1.87||70||21||1.96|
|P. estrellensis (6)||Brazil, Espírito Santo State||N.L.Abreu 253||ERS2203556||28.2||1.79||70||26.3||1.96|
|P. estrellensis (7)||Brazil, Bahia State||N.L.Abreu 252||ERS2203557||25.4||1.91||70||29.3||1.98|
|P. estrellensis (8)||Brazil, Bahia State||N.L.Abreu 252||ERS2203558||14.6||1.79||70||22.4||2.02|
|P. estrellensis (9)||Brazil, São Paula State||N.L.Abreu 256||ERS2203561||27.8||1.61||90||20.4||2.05|
|P. estrellensis (10)||Brazil, São Paula State||N.L.Abreu 256||ERS2203562||23.8||1.81||75||23.1||1.93|
|P. estrellensis (11)||Brazil, São Paula State||N.L.Abreu 256||ERS2203563||22.3||1.68||90||28||1.97|
|P. estrellensis (12)||Brazil, Trindade’s Island||–||ERS2203566||10.1||1.64||130||30||1.88|
|P. estrellensis (13)||Brazil, Trindade’s Island||–||ERS2203567||12.5||1.64||130||30||1.83|
|P. estrellensis (14)||Brazil, Trindade’s Island||–||ERS2203568||10.3||1.48||130||34.7||1.84|
|P. eurychila||Kenya||Kew 17963||ERS2203586||207.6||0.94||85||13.7||1.98|
|P. foliosa (1)||Dominica||Kew 25887||ERS2203590||14.7||1.58||56||31.8||1.96|
|P. foliosa (2)||Venezuela||HBV ORCH 07082||ERS2203574||30.2||1.73||8||18.6||1.79|
|P. golungensis||Kenya||Kew 17966||ERS2203587||104.7||1.14||115||26.1||1.88|
|*P. henrici||Madagascar||Kew 17856||–||22||1.52||56||7.7||2.04|
|P. humbertii (1)||Madagascar||Fischer&Sieder FS2079 (WU)||ERS2203578||116.3||1.82||13||29.4||1.77|
|P. humbertii (2)||Madagascar||Fischer&Sieder FS3017 (WU)||ERS2203580||35.1||1.93||11||34.2||1.89|
|P. melanantha||Kenya||Kew 17954||ERS2203584||207.2||0.94||50||12.2||1.8|
|P. modesta||–||HBV ORCH 05165||ERS2203573||56.8||1.35||13||31.7||1.93|
|P. odorata (1)||Nigeria||Kew 17857||ERS2203581||33.7||1.56||58||46.8||1.8|
|P. odorata (2)||Cameroon||A. Russell 42 (YA)||ERS2203571||36.9||1.91||13||26.2||1.91|
|P. oreocharis (1)||Madagascar||Fischer&Sieder FS2082 (WU)||ERS2203579||58||1.74||13||30.2||1.87|
|*P. oreocharis (2)||Madagascar||Fischer&Sieder FS3152 (WU)||–||20.2||1.82||13||10.5||1.8|
|*P. paniculata (2)||Cameroon||L. Pearce 27 (YA)||–||56.1||1.4||8||8.1||1.98|
|*P. pinicola||Barzil||HBV ORCH 06606||–||31||1.11||7||9.2||1.98|
|P. steudneri||Kenya||Kew 17956||ERS2203585||10.6||1.45||58||33.5||1.87|
|P. tessellata (1) (=P.concreta)||Madagascar||Kew 17859||ERS2203582||176.2||1.1||58||32.7||1.89|
|P. tessellata (2) (=P.concreta)||Madagascar||Kew 17860||ERS2203583||216.1||1.1||49||33.3||1.93|
|*P. tsinjoarivensis (1)||Madagascar||Fischer&Sieder FS3209 (WU)||–||19.3||1.83||13||9.3||1.79|
|P. tsinjoarivensis (2)||Madagascar||HBV FS4182||ERS2203569||22.2||1.75||13||17.2||1.64|
Species analyzed; location of the collection and voucher; DNA concentration and purity before and after genomic library assembly.
Polystachya bicolor Rolfe and P. melanantha Schltr. were chosen as possible more distant relatives of the species in focus, in order to provide additional context for the phylogenetic inference. All studies conducted to date resolve P. melanantha as an outgroup with respect to the Clade III species. Polystachya bicolor has already been treated as a synonym of P. rosea (Mytnik-Ejsmont, 2011), with an uncertain position in the phylogenetic trees generated thus far, being sometimes closely related to P. concreta and other associated species (Mytnik-Ejsmont, 2011) and sometimes closely related to species from other clades (Russell et al., 2010b; Mytnik-Ejsmont, 2011).
Probe design for DNA capture
We used complete chloroplast genome of Phalaenopsis aphrodite subsp. formosana (NC_007499.1) (Chang et al., 2006) as the reference for the design of capture probes, because there is no completely sequenced chloroplast genome of a Polystachya. According to molecular analyses, Polystachya and Phalaenopsis belong to different sub-tribes but are closely related within the Vandeae tribe (Van Den Berg et al., 2005; Górniak, Paun & Chase, 2010; Freudenstein & Chase, 2015). The use of a quite distantly related species is made possible by the DNA capture kit (MYcroarray, Ann Arbor, MI, USA), which is able to support differences larger than 5% between probe sequences and target sequences (Li et al., 2013). The complete sequence of the chloroplast genome of Phalaenopsis aphrodite subsp. formosana was divided into blocks of 360 bp. Every second block was used as the template for probe design; the probes consisting of 120 bp sequences, three to each block without overlap. Given the genomic DNA fragment sizes of between 300 and 400 bp (below) and that fragments can contain complementary sequence anywhere on their length to attach to a probe, fragments can contain up to 200–300 bp of genomic sequence into the flanking regions beyond the probes, or cover the probes with little extent into the flanking sequence, or somewhere in between. In this way, captured sequences produce a series of tiled overlapping sequence for high quality genomic assembly. Additionally, fragments with a base repeated more than seven times in a row were avoided to reduce the capture of repetitive sequences present in many places in the genome. Finally, the reference sequences in blocks used for probe design totaled 63,720 bp and were brought together into a single FASTA file and sent to (MYcroarray, Ann Arbor, MI, USA) to produce the probes.
Sonication and genomic library preparation
Extracted DNA was randomly fragmented by sonication using a Covaris S220 instrument (Covaris, Woburn, MA, USA), in order to evenly cover the full genome. Adapters were incorporated into the fragmented DNA using NEXTflex™ DNA Sequencing Kit and NEXflex™ Barcodes kit (BIOO Scientific, Austin, TX, USA). Uniquely indexed adapters were used for each sample. We selected fragments between 300 and 400 bp using Agencourt AMPure XP magnetic beads kit (Beckman Coulter, Brea, CA, USA). The genomic library was amplified following the program: 98 °C for 2 min; 14 cycles (98 °C for 30 s; 65 °C for 30 s; 72 °C for 60 s); 72 °C for 4 min. The products were purified using a QIAquick PCR Purification Kit (Qiagen, Valencia, CA, USA). The genomic DNA concentrations before and after sonication and the amplification of the library were measured in a NanoDrop 2000c instrument (Thermo Fisher Scientific, Waltham, MA, USA) (Table 1) to ensure that the final concentration exceeded 400 ng/μL.
Enrichment and sequencing
Before the enrichment, equimolar amounts (400 ng/μL) of each amplified library were pooled into six reactions, each one containing eight indexed samples. The enrichment method involves the selection of genomic regions and capture of DNA samples before sequencing (Mamanova et al., 2010). The enrichment was performed with MYBaits target enrichment system (MYcroarray, Ann Arbor, MI, USA), following the manufacturer’s instructions. The probes were recovered using Dynabeads® MyOne™ Streptavidin C1 (Invitrogen Dynal AS, Oslo, Norway).
To increase DNA concentration, 14 cycles of PCR were performed for each hybridization reaction using Herculase II Fusion DNA Polymerase (Agilent, Waldbronn, Germany) and the following program: 98 °C for 30 s; 14 cycles (98 °C for 20 s; 60 °C for 30 s; 72 °C for 60 s), 72 °C for 5 min. Sequencing was performed on the Illumina MiSeq platform (San Diego, CA, USA) by the Genomics Core Facility (University of Gothenburg, Sweden).
Illumina reads were processed using the program CLC assembly cell (CLC Bio, Aarhus, Denmark). Firstly, the Illumina adapter sequences were removed and low-quality reads were excluded. Reads were then mapped against the reference sequence used for probe design (P. aphrodite). Consensus sequences generated for each sample were converted into FASTA format using the SAMTools software (Li et al., 2009) using the mpileup tool with reference sequence option, allowing for the inclusion of indels in the consensus sequences. These sequences were used as a new individual reference sequences for each sample in a second round of mapping. Final consensus sequences were generated using mpileup, without the reference sequence option, to avoid erroneous base calling in low read-depth portions of the read alignment. Sequence alignment was performed using the auto strategy in MAFFT—Multiple Sequence Alignment Software Version 7 (Katoh & Standley, 2013) and later manually refined using Geneious Pro (Biomatters Ltd., Auckland, New Zealand). In the last step we aligned the sequenced samples with the Phalaenopsis aphrodite subsp. formosana (NC_007499.1) chloroplast genome to obtain the sequenced region annotation.
Alignment and phylogenetic analysis
Micro-structural features of chloroplast non-coding sequences can have a profound influence on the multiple sequence alignment, and hence also the phylogeny. Hairpins (short stem-loop structures in single stranded DNA or RNA), for example, can create sites that allow small inversions to occur at a high enough frequency that homoplasious inversions can be observed among sequences from closely related species (Kelchner & Wendel, 1996). Sometimes the inverted sequence is not so short and can disrupt phylogenetic analysis, leading to strongly supported but spurious groupings (Joly et al., 2010). Non-coding sequences, such as group II chloroplast introns, contain many such stem-loop structures (Kelchner, 2002).
We examined the non-coding sequences in our alignment for inverted (reverse complemented) sequences and tested for their effect on the phylogenetic inference. This was done by excluding all but one character of the inversion (to down-weight the inversion to a single event) and rerunning the analysis. The selected character to represent the inversion was arbitrarily chosen. This was done to avoid recoding the inversions as indel characters and creating a new, small partition (with only eight characters) that would have required many additional parameters, in comparison to our approach.
The alignment process can sometimes be confronted with small regions that are difficult to align, probably most often due to overlapping indel events. We identified several such regions and excluded them using the nexus block commands from the Bayesian analysis. Poor alignment excluded sites: 54459-54465, 55220-55224, 55484-55489, 67162-67167, 67426-67431, 68186-68192, 118838-118843, 119266-119364, 119493-119726, 120012-120019, 120867-120872, 121295-121393, 121522-121755, 122041-122048. Inverted loop-associated excluded sites: 74668-74669, 88120-88128, 90474-90481, 94193-94205, 94696-94699, 97241-97245, 104712-104716, 106631-106635.
The aligned sequences were partitioned based on the chloroplast annotation of 116 functional genes, seven pseudogenes and two partitions that concatenated untranslated regions—one partition containing introns with secondary structure and the other concatenating all intergenic sequences. The files that contained the introns were submitted to trimAl v1.2 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009) to test different gap deletion settings.
Faster region assessment
We used the sequences from one sample on GenBank (FS1045) of Polystachya cultriformis and examined two published markers, psbD-trnT (870 length aligned to our samples) and matK (1,521 length aligned to our samples). We compared these sequences pairwise to one of our samples, P. estrelensis8, to check which was the faster-evolving region. We then ran a MrBayes v.3.2.4 (Ronquist et al., 2012) analysis on the faster region, as a representative of a fast part of the chloroplast genome (fast cpDNA hereafter), with one copy of our dataset trimmed to this region alone plus the GenBank sequence.
We then ran successive analyses, using additional copies of the dataset interleaved in the same file, to explore the increase in support with an increase of characters evolving under the same model. This was to discover how much of the fast cpDNA data would be needed to achieve high support on most nodes in the phylogeny (i.e., among species but not necessarily within species). An important assumption we made at this stage was that the single fast region would contain mutations spread across most/all branches of the phylogeny. If this assumption was true, then a single region could carry changes representative of the entire phylogenetic history that we were exploring. The assumption is essentially one of i.i.d. (independently and identically distributed sites)—in that the sites would be representative of many unsampled sites and that double mutations would be rare—coupled with a sufficient dataset size to contain enough changes overall to reflect the history. Although the i.i.d. assumption is rarely true across sites, model-based analysis methods can cope with this, because the frequency of site patterns can be modelled by an i.i.d. process (Steel & Penny, 2000). So we are in effect mainly testing whether the original data set size was sufficient to carry changes reflective of the entire phylogeny under investigation.
Random sample from all non-coding regions
A single region copied many times proved ineffective in recovering most nodes with support (see Results). We therefore explored using random samples of characters without replacement from among all of the non-coding regions in our dataset to test how much data from faster regions would yield supported trees across most nodes. We expected this approach to be less subject to the limitations caused by the stochastic nature of mutations coupled with the limited size of any one region. By sampling across many regions (over 58 kb in this case), even those few characters that have changed on short branches might be sampled occasionally. In contrast, a single region, by chance, may simply not contain any characters changing on a specific short branch.
We sampled without replacement 4%, 8% and 16% of the non-coding data using delete-fraction jackknifing in the seqboot program v3.69 (from http://evolution.genetics.washington.edu/phylip.html), excluding the poorly aligned parts and with down-weighting of the inverted loops (by excluding all but one character of each loop), in 20 replicates each. The approximate average (and range) of posterior probabilities (PP) per clade was taken across the 20 replicates to get an indication of the likely support for selected clades that a non-coding dataset of these sizes would generate. These values were plotted on the whole alignment analysis to compare to the support received when using the whole dataset. Given that the largest dataset we used here (hereafter the “16% dataset,” or ∼9.2 kb) failed to recover support for all nodes found in the whole genome analysis (see Results), we did not end up analyzing the smaller replicates.
MrBayes v.3.2.4 analyses were used for phylogenetic inference. These analyses were run for five million generations (two million for the random sample replicates), using a mixed substitution model (plus gamma and invariant sites) to account for among-site rate variation. Priors on branch lengths were set to unconstrained: exponential (100) to minimize the chance of inferring incorrectly long branches (Marshall, 2010), otherwise with default settings. The paired runs were checked for convergence and high effective samples sizes in the MrBayes output and Tracer v.1.6 (Rambaut et al., 2014), respectively. Burn-in generations were removed by discarding 10% of the samples of parameters and trees, while summarizing in TreeAnnotator v.1.8 (Rambaut & Drummond, 2010) to ascertain clade PP. Trees were rooted using the Phalaenopsis sequence. Analyses using the character partitions were also done, returning nearly identical results to the analysis described above, so they are not reported further.
Our NGS approach allowed the capture of coding and non-coding regions throughout the chloroplast genome. We recovered approximately 132 kb, after the exclusion of gaps, representing 116 genes, seven pseudogenes, as well as regions with intergenic sequences and introns with secondary structure. Compared to the reference annotation, seven genes contained frameshifts that are usually associated with pseudogenization and corresponded to previously reported pseudogenized genes in orchids (Luo et al., 2014).
We excluded eight of the 48 samples due to the low quality of the sequencing results (Table 1). These eight samples showed lower DNA concentrations after the genomic library construction assembly, which may be the cause of low quality sequencing. The remaining 40 samples were submitted to the EMBL/ENA database under accession numbers ERS2203551–ERS2203590. The coding regions have 48,308 polymorphic sites (38.4%). Introns with secondary structure and regions with intergenic sequences have 21,264 (16.9%) and 56,226 (44.7%) polymorphic sites, respectively. The alignment of the concatenated data showed an unbalanced (but fairly typical) mean nucleotide composition of A = 29.9%, C = 19.9%, G = 19.4% e T = 30.8%.
In the non-coding part of alignment, we found evidence for eight putative small inversions (Table 2), based on the presence of inverted repeated motifs that could form stems at least 4 bp long. Stems of this length or longer are part of models of group II structures (Michel, Kazuhiko & Haruo, 1989; Toor, Hausner & Zimmerly, 2001; Kelchner, 2002) and are consistent with sequence patterns observed by one of us in the rpL16 intron sequences in other taxa (Pfeil et al., 2002).
|Loop||Alignment co-ordinates||Sequence (majority)||Sequence (minority)||Substitutions before/after inversion reverted||Sequences (minority)|
|A||74,660–74,678||TCTATCTA-GAA-TAGATAGA||TCTATCTA-TTC-TAGATAGA||3/0||bicolor, concreta7, concreta9, tesselata1, melanantha, Phalaenopsis|
|B||88,106–88,142||GGCCCAATCTTTC-TTTTTTTGAGGA-AAAGATTGGGCC||GGCCCAATCTTTC-TCCTCAAAAAAA-AAAGATTGGGCC||8/0||humbertii2, oreocharis1, tsinjoarivensis2, melanantha|
|E||94,679–94,717||GTATCTAAGGAAGATCC-AAAGG-GGATCTTCCTTAGATAC||GTATCTAAGGAAGATCC-CTTCT-GGATCTTCCTTAGATAC||5/2||estrelensis12, estrelensis13, estrelensis14, eurychila, concreta7, concreta8, golungensis, oreocharis1, melanantha, humbertii1, Phalaenopsis|
|F||97,230–97,258||AACGTCCAGTG-CCAAAGT-CACTGAATGGG||^CCCATTCAGTG-ACTTTGG-CACTGAATGGG||6/0||All 13 sequences of the estrelensis clade (cartooned in Figs. 2 and 4)|
|G||104,698–104,730||ATTCAATTCTTC-ATTCTTTCAA-GAAGAATGAAT||ATTCAATTCTTC-ATGAAAGAAA- GAAGAATGAAT||6/0||foliosa1, foliosa2, odorata1, odorata2, Phalaenopsis|
|H||106,618–106,650||ATTCATTCTTCAT-GAAAGA-ATGAAGAATTGAAT||ATTCATTCTTCAT-TCTTTC-ATGAAGAATTGAAT||6/0||foliosa1, foliosa2, odorata1, odorata2, Phalaenopsis|
The loops that appear to have undergone inversions are underlined. Bolded bases would be treated as phylogenetically informative if the alignment does not take into account the inversion. Bases not involved in stem formation are in subscript.
Faster region assessment
The pairwise identity between Polystachya concreta5 and P. concreta8 (whose common ancestor is relatively old and near the crown of Polystachya) for psbD-trnT was 98.7%. The pairwise identity for these samples for matK was 99.2%. PsbD-trnT was therefore used as the representative fast cpDNA region. The analysis with a single copy of this dataset (870 bp) yielded a MCC tree (Fig. 2A) with only five nodes with high support (>0.95 PP). Increasing the number of copies did not result in much improvement. The analysis using 16 copies of psbD-trnT (∼14 kb) produced a MCC tree (Fig. 2B) containing only eight highly supported nodes.
Random sample from all non-coding regions
The pairwise identity between Polystachya concreta5 and P. concreta8 was 98.7% across 58 236 bp of non-coding region contained in our alignment. This compares to 98.8% identity between the same samples across the coding regions in our alignment.
The replicate datasets that sampled 16% of the original non-coding alignment (excluding poorly aligned parts and down-weighting the inverted loops) failed to return all nodes/clades found in the whole genome analysis. Of 13 selected nodes found in the tree from the whole genome (four subtended by relatively long branches, four by medium length branches, and five by short branches), only five were found with high support across most or all replicates (i.e., at least 16 of 20 replicates had ≥0.95 PP). Three of the selected nodes instead had five or fewer replicates with high support (≥0.95 PP), but only one or no replicates that contained highly supported contradictory nodes (thus the support for the expected node was ≤0.05). Finally, five of the nodes had generally poor support among replicates (i.e., five or fewer replicates had ≥0.95 PP along with six or more replicates with ≤0.05 PP for these nodes). Sixteen percent of datasets (∼9.2 kb) were recovered from seven to 21 highly supported nodes among replicates (mean = 14.8), with more nodes recovered in 19 of 20 replicates than was the case with the larger repeated psbD-trnT dataset (∼14 kb and eight supported nodes). This character sampling strategy was probably more reflective of the underlying support for various nodes than using repeated copies of a single small dataset.
The mutually exclusive foliosa1/concreta2 versus foliosa1/foliosa2 clades (see below) were also examined in the 16% datasets. In the first case (foliosa1/concreta2), just four replicates contained this clade with high or moderate support (≥0.90 PP). The contradictory second grouping (foliosa1/foliosa2) was found with a similar level of support (≥0.90 PP) in only two replicates. The fact that both groupings could be recovered, with high support, in at least some replicates suggests that the original dataset contains the signal of both clades. A NeighborNet analysis (Fig. 3B) confirmed that a mixture of patterns exists in the original dataset involving foliosa1, foliosa2, concreta1 and concreta2.
Analyses with and without the inverted loops (the latter by down-weighting to a single character) returned almost identical trees. The results of only the latter analysis is presented in this section. The tree we recovered was able to resolve the phylogenetic relationships among the groups of the large clade selected for this study, with high support values on almost every node (Fig. 3). The tree was characterized by a large clade with relatively short branches containing only sequences from the Neotropics, with a grade of a few small clades and single sequences containing the remaining sequences (Fig. 3). The large clade contained 21 sequences from Brazil, Dominica and Venezuela, whereas the grade included 19 sequences from tropical central and eastern Africa, as well as Madagascar and the nearby islands (Fig. 3).
The grade recovered include a few geographically identifiable clades (Fig. 3). One of these, attaching fairly deeply within the crown, consists of four Malagasy sequences (Polystachya humbertii1, P. humbertii2, P. oreocharis and P. tsinjoarivensis2) that are sister to a Kenyan sequence (P. eurychila). Another clade comprises a Kenyan sequence (P. golungensis) and one from Reunion (P. concreta8). A third clade contains a pair of central African sequences, one from Cameroon (P. odorata2) and one from Nigeria (P. odorata1). A fourth clade contains sequences from central Africa (P. concreta5 from Cameroon), Madagascar (P. tesselata1), the Comoros (P. concreta9) and two sequences without certain provenance. Finally, a fifth pair of sequences were from samples collected from Mauritius (P. concreta7) and Madagascar (P. tesselata2). Lineages containing only a single sequence in this grade included samples from Kenya (P. melanantha and P. steudneri) and Cameroon (P. dolichophylla).
Sequences from the widely sampled and widely distributed P. concreta did not form a monophyletic group and occurred on different branches of the tree, separated by several well supported nodes (Fig. 3). Similarly, the two P. tessellata sequences from Madagascar did not form a clade. Polystachya estrellensis sequences form a clade with P. concreta sequences collected in Brazil. Although the sequences of P. estrellensis are thus paraphyletic, whether the taxon itself is paraphyletic cannot be established for certain here. The identification of this P. concreta sample could be wrong, given that the identification of these species is confused in Brazil and sometimes they are considered synonymous (see also Discussion).
With versus without loops
The down weighting of the inversions we identified (by excluding all but one character per inversion) resulted in a similar, but not identical, phylogenetic inference. The differences among the maximum clade credibility (MCC) trees involved P. foliosa1, P. foliosa2, P. concreta1, P. concreta2, P. concreta3 and P. concreta4. The analysis with the inversions included returned this tree (Figs. 3A and 4A):
There are several supported differences between these trees, with at least one corresponding to the way the inverted loops are weighted. When the entire loops are analyzed, P. foliosa1 and P. foliosa2 are supported as sisters, with these two sequences appearing to share two loop inversions (if this topology is correct; Fig. 4A). However, down-weighting the inversions produces a tree consistent instead with two independent inversions (Fig. 4 main panel).
Chloroplast genome sequence provides a robust phylogeny
In this work we used the nearly complete chloroplast sequences of 40 Polystachya samples to infer a robust plastid phylogeny. The dataset significantly increased the phylogenetic resolution within the genus. Thus, our results suggest that increasing the number of molecular markers has the potential to solve not only the relationships among species, but also to identify new Polystachya clades and define new sections. The delimitation of new sections will, however, depend upon the inclusion of more taxa than was done by this study—in other words a higher coverage of the genus. Below we highlight some of the clades recovered, their morphological and/or geographical characterization, and a comparison with previous studies.
Polystachya bicolor/rosea position contradicts Russell’s Clade III
Of the two species selected as possible distant relatives to provide more context in the phylogenetic inference, one of them, P. bicolor (=P. rosea), appears in a clade together with samples of P. concreta (from Cameroon and from the Comoros), P. tesselata (= P. concreta) and P. modesta. Not surprisingly, the clade that includes P. bicolor/rosea is deeply nested within the ingroup, thus contradicting the monophyletic Clade III presented by Russell et al. (2010b).
In prior studies, Polystachya bicolor/rosea has an uncertain position in the phylogenetic trees. In an analysis using plastid markers and Bayesian inference, this species appears in a large polytomy with P. concreta and other related species (Mytnik-Ejsmont, 2011), or related to species of other clades (Russell et al., 2010b) depending on the marker used. A phylogeny using nuclear data (ITS sequences) highlighted the lack of monophyly of this species (Mytnik-Ejsmont, 2011), which may be connected to the difficulty in identifying it. Polystachya bicolor/rosea is often mistaken for P. concreta, since differentiation between these is made by subtle differences in the shapes of leaves, and the size and color of the flowers. Unlike P. concreta, which has a pantropical distribution, P. bicolor/rosea is restricted to Madagascar, Comoros and the Seychelles (Mytnik-Ejsmont, 2011).
Brazilian sequences form a clade
The monophyletic nature of the group formed by the Brazilian accessions, contrasting with the paraphyletic group made up of African accessions, is consistent with the hypothesis that Africa is the center of origin with a subsequent (i.e., more recent) dispersal into the Neotropics (Russell et al., 2010a, 2010b).
Hybrid origins of some taxa suggested
The hybrid origin of P. concreta is a possible explanation for this species being found in different positions in the tree (Russell et al., 2010b). P. concreta individuals that have dispersed out of Africa are tetraploid, whereas plants found in continental Africa can be diploid or tetraploid. The sister taxa of African P. concreta are diploid (Russell et al., 2010b), indicating that tetraploidy is a derived state in P. concreta. Allotetraploidy in P. concreta has been confirmed by analysis of low copy nuclear genes (Russell et al., 2010a).
Interspecific hybridization events, as in P. concreta, are considered a source of chloroplast genome exchange via introgression. Chloroplast genome exchange among species is sometimes suggested as an explanation for the inconsistencies between phylogenetic trees based on nuclear and plastid markers in, e.g., Populus (Salicaceae) (Smith & Sytsma, 1990; Tsitrone, Kirkpatrick & Levin, 2003), Nothofagus (Nothofagaceae) and Crassulaceae (Mort et al., 2002; Acosta & Premoli, 2010). In Nothofagus, chloroplast capture results in the association of chloroplast genomes with geographic locations, rather than taxonomic relationships (Acosta & Premoli, 2010). Relationships based on geographic location could be explored as a possible explanation for the proximity of P. concreta (accesses from Brazil) with P. estrellensis (also from Brazil) and not with non-Brazilian accessions of P. concreta. In this case a study of nuclear markers of these taxa would be needed.
Relationships in the group that includes P. concreta, P. foliosa, P. estrellensis and other species are not well resolved due to the low sequence divergence levels between species found in both plastid and nuclear genes (Russell et al., 2010a, 2010b, 2011; Mytnik-Ejsmont, 2011). Generally, the morphological variation observed in this group is identified as P. concreta. Although P. estrellensis is considered a valid species on the official plant list of Brazil (Barros et al., 2010), there is no consensus on synonymy with P. concreta. This can be seen in the herbarium identifications that sometimes consider them as two distinct species, but sometimes as the same species. The same occurs with P. foliosa, a name which would only be correctly applied to plants from the Amazon basin, the Guyana Shield and the West Indies (Peraza-Flores, Fernández-Concha & Romero-González, 2011). This circumscription is not accepted by Mytnik-Ejsmont (2011), who considers P. estrellensis and P. foliosa to be synonymous.
Genetic dissimilarity between African and Neotropical tetraploids was reported by Russell et al. (2010a) and Russell et al. (2011), but the delimitation P. estrellensis, P. concreta and P. foliosa remained uncertain. According to our results, under a molecular perspective, P. estrellensis should be considered distinct from P. concreta. Moreover, our results do not corroborate the placement in synonymy of P. estrellensis and P. foliosa as proposed by Mytnik-Ejsmont (2011). In our tree P. foliosa forms a highly supported group with some P. concreta sequences (from samples collected in Brazil). Finally, although our results indicate a possible separation of Brazilian and African P. concreta, the delimitation of this species remains uncertain, considering that there is no generic taxonomic revision that has rigorously analyzed the morphological variation in this species. Moreover, considering the reticulated evolution by Russell et al. (2010a), further investigation with nuclear markers would be necessary.
Taken together, our analysis suggests that P. estrellensis can be considered a distinct species from P. concreta and P. foliosa, and that Brazilian and African P. concreta should probably be treated as different species. Evidence of hybridization influencing the evolution of P. concreta (Russell et al., 2010a, 2010b) highlights how importance it will be to also consider bi-parentally inherited nuclear DNA when inferring of phylogenetic relationships between this species and other species of the genus. The placement in synonymy of P. estrellensis and P. foliosa proposed by Mytnik-Ejsmont (2011) was not confirmed by this study. In our results, P. foliosa forms a highly supported clade including Brazilian samples of P. concreta.
Implications for data requirements
The entire chloroplast is more useful than a fast subset
By using a relatively large number of chloroplast sequences we were able to resolve the polytomy involving the Neotropical species. But, if on one hand this dataset is promising in the formulation of more robust phylogenetic hypotheses, on the other hand, the complete chloroplast genome sequencing may be costly for systematic projects that consider genera with many species (Särkinen & George, 2013), such as Polystachya, which has about 250 species. This was the main motivation for testing how well a smaller sample of the chloroplast genome would perform in phylogenetic inference. This was done in two ways: by duplicating a fast region and analyzing multiple copies of this dataset, and by sampling without replacement from all non-coding regions in our alignment.
We found that sampling without replacement up to ∼9 kb of non-coding sequence (16% of our alignment) was not sufficient to return a robust inference across all nodes. This was in contrast to the analysis of the entire chloroplast and showed that in the case of these samples of Polystachya, more data were needed to resolve their relationships. The cost of primers, amplification and Sanger sequencing of only three or four regions begins to exceed that of gene capture of the entire chloroplast. It is therefore more cost effective and produces a more robust result to undertake the collection of the entire chloroplast genome. That said, our 16% sample did resolve some nodes with high support, and other nodes obtained moderate to high support from a few of the replicates. This suggests that these data are on the way to resolving most nodes, but a gradual increase in resolving power occurs as characters are added.
Duplicating a single fast region even 16 times, in this case psbD-trnT copies totaling ∼14 kb, failed to achieve a robustly resolved phylogeny. The results for the psbD-trnT duplicated analysis was poorer even than that of sampling fewer but more representative characters across the non-coding region (above). It appears that a small sample size (only 870 bp of independent sequence sites) is a serious source of stochastic error in this case. Sampling one versus 16 copies of the same dataset only slightly increased the number of resolved nodes (but still falling short of the number of nodes usually resolved with support by the smaller 16% sample), confirming the limitations of the original dataset. It is likely that the original dataset simply did not contain sites that changed on most branches of the phylogeny during the span of history that we investigated.
Numerous previous studies have also examined which regions of the plastid genome are typically evolving faster than others (Small et al., 1998; Shaw et al., 2005, 2007, 2014). Prior to NGS methods, the aim was to identify the “best” regions, when sequencing only a limited number could be afforded in most projects. However, given current technology, we should shift our focus to whether a few of the “best” regions are cost effective compared to using the entire genome, as the latter has become affordable for even small phylogenetic projects.
Homoplastic hairpin inversions affect phylogenetic analysis
One issue raised here that is rarely taken into account in analyses of whole chloroplasts is that sequence patterns at the small scale, namely hairpin inversions of loops, can still have an effect on phylogenetic inference, despite using very large data sets. Our results indicate that at least some of the differences between the trees inferred using entire loops versus down-weighted loops were driven by these hairpin loop inversions. This kind of phylogenetic effect has been observed in other cases, although with smaller data sets (Kim & Lee, 2005; Joly et al., 2010). If loops invert in a single molecular event (as is currently believed: Kelchner & Wendel, 1996; Kim & Lee, 2005), such as an intra-molecular recombination, then there is no good reason to use each character state difference found between sequences in the entire loop in an analysis. This simply inflates the phylogenetic impact of a single event, treating it instead as many independent events (corresponding to the number of character state differences in the inversion), as also noted by Kim & Lee (2005). As shown here, a larger data set simply does not give license to ignore known analytical pitfalls.
Together, these findings show that sampling the entire chloroplast, analyzed carefully, is a better option than sampling a few (even a dozen or more) fast regions. This is true, at least in Polystachya, but a similar result has also found by other studies, such as Parks, Cronn & Liston (2009) for Pinus. Based on cost alone, it seems there is no benefit to be gained by screening the chloroplast for faster markers when there are many short branches in the particular tree, as there are here. Whole chloroplast analyses are likely to be a better way forward than sampling individual chloroplast markers in addressing many phylogenetic questions. If gene capture is used, as it was here, it is also very easy to add probes to unlinked nuclear regions, further increasing the power of this approach as a general solution to the issue of data sampling.
Our results show that significantly increasing the number of nucleotides can be an effective option in the phylogenetic inference of taxonomic challenging taxa, such as the orchid genus Polystachya. We generated complete chloroplast sequences of 40 Polystachya specimens using a combination of Illumina NGS sequencing and a sequence capture, which solved a notorious polytomy for Neotropical species. Our tests on how well a smaller sample of the chloroplast genome would perform in phylogenetic inference shows that the whole chloroplast is a better option than selecting just a few highly variable markers. Full plastid genomes appear particularly powerful when the expected tree is likely to contain many short branches, but nonetheless need to be analyzed with care.
DNA sequence alignment with loops.
The file contains the chloroplast sequence alignment and a Mr Bayes analysis block that excludes poorly aligned sites, but includes the entire inverted loop regions.
DNA sequence alignment with loops down-weighted.
The file contains the chloroplast sequence alignment and a Mr Bayes analysis block that excludes poorly aligned sites, as well as down-weighting the inverted loop regions by excluding all but one character in each inverted loop.