Introduction

There are some rare ecosystems on earth that demand extreme adaptations from their inhabitants. One of these places is the Atacama Desert. It is known to be the world's oldest and driest non-polar desert, which has been predominantly hyperarid (at least in the core zone) since the Miocene or even longer (Hartley et al. 2005; Ritter et al. 2018). Despite this severe aridity, the Atacama Desert is habitat for specialized plants building local vegetation that can only exist because of unique climatic conditions occurring in this area (Schulz et al. 2012; Bull et al. 2018). The Atacama Desert is stretching along the southern Peruvian and northern Chilean coast with its core zone of hyperaridity between 15 and 30°S (Houston and Hartley 2003). The eastward extent of the desert is limited by the Andean Cordillera. It only stretches from 20 to around 100 km inland (Rundel et al. 2007) with altitudes between sea level and 3500 m above sea level (m a. s. l.) (Houston and Hartley 2003). In these arid and hyperarid regions, fog is often the most important source of water for plant life. Fog zones occur mainly in the cliffs close to the coastline above 400 m and in mountain slopes facing the predominant winds (S and SW) that advect the stratocumulus clouds. Concentrated against the hillsides they allow the development of plant communities, which are called "lomas" (Spanish for small hills) formations or fog oases. The mountainous coast, especially in northern Chile, tends to trap the coastal fog and further enhance the aridity of the continental interior. Only fog corridors connecting the coast with the intermediate depression enable loma formation also in the inland (Alpers and Brimhall 1988; Rundel et al. 1991; Cereceda et al. 2004). The existence of vegetation depends heavily on the altitude and where the top and bases of the stratocumulus clouds occur in the region (Cereceda et al. 2004).

Frequent members of loma formations are from the family of Bromeliaceae, most commonly known from tropic or subtropic climates, and among Bromeliaceae members of the genus Tillandsia often contribute to loma vegetation. Tillandsia L. is the largest genus within the family of Bromeliaceae and respective subfamily Tillandsioideae. The subfamily is composed by four tribes: Catopsideae, Glomeropitcairnieae, Vrieseeae and Tillandsieae including 22 genera in total and comprising c. 1400 species (Barfuss et al. 2005, 2016; Gomes-da-Silva and Costa 2013; Leme et al. 2017). The more than 600 epiphytic and sometimes terrestrial species of the genus Tillandsia are distributed from southern United States to central Argentina and Chile (Smith and Downs 1977). The genus was traditionally divided into seven subgenera: Allardtia A. Dietr., Anoplophytum (Beer) Baker, Diaphoranthema (Beer) Baker, Phytarrhiza (Vis.) Baker, Pseudalcantarea (Beer) Baker, Pseudo-Catopsis Baker and Tillandsia (Smith and Downs 1977). These subgenera turned out to be mostly non-monophyletic; therefore, a new subdivision was suggested: Aerobia Mez in C.DC., Anoplophytum, Diaphoranthema, Phytarrhiza, Pseudovriesea Barfuss & W. Till, Viridantha (Espejo) W. Till & Barfuss, Tillandsia and several unclassified species complexes (Barfuss et al. 2016). Subgenus Allardtia is treated now in synonymy of subgenus Tillandsia, and Pseudalcantarea was elevated to generic rank (Barfuss et al. 2016). Several Tillandsia species invaded sandy soil habitats with adaptations like epiphytic or unrooted (epiarenitic) growth, utilization of crassulaceaen acid metabolism (CAM) and specialized leaves with water absorbing trichomes (Rundel et al. 1997; Zizka et al. 2009). In fact, they possess rudimental roots but without function regarding water uptake or solidification on the surface (Rauh 1985). Terrestrial communities with species from the genus Tillandsia, also called "tillandsiales," are quite abundant in the Atacama Desert but their existence follows a distinct pattern and is associated to the occurrence of above-mentioned fog corridors (Pinto et al. 2006). In Chile tillandsiales are distributed in altitudes from 900 to 1200 m in a linear area of 470 km from Arica (18°20′S) to the Loa River (21°25′S) with a large gap from Camarones (19°15′S) to Iquique (20°12′S). Furthermore, they can be found in ranges from 3 to 45 km inland. Lomas show different level of degradation according to their geographic position. They prosper well in belts near the coast growing on slopes toward SW and W and become sparse toward inland localities. They also show evidence of dieback in their lower limits showing the struggle and dependence of atmospheric moisture (Cereceda et al. 2004; Pinto et al. 2006; Schulz et al. 2011). Pinto et al. (2006) recorded three different terrestrial Tillandsia species in northern Chile: Tillandsia landbeckii Phil., Tillandsia marconae W.Till & Vitek and Tillandsia virescens Ruiz Pav. Tillandsia marconae only occurs close to Arica and T. virescens is found in small patches only with close proximity to T. landbeckii lomas. Most fog oases observed by Pinto et al. (2006) were monospecific stands of T. landbeckii (Koch et al. 2019, 2020), which is the main object of this study. Tillandsia landbeckii belongs to the subgenus Diaphoranthema, which may include a total of 29 species (Donadío et al 2015; but see also Barfuss et al. 2016). They are characterized by few (1–10) reduced flowers and grow on trees, bushes, rocks or bare sand in altitudes from sea level to 4300 m (Donadío et al. 2015). Some species within the subgenus Diaphoranthema are facultatively autogamous or even cleistogamous, which makes them independent from pollinators (Till 1992). This seems to be advantageous, considering the extreme habitat with possibly rare occurrence of potential pollinators as shown for T. landbeckii (Koch et al. 2019).

