Integrating restriction site-associated DNA sequencing (RAD-seq) with morphological cladistic analysis clarifies evolutionary relationships among major species groups of bee orchids

Abstract

Background and aims: Bee orchids (Ophrys) have become the most popular model system for studying reproduction via insect-mediated pseudo-copulation and for exploring the consequent, putatively adaptive, evolutionary radiations. However, despite intensive past research, both the phylogenetic structure and species diversity within the genus remain highly contentious. Here, we integrate next-generation sequencing and morphological cladistic techniques to clarify the phylogeny of the genus.

Methods: At least two accessions of each of the ten species groups previously circumscribed from large-scale cloned nuclear ribosomal internal transcibed spacer (nrITS) sequencing were subjected to restriction site-associated sequencing (RAD-seq). The resulting matrix of 4159 single nucleotide polymorphisms (SNPs) for 34 accessions was used to construct an unrooted network and a rooted maximum likelihood phylogeny. A parallel morphological cladistic matrix of 43 characters generated both polymorphic and non-polymorphic sets of parsimony trees before being mapped across the RAD-seq topology.

Key results: RAD-seq data strongly support the monophyly of nine out of ten groups previously circumscribed using nrITS and resolve three major clades; in contrast, supposed microspecies are barely distinguishable. Strong incongruence separated the RAD-seq trees from both the morphological trees and traditional classifications; mapping of the morphological characters across the RAD-seq topology rendered them far more homoplastic.

Conclusions: The comparatively high level of morphological homoplasy reflects extensive convergence, whereas the derived placement of the fusca group is attributed to paedomorphic simplification. The phenotype of the most recent common ancestor of the extant lineages is inferred, but it post-dates the majority of the character-state changes that typify the genus. RAD-seq may represent the high-water mark of the contribution of molecular phylogenetics to understanding evolution within Ophrys; further progress will require large-scale population-level studies that integrate phenotypic and genotypic data in a cogent conceptual framework.

Keywords: Biogeography; Mediterranean; Ophrys; RAD-seq; character mapping; coalescence; convergence; evolution; internal transcribed spacer; macrospecies; microspecies; morphology; paedomorphosis; phylogenetics; plastid; pseudo-copulation; systematics.

Fig. 1.

Flowers of 13 microspecies representing the nine molecularly circumscribed macrospecies (groups) of Ophrys discussed in the present phylogenetic study, together with four further subgroups created for specific use in our morphological cladistic analysis. (A) O. episcopalis, Crete (fuciflora group, fuciflora subgroup: H’1), (B) O. insectifera, UK (insectifera group: A), (C) O. regis-ferdinandii, Chios (speculum group: D), (D) O. grandiflora, Sicily (tenthredinifera group: B), (E) O. bombyliflora, Sardinia (bombyliflora group: C), (F) O. apifera, Sicily (apifera group: F), (G) O. oestrifera, Chios (fuciflora group, scolopax subgroup: H’2), (H) O. bornmuelleri, Cyprus (umbilicata group, bornmuelleri subgroup: J2), (I) O. kotschyi, Cyprus (umbilicata group, umbilicata subgroup: J1), (J) O. israelitica, Cyprus (fusca group), (K) O. spruneri, Crete (sphegodes group, sphegodes subgroup: G2), (L) O. argolica, Peloponnese (sphegodes group, argolica subgroup: G1), (M) O. bertolonii, Sicily (sphegodes group, bertolonii subgroup: G3). Labels on (A): la, labellum (lip); lp, lateral petal; ms, median sepal; ls, lateral sepal; g, gynostemium (column); sc, stigmatic cavity; tc, temporal callosity (pseudoeye); bf, basal field; sp, speculum; a, appendix.

Fig. 2.

Unrooted SplitsTree network based on 4060 RAD-seq-derived SNPs for 32 plants that together represent the ten putative Ophrys macrospecies (A–J) illustrated in Fig. 1. Inset: magnified view of topology for representatives of groups G–I. Details of samples are given in Table 1.

Fig. 3.

Rooted RAxML tree of RAD-seq data for the same 32 plants that formed the basis of Fig. 2, plus two outgroup accessions. The tree is based on 4159 high-quality, filtered SNPs. Values above the branch are bootstrap values, and groups A–J of Devey et al. (2008) are labelled. Details of samples are given in Table 1.

Fig. 4.

Comparison of topologies obtained in previous phylogenetic studies of Ophrys, reduced to the ten macrospecies (labelled A–J) recognized by Devey et al. (2008). (A) Devey et al., 2008, fig. 2; ITS, MP. (B) Devey et al., 2008, fig. 2; three plastid regions, MP. (C) Soliva and Widmer, 2001, fig. 1; ITS + one plastid region, MP. (D) Breitkopf et al., 2015, fig. 1; six low-copy nuclear genes, ML. (E) Present study; RAD-seq, ML. MP, maximum parsimony; ML, maximum likelihood. Numbers associated with branches are bootstrap values. Dashed branches represent only a single analysed sample and so do not test monophyly of the relevant macrospecies. It was necessary to shift the horizontal position of the fusca group, E, in (B), and to interpolate a putative additional group based on O. heldreichii, K, in (D). In (B), the basal position of the insectifera group, A, was dictated by its use as the de facto outgroup. Question marks denote indistinguishable groups.

Fig. 5.

Morphological cladograms generated via maximum parsimony from a matrix of 13 ingroup plus three outgroup species. (A) One of the nine most-parsimonious trees generated from the initial matrix that included polymorphic cells. (B) One of three most-parsimonious trees generated from the present morphological cladistic matrix after all polymorphic cells had been resolved in favour of the most frequent character state within each. Arrowed nodes collapsed in the respective strict consensus trees. Branch lengths reflect Acctran optimization. Numbers on branches are bootstrap support values.

Fig. 6.

(A) Integrated phylogeny generated by constraining the 43 characters (33 informative) of the non-polymorphic morphological matrix to the topology dictated by the RAD-seq tree illustrated in Fig. 3. Character numbers and states reflect those given in Appendix 1. Acctran optimization; non-homoplastic character-state transitions shown in bold, homoplastic character-state transitions shown in italics. (B) Plot of standardized branch lengths for molecularly constrained morphology against the RAxML RAD-seq tree; both sets of branch lengths have been standardized to unit variance.