First Guatemalan record of natural hybridisation between Neotropical species of the Lady’s Slipper orchid (Orchidaceae, Cypripedioideae)

Ecological niche modeling analysis

Ecological niche modeling (ENM) analysis was used to define areas of potential hybridization between C. dickinsonianum, C. irapeanum and C. molle Lindl. The modeling was based on the maximum entropy method implemented in Maxent version 3.3.2 (Phillips, Dudík & Schapire, 2004; Phillips, Anderson & Schapire, 2006; Elith et al., 2011) based on species presence-only observations. The area of the analysis extended from −119.940 to −81.155 longitude and from 29.786 to 8.871 latitude. The list of localities was prepared based on information provided in herbarium specimens deposited in the following herbaria: AMES, AMO, BIGU, CAS, DS, herb. Hinton, K, LL, MEXU, MO, MSC, and WTU. In total, 27 georeferenced records for C. irapeanum, nine for C. molle and only three for C. dickinsonianum were gathered. Two datasets of localities were created (Fig. 1). The first included all the assembled data. To reduce sample bias, we applied spatial filtering in the second data set (Boria et al., 2014) and randomly removed localities that were within 25 km of one another, while retaining the most localities possible. With this approach the second dataset included 14 records for C. irapeanum, four for C. molle and three for C. dickinsonianum. Both datasets are provided as Supplemental Information.

Two groups of bioclimatic variables in 2.5 arc minutes (±21.62 km² at the equator) developed by Hijmans et al. (2005) were used, together with the altitudinal data (Table 2). The first group included all 19 variables. From the second dataset, we removed seven “bioclims”, due to their significant and mutual correlation (above 0.9) as evaluated by the Pearson correlation coefficient calculation computed using ENMTools v1.3. The following variables were excluded from the dataset: bio6, bio7, bio9, bio10, bio11, bio16 and bio17.