Drastic environmental changes from wet and cold to dry and warm periods within Pleistocene glacial–interglacial cycles might have been a strong impulse for population divergence, species distribution and speciation by influencing gene flow, demographic expansion or contraction and genetic bottlenecks. Unfortunately, there are only very few studies about how these fluctuations in climatic conditions during the Quaternary impacted plant species and especially loma formations in the Atacama Desert. Gilmartin (1983) for example hypothesized that there was a rapid speciation of xeromorphic groups of Tillandsia during the Pleistocene. However, since hyperarid conditions exist for much longer at this area, it seems to be possible that some of the adapted species are much older, as well. According to Granados Mendoza et al. (2017) migration of members of the subgenus Diaphoranthema, which includes T. landbeckii, might have started from the Andes and Southern Chile toward the Brazilian shield. Furthermore, the origin of the lomas in Peru show subtropical connections, while the flora of the Chilean lomas show more affinities to the flora of Central Chile and the adjacent Andes. In Peru, the most widely distributed loma forming Tillandsia species is T. purpurea Ruiz Pavon, and in the respective contact zone with T. landbeckii, T. marconae has been found (Zizka and Munoz-Schick 1993). There is convincing evidence that T. marconae is a hybrid between both species, T. landbeckii and T. purpurea (Bratzel et al. 2020). This coincides with the assumption of a strong phytogeographic boundary existing between Peru and Chile (Rundel et al. 1991, 2007) and independent evolutionary histories of the various species can be hypothesized since T. purpurea has been placed in subgenus Phytarrhiza, which also appeared to be a polyphyletic subgenus (Barfuss et al. 2016). Among tribe Tillandsieae terrestrial growth forms have evolved on species-group/clade-level at least five times independently (Barfuss et al. 2016).

Herein we present a time calibrated subfamily-wide plastome phylogeny demonstrating stem and crown group age of Tillandsia landbeckii. Since T. marconae has been identified as a hybrid between T. purpurea and T. landbeckii, we also tested the hypothesis that as a consequence T. marconae plastome transmission and/or plastome capture must be younger than the crown group age of T. landbeckii. The onset of T. marconae evolutionary history might denote a major environmental change shifting distribution areas of Peruvian T. purpurea lomas southward and Chilean T. landbeckii lomas northward and consequently resulting in a sympatric narrow distribution range.

Material and methods

Study design

There is limited information on divergence time estimates and calibration points in Bromeliaceae in general and subfamily Tillandsioideae in particular. The most comprehensive study has been presented by Givnish et al. (2011). From this study, we selected a set of taxa representing crucial nodes being used as calibration points in our analysis. Previous taxonomic and phylogenetic results and conclusions are almost exclusively based on plastid DNA information; therefore, we expanded earlier few gene analyses (Givnish et al. 2011; Barfuss et al. 2016) to a whole plastome analysis. The respective largely expanded amount of DNA sequence information can then be used to resolve also divergence time estimate for shorter periods in time of few thousands or hundred-thousands of years (e.g., Hohmann and Koch 2017; Hohmann et al. 2018) and, thereby, resolving also node ages of herein investigated T. landbeckii.

Plant material

Taxon sampling from Bromeliaceae is representing subfamily Tillandsioideae and representatives from its sister clade (selected taxa from subfamilies Pitcairniodeae, Bromeliodieae, Puyoideae and Hechtioideae) (Online Resource 1). We also selected a representative Tillandsia landbeckii accession from Chilean terrestrial loma vegetation and Peruvian accessions from higher altitudes (T. landbeckii ssp. andina – MHJB-B1687/HEID113097, 2300 m a. s. l.; T. landbeckii var. rigidior – MHJB-B1684/HEID113094, 3100 m a. s. l.) allowing to include the species´ geographic, ecological and climatic disjunction (Peru versus Chile). Furthermore, two accessions from T. purpurea and T. marconae have been investigated, because T. marconae was proposed hybrid between T. purpurea and T. landbeckii with a very restricted intervening distribution range. Two of these accessions (T. purpurea – P06, and T. marconae – P07) are growing sympatrically in the same loma vegetation near Tocyasca, Peru. The other two accessions are also from Peru (T. purpurea: HEID104854/Rauh38703; T. marconae: HEID103591/Rauh20899). Plants are cultivated at Heidelberg Botanical Garden as permanent living collections (HEID). This selection of accessions allows testing some key hypotheses on the timeline of evolutionary history of T. landbeckii in Peru and Chile.

For genome size estimates we used 15 accessions grown at Heidelberg Botanical Garden representing 3 species (T. landbeckii, T. marconae, T. purpurea) (Online Resource 1).

Two additional plastomes were selected from NCBI representing T. usneoides (KY293680) and Ananas comosus (AP014632 [= NC026220]) with detailed information on its respective structure and function (Poczai and Hyvönen 2017; Nashima et al. 2015). Additional plastomes from Ochagavia carnea (NC045386), Puya mirabilis (NC045380), and Fasciluaria bicolor (MN563795) (Paule et al. 2020a, b) have been used for further comparisons of gene content and selected sequence characteristics.

DNA extraction, plastome sequencing, assembly and high-quality reference plastomes