Initially, four models were created for each studied species. The first one was created based on all known localities of the studied species and all climatic variables. In the second, a reduced dataset of variables was used. The third included records with a reduced sampling bias and all variables. The last model was created based on the same locality dataset and a reduced variable dataset. In all these analyses the maximum number of iterations was set to 10,000 and the convergence threshold to 0.00001. The “random seed” option, which provided a random test partition and a background subset for each run, was applied. The run was performed as a bootstrap with 1,000 replicates, and the output was set to logistic. All operations on GIS data were carried out using ArcGis 9.3 (ESRI; https://www.esri.com/en-us/home) and QGIS applications.

Furthermore, to reduce overfitting (Radosavljevic & Anderson, 2014) of the models resulting from the small sample size, two additional analyses were made. In both these experiments the reduced locality and variables datasets were used with the same setting as described above. In the first study the regularization multiplier was set at 2 and in the second study it was set at 4.

The evaluation of the models was performed using the most common metric - area under the curve (AUC), which was automatically calculated by the MaxEnt application. The niche overlap between the three studied species was calculated using ENMTools v1.3.

Macromorphological features

Observations in situ have been conducted since 2008. The herbarium material was prepared according to standard classical taxonomy procedure and studied using a stereomicroscope. The comparative research was conducted at the following herbaria: AMO, BIGU, MA, P, W and UGDA. The following vegetative characters of individual plants were analyzed: stem (height, surface), leaves (number, size, shape), inflorescence (length, number of flowers), floral bracts/pedicellate ovary ratio, perianth segments (size and surface of tepals and lip), as well as gynostemium (size and shape of the staminodial shield).

DNA extraction, amplification, sequencing and sequences analysis

Total genomic DNA was extracted from 20 mg of silica-dried petals (Chase & Hills, 1991) from C. irapeanum (two specimens), the putative hybrid (two specimens), C. molle (one specimen), and C. dickinsonianum (two specimens) using a DNA Mini Plant (A&A Biotechnology, Gdynia, Poland), following the manufacturer’s protocol. The voucher for all specimens is Fred Muller s.n., Guatemala, BIGU. The nuclear ribosomal region spanning the internal transcribed spacers (ITS1 and ITS2) and the 5.8S rRNA gene (ITS), nuclear low copy gene Xdh and plastid gene matK were used for detection of the hybrid origin of specimens from Guatemala. ITS was amplified using the primers 17SE and 26SE (Sun et al., 1994). Xdh was amplified using the primers Xp551F and Xp1590R (Górniak et al., 2014). The gene matK was amplified with the following two primers: - 19F (Molvray, Kores & Chase, 2000) and 1326R (Cuénoud et al., 2002). Polymerase chain reactions (PCR) were carried out in a total volume of 25 µl, containing 5 µl 5× buffer, 1 µl 50 mM MgCl₂ (only plastid markers), 1 µl 5mM dNTPs, 0.5 µl of 10 µM of each primer, 1 µl dimethyl sulfoxide (DMSO) (only ITS and Xdh) and 1.0 unit of Blue Perpetual DNA polymerase (Eurx, Gdansk, Poland). Amplification conditions for ITS and matK were 94 °C for 4 min; 30× (94 °C, 45 s; 52 °C, 45 s; 72 °C, 1 min/2 min, respectively); and 72 °C, 7 min. A touchdown protocol was used for PCR amplification of the Xdh: the initial denaturation step (94 °C for 2 min) was followed by six cycles of 94 °C for 45 s, 55 °C (reducing 1 °C per cycle) for 45 sand 72 °C for 90 s. The next 28 cycles used 94°C for 45 s, an annealing step at 49 °C for 45 s, and 72 °C for 90 s. The final extension step used 72 °C for 5 min. PCR products were purified using a High Pure PCR Product Purification Kit (Roche Diagnostic GmbH, Mannheim, Germany). Cycle sequencing was performed using a Big Dye Terminator v 3.1 Cycle Sequencing Kit (Applied Biosystems Inc., ABI, Warrington, Cheshire, UK) with the same primers as were used for PCR amplification: 2.0 µl of 5× sequencing buffer, 1.0 µl of Big Dye Terminator with 1.5 µl of 1 µM primer, 1–4 µl of amplified product (30–90 ng/µl), and 0.5 µl DMSO and H₂O in a total of 10 µl reaction volume. Cycle sequencing conditions were as follows: 25 cycles each with 15 s denaturation (94 °C), 5 s annealing (52 °C) and 4 min elongation (60 °C). The sequences were generated on an ABI 3720 automated capillary DNA sequencer from Genomed LLC (Warszawa, Poland). Both strands were sequenced to assure accuracy in base calling. Finch TV (Geospiza) was used to edit the sequences, and the two complementary strands were assembled using AutoAssembler (ABI). Representatives of the sections of the genus Cypripedium (gene matK) were downloaded from GenBank: JQ182208 Cypripedium molle, JQ182205 Cypripedium debile, JQ182207 Cypripedium irapeanum, AF263649 Cypripedium calceolus, AY557208 Cypripedium calceolus, JQ182204 Cypripedium acaule, JQ182203 Cypripedium palangshanense, JQ182202 Cypripedium margaritaceum, JQ182206 Cypripedium subtropicum, JQ182201 Cypripedium californicum, JQ182200 Cypripedium passerinum, JQ182199 Cypripedium candidum, JQ182198 Cypripedium farreri, JQ182197 Cypripedium tibeticum, JN181460 Cypripedium fasciculatum, JN181459 Cypripedium bardolphianum, JN181458 Cypripedium japonicum, JN181457 Cypripedium flavum and EF079360 Selenipedium aequinoctiale. All sequences were aligned by eye using SeaView v. 4 (Gouy, Guindon & Gascuel, 2010). For detection of seed parent plastid data (matK) the matrix was analyzed using the PAUP* heuristic search method (Phylogenetic Analysis Using Parsimony *and Other Methods) version 4.0b10 (Swofford, 2002). The optimality criterion was the likelihood of tree-bisection-reconnection (TBR) branch swapping and the MULTREES option was in effect. The internal support of clades was evaluated by the bootstrap (Felsenstein, 1985) method with 500 replicates. The General Time Reversible model of substitution with gamma distribution (GTR+G) was selected as the best fitting model by Akaike information criterion in ModelTest v. 3.7 (Posada & Crandall, 1998). To show hybridization visual pairwise comparisons were made (ITS, XDH).

Embryological study

Four capsules from dry material were tested to assess the developmental stages of the ovules/seeds. The procedure of staining in tetrazolium chloride was used (TTC; Van Waes & Debergh, 1986, modified; M Rykaczewski, pers. comm., 2017). After pretreatments (10% glucose, 24 h; then 1% of sodium hypochlorite solution, pH 7.5, 30 min;) the pieces of placenta with ovules/seeds were incubated in 1% TTC in phosphate buffer (pH 7.5) at 40 °C for 24 h. The analyses of pieces were performed firstly under a stereomicroscope (Nikon SMZ 1500) and then examined under a Nikon Eclipse E 800 microscope equipped with differential interference contrast (DIC) optics.

The developmental stages were assessed for approx. 500 ovules/ seeds of each capsule (100 randomly selected ovules, 5 repeats).

Journal nomenclatural statement

The electronic version of this article in Portable Document Format (PDF) will represent a published work according to the International Code of Nomenclature for algae, fungi, and plants (ICN), and hence the new names contained in the electronic version are effectively published under that Code from the electronic edition alone. In addition, new names contained in this work which have been issued with identifiers by IPNI (International Plant Names Index) will eventually be made available to the Global Names Index. The IPNI LSIDs can be resolved and the associated information viewed through any standard web browser by appending the Life Science Identifier (LSID) contained in this publication to the prefix “http://ipni.org/”. The online version of this work is archived and available from the following digital repositories: PeerJ, PubMed Central, and CLOCKSS.

Ecological niche modeling analysis

The calculated AUC values for all the created models received high scores of over 0.9 (Tables 3 and 4). Based on this test, the most reliable models were created using all available occurrence and climatic data with default regularization multipliers (1). These are presented in Fig. 2. All other models are provided in Figs. 3–5. According to the most reliable models, the factors limiting the distribution of the three studied species are related to the altitude and temperature (temperature seasonality and mean temperature of the warmest quarter). However, in the models created with a reduced climatic variable dataset some additional factors were indicated as influencing the analysis, e.g., bio2, bio12, bio13, bio1, bio8 and bio19 (Table 5). In addition, their contribution in particular models varied between the species. The niche overlap statistics (Table 6) indicated that the highest probability of co-occurrence between the studied Cypripedium species is observed within C. dickinsonianum and C. irapeanum (I = 0.721, D = 435) and this was also confirmed in the same statistics calculated for three other datasets (Table 6).

The ENM analysis indicated several regions characterized by bioclimatic conditions suitable for the studied species located outside their known geographical ranges (Cribb & Soto-Arenas, 1993). For C. dickinsonianum, such areas may be found in the Mexican Volcanic Axis, the Southern Sierra Madre and the Chorotega volcanic front (Fig. 2A). The model of suitable niche distribution created for C. irapeanum is quite consistent with its known geographical range (Fig. 2B) with additional potential habitats in the eastern Sierra Madre del Sur. In the Chorotega volcanic front, the area indicated in ENM analysis as suitable for C. molle (Fig. 2C), no populations of this species have been found thus far.

The ENM analysis indicated two areas characterized by habitats suitable for all three studied species: the Sierra Madre de Chiapas and the Cordillera Neovolcánica. Within these regions the potentially available habitats for C. dickinsonianum, C. irapeanum and C. molle are separated by less suitable zones. The potential hybrid zones of C. irapeanum and C. molle are located in the eastern Sierra Madre del Sur.

Molecular analysis

Results from phylogenetic analyses based on the plastid matK gene are presented in a phylogram (Fig. 6). Bootstrap support (BS) above 50% is given for supported clades above branches. The matK tree can be divided into two highly supported clades (A = 99 and B = 100). Clade A consists of species represented by various sections of Cypripedium. The base of the tree (clade B) comprises three species from section Irapeana. C. irapeanum together with the putative hybrid composing one clade which is a sister to C. molle. C. dickinsonianum is a sister to them. Pairwise alignment of nuclear ITS, Xdh and the plastid sequence comprising the 5′end of the intron trnK and matK gene revealed significant differences between Cypripedium irapeanum and C. dickinsonianum. Within the ITS sequence, four substitutions were observed—two transversions and two transitions. In addition to the sequence of ITS2, an indel of 15 base pairs in length occurred. Sequences (chromatograms) of putative hybrids are noisy (weak) from that site (sequences from two different alleles overlap each other making chromatograms unreadable). This feature was observed in both forward and reverse strands (see the chromatogram file provided as Supplemental Information). Within the Xdh sequence, seven substitutions were observed, of which five were transitions. Two specimens of Cypripedium which exhibited characteristics of hybrids in polymorphic sites have double peaks corresponding to nucleotides found in both species (Table 7). Comparison of the plastid sequence between C. irapeanum and C. dickinsonianum showed an indel of seven base pairs in length at the 3′trnK intron and five transitions and one indel in the matK gene. C. irapeanum and the putative hybrid have identical sequences of the matK gene. Comparison of molecular markers identified three substitutions between C. irapeanum and C. molle, one in each of the analyzed markers. DNA data matrices are provided as Supplemental Information.

Taxonomic treatment

Due to the detection of gene flow between C. dickinsonianum and C. irapeanum and mixed morphological characters of the population discovered by Mr Muller in Guatemala we decided to describe it as the first, natural hybrid in the section Irapeana under the name Cypripedium × fred-mulleri.

Cypripedium × fred-mulleri Szlach., Kolan. & Górniak, hybr. nov.

Diagnosis: Cypripedium × fred-mulleri is characterized by having flowers 5.2–7 cm across, elliptic, acute dorsal sepal, oblong-elliptic, obtuse petals, deeply saccate, obovoid-globose lip and trullate, acute staminode. It differs from C. irapeanum in its smaller flowers, deeper color (closer to C. dickinsonianum), density of windows on the lip, and form of dorsal sepal and petal apex. From C. dickinsonianum it is distinguished, inter alia, by the shape of the staminode and lip as well as by the petal form.

Type: Guatemala, Alta Verapaz. South of Cobán. 30 May 2013. F. Muller s.n. (BIGU! 309 holotype). UGDA-DLSz! - drawing of type, photos.

Description: Plants up to 75 cm tall, densely and softly hairy throughout. Stem erect, rather stout. Leaves up to 15 cm, distributed along the stem, 3–8 cm long, 2.8–3.8 cm wide, ovate to ovate-lanceolate, acute. Inflorescence 15–33 cm long, loosely 5–8-flowered. Flowers showy, large, yellow. Floral bracts 4–6.6 cm long, ovate-lanceolate, acute. Pedicel up to 1 cm long, pubescent. Ovary up to 2.5 cm long, pubescent. Dorsal sepal 3–3.8 cm long, 1.7–2 cm wide, elliptic, acute, margins pilose. Petals 3.4–4.3 cm long, 1.6–2.1 cm wide, oblong-elliptic, obtuse, pilose, especially near the base. Synsepal 2.2–3.2 cm long, 1.6–2 cm wide, elliptic, obtuse, bifid or occasionally free to the base, margins pilose. Lip 3.5–4 cm long, 2.8–3 cm wide, deeply saccate, obovoid-globose, margins incurved around the lip opening, with translucent windows all over the surface. Staminode 1.2 cm long, 0.7–0.8 cm wide, trullate, acute. Capsule 2.2–2.6 cm long. Figure 7.

Paratypes: Guatemala. Alta Verapaz, South of Cobán. 25 Jun 2009. (Muller - photo!); The same location 26 Jun 2010. (Muller - photo!).

Etymology: Dedicated to the discoverer of this hybrid, Fred Muller.

Distribution: Known so far to be exclusively from the Guatemalan department of Alta Verapaz. Due to the vulnerability of populations of C. irapeanum, C. dickinsonianum and C. × fred-mulleri to illicit harvesting, the exact locality is not given. The known localities of C. irapeanum are distributed from Central Mexico to Guatemala and Honduras while the currently known range of C. dickinsonianum is discontinuous, extending from eastern Chiapas (México), through the Sierra de los Cuchumatanes and the Sierra de Chamá to the central Honduran uplands (although herbarium vouchers are currently lacking Dix & Dix, 2000). Figure 8.

Ecology: The hybrid population was found on a south-oriented limestone hillside at an altitude of about 1,500 m. The plants grow in an open, seasonally dry pine-oak forest with Brahea dulcis (Kunth) Mart. (Arecaceae) and species of Agave L. (Asparagaceae). Other terrestrial orchid species occurring in this area are: Cyrtopodium punctatum (L.) Lindl., Stenorrhynchos pubens (A. Rich. & Galeotti) Schltr. and Dichromanthus cinnabarinus (La Llave & Lex.) Garay. Moreover, two species of Bletia Ruiz & Pav. have been reported from this location. The hybrid plants begin blooming in mid-May, at the beginning of the rainy season. The flowers have been observed as late as at the end of July, which is the beginning of the flowering season for both C. irapeanum and C. dickinsonianum in nearby colonies. Field observations in 2013 suggested that the population might have benefited from a recent wild fire, as a significant increase in the number of flowering specimens had previously been recorded in the season following a fire at the locality.

Notes: Morphologically, C. × fred-mulleri is transitional between its parental species in many respects, as we describe in Table 1. The hybrid occupies an intermediate position in the general size of the plant, number of leaves and inflorescence length. For example, according to literature data and our own study, the inflorescence of C. irapeanum reaches up to 40 cm in length, whereas in C. dickinsonianum it is less than 9 cm. The length of inflorescence of C. × fred-mulleri is between 15 and 33 cm. This inflorescence can bear five to eight flowers. The reported number of flowers per inflorescence in C. irapeanum is up to 12, and in C. dickinsonianum it is between two and six (Figs. 9 and 10). Even the number, size and distribution pattern of diaphragmatic windows on the lip is manifestly transitional between both parental species. In C. dickinsonianum the diaphragma is outspread between somewhat thickened, dendritic veins and cover ca 40% of the total lip surface. On the other hand, in C. irapeanum the windows are relatively small and occupy less than 10% of the lip surface. The lip of C. × fred-mulleri, although in form similar to the ovule parent, is covered by diaphragma in a similar pattern as in its pollen parent which cover ca 30% of whole lip surface. In some respects, however, C. × fred-mulleri is more similar to its pollen parent (densely hairy stem and leaves, length of the leaf blade), but in some others to its ovule parent. This set of characteristics concern the form and width of the leaf blade, general flower architecture, and the length and general form of staminode. It is noteworthy that in C. × fred-mulleri all taxonomically important characters useful in determination of Neotropical Cypripedium species, i.e., number of flowers per inflorescence, size of the flower segments and generative parts, are intermediate between parental species.

Key to the taxa of Cypripedium sect. Irapeana

Ovule and seed development

The seeds inside four open-pollinated flowers did not react with TTC as was indicated under stereomicroscope (Figs. 11A–11E). Deeper analysis revealed that the enlarged capsules contained a mix of unfertilized (Figs. 12A–12G) or embryo-bearing ovules (9.2–26.2%; Fig. S1). In the two ovaries, a small number of embryo-bearing ovules was accompanied by many unfertilized ovules that were at maturity (Figs. 12F–12G) or aborted (Fig. 11E). Very early stages (from the zygote to a few-celled proembryo) of embryo development were detected inside ovules/young seeds (Figs. 11E, 12H–12J). The ovaries without embryos contained ovules that were mostly at the bisporic stage (Figs. 12B–12C), sporadically at modified monosporic megasporogenesis stages (Fig. 12D), or at megagametogenesis (Fig. 12E).