DNA was extracted from fresh leaf material of all samples using the Invisorb Spin Plant Mini Kit (STRATEC Biomedical, Birkenfeld, Germany) following the manufacturer’s protocol or using a modified CTAB protocol optimized for Bromeliaceae (Bratzel et al. 2020). DNA integrity was checked on 1% agarose gels and DNA concentration was measured using the Invitrogen Qubit 2.0 fluorometer kit (Thermo Fischer, Freiburg, Germany). Preparation and sequencing of total genomic DNA libraries were performed at the CellNetworks Deep Sequencing Core Facility (Heidelberg, Germany). Library insert size was between 200 to 400 bp and the SMARTer® ThruPLEX® Tag-seq kit (Takara Bio Inc, Saint-Germain-en-Laye, France) was used for library preparation. Sequencing was performed on an Illumina NextSeq 550 sequencing system in 150-bp paired-end mode. Raw Illumina reads were filtered for quality by retaining only reads with a Phred score > 30 and longer than 50 bp and adapter sequences were removed both by using the program trimmomatic (Bolger et al. 2014). For further analyses, only reads were used where both paired reads passed the filtering criteria. For generating reference-based assemblies, the trimmed reads were mapped to an earlier published T. usneoides plastid genome (KY293680) using BWA (mem algorithm; Li 2013) setting the penalty for unpaired reads to 15. Prior to mapping one of the copies of the inverted repeat region was removed from the reference sequence as duplicate sequences in the reference produce secondary alignments. To increase mapping quality further, the mappings were filtered for reads with a mapping quality > 1 and duplicate reads were removed both using samtools (Li et al. 2009). Gatk3 tools (McKenna et al. 2010) were used for enhancing alignment quality in case of presence of indels and for variant calling (SNP as well as indel polymorphisms). The gatk3 tool FastaAlternateReferenceMaker was used for generating sequences including the detected SNPs and indels and regions with a coverage < 5 and a mapping quality < 30 were masked using samtools maskfasta (Li et al. 2009). All plastid sequences obtained from the reference-based mapping approach were aligned to the reference sequence and 120 gene annotations were transferred using the program Geneious (version 7.1.7, Biomatters Ltd., Auckland, New Zealand). Coding regions as well as regions encoding tRNAs or rDNAs were extracted using bedtools (Quinlan and Hall 2010) except for ycf1 with sequence overlapping with inverted repeat and aligned using mafft (Katoh et al. 2002). All resulting 112 gene alignments were then concatenated using the script catfasta2phyml (https://github.com/nylander/catfasta2phyml). The alignment generated from 112 coding sequences of the plastid genome (including rDNA and tRNA encoding sequences) had a total length of 71,429 bp and was used for phylogenetic reconstruction and divergence time estimates (see below; DRYAD submission https://doi.org/10.5061/dryad.kh189324d).

In addition to the reference-based approach Unicycler v0.4.8 (Wick et al. 2017) was used for generating de novo assemblies of four accessions (T. usneoides (HEID109235), T. purpurea (HEID104854), T. espinosae (HEID130781), and T. malzinei (HEID108286)). The de novo assemblies were indexed as blast libraries and searched for contigs matching the plastid genome by using T. usneoides complete plastome (NCBI KY293680) as query sequence. By aligning the matching contigs to the query sequence order and orientation of contigs were determined, and primers were designed for amplifying regions between contigs. PCR products were Sanger sequenced and used for filling gaps between contig ends. Detailed primer information and PCR protocols are provided with Online Resource 2.

Phylogenetic analysis and divergence time estimates

In order to determine the most suitable molecular evolutionary model as well as partitioning schemes for phylogenetic reconstruction PartitionFinder 2.1.1 (Lanfear et al. 2017) was used. All models implemented in RAxML-ng (Kozlov et al. 2019) as well as all models implemented in BEAST2 (Bouckaert et al. 2014) were tested using the corrected Akaike Information Criterion (AICc) for model selection in a greedy search, and branch lengths were allowed to be linked. RAxML-ng (Kozlov et al. 2019) was used for phylogenetic reconstruction taking the alignment described before as well as the list of partitions and according most suitable models as input. Calculations started from 20 parsimony trees generated by the program, branch length was set to linked and 1000 bootstrap replicated were run. The resulting tree with added bootstrap support values was visualized using the program FigTree v1.4.2 (Rambaut and Drummond 2012).

Divergence times estimates were performed using the program BEAST2 (Bouckaert et al. 2014) using the alignment as well as the partition scheme and according models described above and in the previous section and trees and clock models were linked. In respect to the best-fitting models K81 had to be excluded as BEAST2 cannot handle those models. In the six cases where K81 had been determined by PartitionFinder 2.1.1 as best-fitting model the second best model was chosen instead (subset 4: TIM + I, subset 7: TIM + I + G, subset 12: TIM, subset 15: TVM + I, subset 18: TIM and subset 33: K80). A relaxed log normal molecular clock was applied allowing for rate heterogeneity across branches, the tree resulting from the RAxML-ng analysis was used as starting tree and the Birth Death model was set as prior for all partitions. Secondary divergence time estimates were constrained using estimates from a plastid 8-locus dataset (Givnish et al. 2011): (A) crown group age of subfamiliy Tillandsioideae (T. usneoides versus Catopsis floribunda; split time 14 mya with a SD of ± 0.5 mya); (B) crown group age of Puyoideae being a member of the sister group of subfamily Tillandsioideae (Puya laxa versus P. alpestris; split time 9.8 mya with a SD of ± 0.7 mya). The Markov Chain Monte Carlo (MCMC) algorithm was run for 100 million generations and every 1000th tree was sampled, and six independent runs using identical settings were performed. Tracer v1.7.1 (Rambaut et al. 2018) was used for determining whether the stationary phase was approached, whether the effective sampling size was > 200 and whether all runs had converged. Using LogCombiner v2.5.1 trees from all 6 runs were combined and 20% of the data were excluded as burnin. TreeAnnotator v2.5.1 was used for adding age estimates as Median heights to the nodes before visualizing the tree by FigTree v1.4.2. The tree is provided with DRYAD (https://doi.org/10.5061/dryad.kh189324d).

DNA barcoding of Tillandsia purpurea and T. marconae

The herein selected accessions of Tillandsia, and in particular T. purpurea and T. marconae, have been also used to analyze their phylogenetic position in a genus-wide phylogeny of Tillandsia and subfamily Tillandsioideae presented by Barfuss et al. (2016). In this earlier study, a preliminary analysis of plastid trnK-matK-trnK (Large Single-Copy region) (alignment length of 1887 bp.), rpoB-trnC-petN (LSC-Region) (alignment length of 3225 bp.) and ycf1 (Small Single-Copy Region) (alignment length of 4722 bp.) provides a comprehensive phylogenetic overview, and it was indicated that T. purpurea represents a species complex. However, they did not include T. marconae, which our research group has shown earlier to be of a putative hybrid origin (Bratzel et al. 2020). We extracted the respective marker regions from our data and re-run the large-scale phylogenetic tree (Barfuss et al. 2016) including our herein analyzed Tillandsia accessions using RaxML-ng with the following settings: raxml-ng —all —msa —brlen linked —model GTR —tree pars{20} —bs-trees 200 (Kozlov et al. 2019).

Genome size estimates

Genome size estimates were performed for selected individuals to analyze potential ploidal level variation among the three species T. purpurea, T. landbeckii and T. marconae to collect further information on hybrid speciation processes (e.g., homoploid hybrid speciation). Nuclear DNA content was determined using flow cytometry following a simplified protocol (Doležel et al. 2007). Approximately, 10 mm2 of fresh and young leaf tissue from each plant to be analyzed was chopped together with an appropriate volume of respective internal reference standards (RaphanuItas sativus cv. Saxa, 1.11 pg/2C; Solanum lycopersicum cv. Stupicke, 1.96 pg/2C; Zea mays cv. CE-777, 5.43 pg/2C) (Doležel et al. 1992) using a sharp razor blade in a Petri dish containing 0.5 mL of ice-cold Lysis buffer LB01 (Doležel et al. 1989) supplemented with 2 µg mL − 1 β-mercaptoethanol. The suspension was filtered through a 30-μm CellTrics® filter (Sysmec-Partec GmbH, Münster/Görlitz, Germany) and 1.5 mL of LB01 buffer was added as well as propidium iodide and RNase (both to a final concentration of 50 μg mL − 1). After 90 min incubation on ice, the relative fluorescence intensity of 5000 particles was recorded using a flow cytometer (CyFlow Space; Sysmex-Partec GmbH, Münster/Görlitz, Germany) equipped with a green (532-nm) solid state laser. We applied the following stringent criteria to get precise and stable flow cytometric results: (i) only analyses where the coefficient of variation of the sample peak was below 5% were taken into account, (ii) each sample was measured twice on different days to minimize potential random instrumental drift (Doležel and Bartoš 2005), and (iii) if the between-day variation exceeded a 5% threshold another measurement was made and the most remote value was discarded when the sample was re-analyzed. The histograms were evaluated with the FloMax FCS 2.0 program (Sysmex-Partec GmbH, Münster/Görlitz, Germany).

Results and discussion

Bromeliaceae high-quality reference genomes

Bromeliad phylogenetics and evolutionary research is largely relying on plastid markers, therefore there is some need of respective reference plastomes. Complete chloroplast genomes from the family of Bromeliaceae are published for six species only (Ananas comosus—AP014632, KR336549, KU598872; Tillandsia usneoides—KY293680; Ochagavia carnea—MN563799; Ochagavia elegans – MN563796 and MN563798; Puya mirabilis—MN563797; Fascicularia bicolor ssp. bicolor—MN563795), and only one represents the genus Tillandsia. Only the plastomes of Ananas (Nashima et al. 2015; Redwan et al. 2015) and Tillandsia usneoides (Poczai and Hyvönen 2017) have been discussed in detail, the other GenBank submissions remained largely unexplored in earlier studies (Liu et al. 2016; Paule et al. 2020a, b).

Within this study four additional plastomes representing four Tillandsia species (T. usneoides, T. malzinei, T. purpurea, T. espinosae) were de novo assembled from genome skimming data nearly doubling the number of available plastome data. Herein, newly presented plastomes had an average size of 157,163 ± 847 bp (Online Resource 3), and their structure was the same quadripartite structure as in already published plastid genomes [Large Single-Copy region (LSC) neighboured by Inverted repeat b (IRb), followed by Small Single Copy region (SSC) and Inverted Repeat a (IRa); length of segments is shown in Table 1) (Gitzendanner et al. 2018). The two inverted repeat regions were identical so they assembled into one contig, which was inserted as IRa into the scaffolded plastomes, and connections between contigs as well as short sequences between contigs were confirmed and added by Sanger sequencing. This led to slight size differences between IRa and IRb (3 to 266 bp). In all four plastomes, 94 genes could be annotated in LSC and SSC (72 protein coding, 22 tRNAs, 0 rRNAs) and 20 genes in the IRa/IRb regions (8 protein coding (including ycf15), 8 tRNAs, 4 rRNAs) making up a total of 134 genes (gene list is found with DRYAD submission https://doi.org/10.5061/dryad.kh189324d). However, two of these genes were found to be partially deleted in two accessions: In T. usneoides 106 bp were missing from rpl23 (usually 174 bp full length) and in T. purpurea 23 bp were missing from psaI (usually 111 bp full length). In addition, ycf1 is located in all four generated plastomes on the border of SSC and IRa and creates a pseudogene of roughly 1000 bp (1089 bp in T. purpurea and T. espinosae, 1077 bp in T. malzinei and 1091 bp in T. usneoides). The partial duplication is due to the IR expansion (Poczai and Hyvönen 2017). Hence, difference in plastome size of T. usneoides of KY293680 compared to our herein presented data of 2.4% are largely because of rpl32 and ycf1 (pseudogene) sequence information not present in our accession. Furthermore, more than 1000 bps (position 6000–7100 spacing rps16 and trnQ) have not been found in our accession of T. usneoides, but BLAST searches demonstrate that this sequence information is distributed with very high sequence identity in the Ananas comosus nuclear genome (chromosomes 1, 2 and 7). However, this may be validated in future studies, e.g., via long-read sequencing approaches or fragment-specific and individual PCR-based tests relying on unique primer combinations.

Table 1 Total length and region lengths (in basepairs) of herein presented reference plastomes including large single copy (LSC), small single copy (SSC) and inverted repeat (IR) regions a and b
Full size table

Both, the already published plastomes and the ones generated in this study contain overall the same set of genes, with the exception of ycf15, which has been only described for the published T. usneoides plastome (Poczai and Hyvönen 2017) as well as the four Tillandsia plastomes generated in this study. Despite this, we were able to find the sequence of 256 bp in Ananas comosus, Puya mirabilis, Ochagavia carnea and Fascicularia bicolor ssp. bicolor, as well. This was verified with a BLAST search against Arabidopsis thaliana. This gene is considered to be a pseudogene (Poczai and Hyvönen 2017).

In comparison to already published plastid genomes we identified 134 genes in our newly assembled plastids. This means there is a minor discrepancy in gene content of four genes to Ochagavia carnea and to Puya mirabilis and 1 gene to Fascicularia bicolor subsp. bicolor (Table 2). Order of genes and annotations among studied plastomes are identical (gene list provided with DRYAD submission https://doi.org/10.5061/dryad.kh189324d). Accession codes for raw sequence reads are provided with Online Resource 1 and can be found with NCBI under BioProject ID PRJNA701548.

Table 2 Comparison of major features among plastomes from bromeliad species deposited with GenBank
Full size table

Phylogenetic reconstructions are congruent with traditional concepts

The results from maximum likelihood analysis of the plastome data (Fig. 1) are congruent with recent phylogenetic analyses of the Bromeliaceae in general and subfamily Tillandsiodeae in particular (e.g., Givnish et al. 2011; Barfuss et al. 2016).

Fig. 1
figure 1

Maximum likelihood-based phylogeny of plastome coding regions. Bootstrap support is indicated for values > 95%. Tillandsia species in focus of this study are marked with gray boxes (T. landbeckii, T. purpurea, T. marconae). The tree is rooting according to Givnish et al. (2011) and Barfuss et al. (2016), and distance scale represents mean expected rate of substitutions per site

Full size image

Bootstrap support in our analysis is most often very high. In case of Puya laxa (Puyoideae) we obtained a similar result as Givnish et al. (2011) demonstrating in ML analysis an unresolved subfamily Puyoideae with deeply branching Puya taxa (bootstrap support for paraphyly 51%). BEAST analysis (Fig. 2) recovered the same phylogenetic topology and, thereby, confirming a robust phylogenetic signal of the dataset. Tillandsia marconae, which is of hybrid origin (Bratzel et al. 2020), appears at two different positions in the phylogenetic tree indicating a putative multiple and polytopic origin. The T. purpurea accession grouped together with one T. marconae and all T. landbeckii accessions might have captured its plastome via reticulate evolutionary processes. The second T. marconae accession (HEID103591/Rauh20899) has been analyzed earlier in a broader context of the evolution of CAM metabolism in Bromeliaceae (Hermida-Carrera et al. 2020). For this taxon and respective accession CAM metabolisms has been confirmed and maternal phylogenetic affinity with T. purpurea has been documented sequencing the rbcL gene (Hermida-Carrera et al. 2020). Our second herein analyzed T. purpurea accession (HEID104854) unexpectedly did not group with this T. marconae (Fig. 1). However, in an earlier barcoding study using the nuclear gene agt1 this accession clearly shows a signature of the T. purpurea agt1 gene (Bratzel et al. 2020). This leads to the conclusion that this T. purpurea accession might be also of hybrid origin or have been affected by introgression and/or plastid capture of an unknown species. And indeed, Werner Rauh, who collected the plants in 1975 (http://scriptorium.cos.uni-heidelberg.de/entry/?20233; Koch et al. 2013) added a notice, that it is an unusual “microform” (see next paragraph).

Fig. 2
figure 2

BEAST chronogram with divergence time estimates indicated (blue bars represent 95% confidence interval). Tillandsia species in focus of this study are marked with gray boxes (T. landbeckii, T. purpurea, T. marconae). Divergence time estimates relevant to our understanding of T. landbeckii evolutionary history are indicated in red

Full size image

DNA barcoding indicates a putative reticulate origin of Tillandsia marconae

Maximum-likelihood analysis of herein investigated Tillandsia accessions integrated into the large-scale phylogeny provided by Barfuss et al. (2016) should have guided us to an explanation for the phylogenetic positions of T. marconae and T. purpurea. The ML phylogeny is shown in Online Resource 4. This analysis provided some valuable information on T. purpurea (HEID104854). This accession has been also re-collected later in 1986 by W. Till (Vienna, coll. WT2144) and was named as dwarf form of T. kunthiana Gaudich. (T. latifolia group) with similar plants resembling true T. purpurea found nearby (W. Till, coll. WT2146). Therefore, we are confident that the accession analyzed herein is a result of hybridization processes between the two species T. purpurea and T. latifolia, both of them are well-known to build up Peruvian tillandsia lomas (Aponte and Flores 2013; Hesse 2012, 2014; Rundel and Dillon 1998). Sympatric occurrence of the various terrestrial Tillandsia has not been documented in a systematic way yet. Fortunately, there is few documentation as illustrated with Fig. 3. In Peru there is T. latifolia growing sympatrically with T. purpurea and T. paleacea (Fig. 3b). In Chile near Iquique (Cerro La Isla) there are few T. virescens individuals at the top of the hill growing sympatrically with T. landbeckii (Fig. 3a). The study site at Tocyasca (Peru) analyzed herein has sympatrically T. purpurea and T. marconae growing and near Ancon (Peru) there is sympatrically growing T. landbeckii and T. purpurea (Fig. 3c). Another recent survey of various locations near Tacna (Peru) exemplified sympatric occurrence of T. purpurea, T. landbeckii and T. werdermannii Harms in various combinations (Lazo Ramos 2011). In summary this may provide evidence that all these seven terrestrial and on sand growing (epiarenic) Tillandsia species from Chile and Peru may have genetic contact among each other.

Fig. 3
figure 3

a Sympatrically growing Tillandsia virescens and T. landbeckii (background) close to Iquique (Cerro La Isla, Chile) (image taken by M. Koch, Nov 2016). b Sympatrically growing T. latifolia, T. purpurea and T. paleacea at Reserva Nacional de Ancon (Peru) (image taken by M. Douglas, 9 Dec 2017). c Sympatrically growing T. landbeckii and T. purpurea 15 km NE of Ancon (Peru) (image taken by M. Douglas, 9 Dec 2017). d Monospecific Loma vegetation of T. landbeckii at Cerro Guanacos (Iquique, Chile) (image taken by M. Koch, 22 Nov 2016)

Full size image

Genome size estimates indicate homoploid and reticulate evolution among Tillandsia landbeckii and T. purpurea

Estimated genome sizes for the three species T. landbeckii, T. purpurea and T. marconae are provided with Suppl. Mat. Table 1. For T. landbeckii genome sizes vary from 2.23–2.61 pg (n = 10, mean 2.41 pg, SD 0.11), T. purpurea varies from 2.11–2.96 pg (n = 3, mean 2.40 pg, SD 0.39), and T. marconae showed a mean of 2.28 pg (n = 2, SD 0.01). We found no difference among measurements from runs with different standards, and here we present the values referring to the Raphanus standard. Individual chromosome counts of these accessions are not available. However, for T. purpurea a haploid chromosome number (n = x = 25) has been reported (Brown and Gilmartin 1989), confirming the predominant number of 2n = 2x = 50 found within Bromeliaceae (e.g., Cotias de Oliveira et al. 2000; Paule et al. 2020a, b). These results allow us to assume that all three species are on the same ploidal level. Considering that two analyzed T. purpurea and T. marconae accessions from the same loma vegetation carry a T. landbeckii plastome type, it can be hypothesized that not only T. marconae evolved with varying parental contribution (maternal contribution from T. landbeckii and T. purpurea), but also T. purpurea is subjected to reticulate evolutionary processes and thereby capturing T. landbeckii plastomes, which is best explained by additional (back)crossing events toward T. purpurea.

Since herein investigated T. purpurea from Tocyasca (Peru) shows genetic contact with T. latifolia, and since also genetic contact among the other epiarenic species is likely, we hypothesize a common diploid base chromosome number of 2n = 50 among all these species not hindering respective gene flow (Paule et al. 2020a, b).

BEAST analysis indicates Pleistocene crown ages of Tillandsia landbeckii evolutionary history

Divergence time estimates (crown and stem group ages) among the five selected subfamilies do not differ between our study and Givnish et al. (2011) (Fig. 2, Table 3).

Table 3 Crown and stem group ages of subfamilies included in the chronogram
Full size table

This has to be expected, because we used divergence times from this study to calibrate our phylogeny. However, it demonstrates the robustness of the datasets, and it further allows to consider our herein presented divergence time estimates for T. landbeckii with confidence. Considering T. mollis and T. usneoides as closest relatives of T. landbeckii (Barfuss et al. 2016) stem and crown group ages of T. landbeckii are 2.038 and 0.570 mya, respectively. This takes into consideration that T. purpurea and T. marconae from the same clade (here growing in the same Tillandsia field) carries respective T. landbeckii plastome types transmissed maternally from T. landbeckii into these species. Accordingly, the origin of hybrid T. marconae indicating past genetic contact between Chilean T. landbeckii and Peruvian T. purpurea is estimated with a minimum of 0.112 mya (Fig. 3). This fits well with a rainy northern Atacama Desert in Peru during the last Interglacial (appr. 125 kya) (Contreras et al. 2010), which may have enforced Peruvian T. purpurea to migrate southward and hybridized with T. landbeckii. The onset of genetic differentiation in T. landbeckii (crown group age) 0.570 mya is not coinciding with any obvious documented environmental/bioclimatic impact in the Chilean Atacama Desert and, therefore, a trigger of diversification remains unknown.

It is remarkable that within subgenus Diaphoranthema (mostly epiphytic species) T. landbeckii belongs to a clade of few Tillandsia species comprising mostly epiarenic or epilithic species, with epiarenic species being quite rare in the genus (Suppl. Material Fig. 1). This supports the view that epiarenic life form is ancestral to this clade, which originated in the early Pliocene appr. 4.67 mya. In the Atacama Desert a semiarid climate persisted from 8 to 3 mya, punctuated by a phase of increased aridity at ca. 6 mya (Hartley and Chong 2002). Therefore, onset of the evolution of this particular clade might have followed this punctuated phase. The stem group age of T. landbeckii coincides with the onset of the ENSO climate system (Hartley and Chong 2002), which also denotes the radiation of the most diverse genus of angiosperm plants in the Atacama Desert, Nolana (Guerrero et al. 2013). The El Niño-Southern Oscillation (ENSO) is a recurring climate pattern on periods ranging from 3–7 years involving changes in the temperature of waters in the central and eastern tropical Pacific Ocean. The resulting oscillations of warming and cooling ocean water (ENSO cycle) directly affects rainfall and fog distribution in the tropics in general and in the Atacama Desert in particular.

Conclusions and outlook

Tillandsia landbeckii is a fascinating monospecific loma forming species in the Atacama Desert (Fig. 3d) starting its evolutionary history nearly 2 mya. Despite extreme environmental conditions in a hyperarid region and a presumably very narrow ecological niche we found evidence for ample gene flow among the various epiarenic and loma forming species including T. landbeckii and also ongoing hybrid speciation processes. There is some evidence that these processes might be triggered by major environmental change and significant distribution range shifts (T. marconae), but we are actually lacking detailed knowledge to understand the spatio-temporal dynamics behind. We aim in our future studies to compare phylogeographic and temporal evolutionary patterns among species to evaluate in some more detail potential drivers of speciation processes in epiarenic and loma forming Tillandsia species growing at the dry limits since latest from the Early Miocene onward.