# The Explosive Radiation of the Neotropical *Tillandsia* Subgenus *Tillandsia* (Bromeliaceae) Has Been Accompanied by Pervasive Hybridization

GIL YARDENI<sup>1,2,*</sup>, MICHAEL H. J. BARFUSS<sup>1,†</sup>, WALTER TILL<sup>1,†</sup>, MATTHEW R. THORNTON<sup>1,‡</sup>, CLARA GROOT CREGO<sup>1,3,§</sup>, CHRISTIAN LEXER<sup>1,†</sup>, THIBAULT LEROY<sup>1,4,§</sup>, AND OVIDIU PAUN<sup>1,*,§</sup>

**Abstract.**—The recent rapid radiation of *Tillandsia* subgenus *Tillandsia* (Bromeliaceae) provides an attractive system to study the drivers and constraints of species diversification. This species-rich Neotropical monocot clade includes predominantly epiphytic species displaying vast phenotypic diversity. Recent in-depth phylogenomic work revealed that the subgenus originated within the last 7 myr, with one major expansion from South into Central America within the last 5 myr. However, disagreements between phylogenies and lack of resolution at shallow nodes suggest that hybridization may have occurred throughout the radiation, together with frequent incomplete lineage sorting and rapid gene family evolution. We used whole-genome resequencing data to explore the evolutionary history of representative ingroup species employing both tree-based and network approaches. Our results indicate that lineage co-occurrence does not predict relatedness and confirm significant deviations from a tree-like structure, coupled with pervasive gene-tree discordance. Focusing on hybridization, ABBA-BABA and related statistics were used to infer the rates and relative timing of introgression, whereas topology weighting uncovered high heterogeneity of the phylogenetic signal along the genome. High rates of hybridization within and among subclades suggest that, contrary to previous hypotheses, the expansion of subgenus *Tillandsia* into Central America proceeded through several dispersal events, punctuated by episodes of diversification and gene flow. Network analysis revealed reticulation as a plausible propeller during radiation and establishment across different ecological niches. This work contributes a plant example of prevalent hybridization during rapid species diversification, supporting the hypothesis that interspecific gene flow facilitates explosive diversification. [bromeliad; gene-tree discordance; hybridization; phylogenomics; Neotropical diversity; species network; *Tillandsia*.]

Rapid evolutionary radiations are characterized by accelerated and substantial diversification, usually following dispersal into a novel geographic area. As a myriad of diverse ecological spaces become available, a lineage can opportunistically adapt to quickly occupy them in a process termed adaptive radiation (Linder 2008; Rundell and Price 2009; Hughes et al. 2015; Stroud and Losos 2016; Naciri and Linder 2020). Rapid radiations provide ample opportunities to investigate the mechanisms driving diversification. Although not all rapid radiations are adaptive, that is, not necessarily associated with an increase in ecological occupancy (Schluter 2000; Givnish 2015; Stroud and Losos 2016), evolutionary radiations are generally connected with habitat heterogeneity (Paun et al. 2016; Soltis et al. 2019), climatic fluctuations (Slovak et al. 2023), landscape fragmentation (Hughes and Atchison 2015), and/or orogenesis (Boschman and Condamine 2022). Ubiquitous among angiosperms in general and in the Neotropics in partic-

ular, plant radiations remain to date understudied and elusive, especially when compared with animal systems (but see Richardson et al. 2001; Linder 2008; Drummond et al. 2012; Lagomarsino et al. 2016; Paun et al. 2016; Pérez-Escobar et al. 2017).

Phylogenomic studies of young, rapid radiations face multiple challenges. Low sequence divergence and incomplete reproductive barriers increase the likelihood of short internal branches and impede phylogenomic resolution (Whitfield and Lockhart 2007; Straub et al. 2014; Giarla and Esselstyn 2015). Furthermore, population-level processes during rapid radiations—such as fluctuations in population size and incomplete lineage sorting (ILS)—generate high rates of gene-tree conflict, potentially reducing the phylogenetic signal (Oliver 2013; Parins-Fukuchi et al. 2021). Such processes were traditionally deemed as analytical noise when aiming to reconstruct bifurcating species trees. Yet an increasing amount of research indicates that episodes of phyloge-

<footnote n="1"><sup>1</sup>Department of Botany and Biodiversity Research, University of Vienna, 1030 Vienna, Austria</footnote>

<footnote n="2"><sup>2</sup>Department of Biotechnology, Institute of Computational Biology, University of Life Sciences and Natural Resources, 1190 Vienna, Austria</footnote>

<footnote n="3"><sup>3</sup>Vienna Graduate School of Population Genetics, 1210 Vienna, Austria</footnote>

<footnote n="4"><sup>4</sup>GenPhySE, Université de Toulouse, INRAE, ENVT, 31326 Castanet Tolosan, France</footnote>

<footnote>*Correspondence to be sent to: Gil Yardeni, Department of Botany and Biodiversity Research, University of Vienna, 1030 Vienna, Austria; E-mail: gil.c.yardeni@gmail.com; Ovidiu Paun, Department of Botany and Biodiversity Research, University of Vienna, 1030 Vienna, Austria; E-mail: ovidiu.paun@univie.ac.at</footnote>

<footnote>†Deceased.</footnote>

<footnote>Thibault Leroy and Ovidiu Paun contributed equally to this article.</footnote>

<footnote>Received 21 November 2023; reviews returned 5 May 2025; accepted 25 June 2025</footnote>

<footnote>Associate Editor: Michael Matschiner</footnote>
<figure><img src="image_1.png" /><figcaption>FIGURE 1. A simplified phylogeny of the subclades sampled in this study, as generally hypothesized in previously published studies (Barfuss et al. 2016; Pinzón et al. 2016; Vera-Paz et al. 2023). The shared grouping of the South American (SA) subclades denotes disagreements in the literature. Subclade names correspond to naming in this study, following previously proposed names.</figcaption></figure>

nomic conflict represent typical signatures of microevolutionary processes and are often associated with rapid phenotypic shifts (Filiault et al. 2018; Parins-Fukuchi et al. 2021). Hence, nonbifurcating relationships should be regarded and studied as fundamental biological phenomena of diversification, underlying the formation of species (Malinsky et al. 2018; Slovak et al. 2023; Wogan et al. 2023).

Historically, hybridization and introgression were viewed as processes that impede speciation by homogenizing distinct genomes, introducing potentially maladaptive alleles, and disrupting locally adapted gene networks. However, they also recently emerged as key drivers of species diversification (e.g., Abbott et al. 2013; Seehausen 2014), with gene flow now regarded as an important process during speciation (Mallet et al. 2016; He et al. 2019; Linck and Battey 2019). Gene flow between recently diverged species can increase genetic diversity, offering substrate for natural selection. Interspecific hybridization can prompt ecological diversification through heterosis and adaptive introgression, a process in which advantageous alleles are transferred between gene pools, enhancing adaptation. Furthermore, introgression can introduce advantageous alleles upon which selection has already acted or create novel allele combinations, thus catalyzing evolutionary radiations (Abbott et al. 2013; Harrison and Larson 2014; Suarez-Gonzalez et al. 2018; Edelman et al. 2019; Meier et al. 2019).

The rapidly radiating and highly diverse Neotropical *Tillandsia* subgenus *Tillandsia* provides a particularly relevant study system to investigate the drivers of radiations. Subgenus *Tillandsia* is the most diverse of 7 subgenera within genus *Tillandsia* (family Bromeliaceae), comprising more than 250 predominantly epiphytic species. *Tillandsia* species are distributed across a wide range, from the southeastern United States to central Argentina (Barfuss et al. 2016). Members of the subgenus exhibit tremendous morphological diversity of adaptive traits, which allows them to occupy dis-

parate habitats, from tropical rainforests to deserts, and from lowlands to highlands (Benzing 2000; Givnish et al. 2011). Adaptations such as adventitious roots, trichomes modified for water and nutrient absorption, a water-impounding leaf rosette tank, and Crassulacean acid metabolism (CAM) appear to have undergone in some cases correlated and in others contingent evolution, resulting in adaptive syndromes (Givnish et al. 2014).

Early efforts employing conserved plastid and nuclear markers circumscribed the monophyly of the subgenus and its clades, yet recovered largely inconsistent relationships between subclades (Barfuss et al. 2005; Chew et al. 2010; Barfuss et al. 2016; Pinzón et al. 2016; Granados Mendoza et al. 2017). Using plastome phylogenomics, Vera-Paz et al. (2023) confirmed the monophyly of the subgenus and the presence of 3 main clades within the subgenus, including a monophyletic, radiated "clade K" with mostly epiphytic and ornithophilous members distinctly distributed in North and Central America (Fig. 1; Barfuss et al. 2016; Granados Mendoza et al. 2017). They further identified K.1 and K.2 as subclades of clade K. Within this study's sampling, the *T. punctulata* (K.1) and *T. utriculata* (K.2.1) subclades include widespread or endemic species with intermediate CAM-C3 syndromes or various CAM-like photosynthetic syndromes, respectively. The *T. guatemalensis* subclade (K.2.2) is a small group of widespread and predominantly C3 plants, exhibiting smooth leaves and a prominent tank habit. Finally, the *T. fasciculata* subclade (previously K.2.3) includes species largely exhibiting CAM photosynthetic syndrome with typical water-absorbing trichomes. In addition, species within this subclade present varying distributions, from endemic to widespread, some extending into the Antilles (Till 2000; Espejo-Serna et al. 2004).

The subgenus *Tillandsia* was estimated to be phylogenetically young: Existing Bromeliaceae phylogenies set the crown age of the subfamily *Tillandsioideae* at 15.2 ± 0.4 Ma and of "core tillandsioids" at 9.6 ± 0.7
Ma (Givnish et al. 2014). Several studies (Winkler 1990; Till 2000; Givnish et al. 2014; Vera-Paz et al. 2023) proposed that the ancestor of “clade K” was most likely South American, expanded to Central America approximately 4.9 Ma and subsequently dispersed multiple times to North America. However, the relationships within “clade K” remained elusive and disagreements between gene trees are common, calling to consider the role of different evolutionary processes in the subgenus’ radiation. As a result, hybridization was put forward as especially prevalent in the subgenus (Till 2000; Vera-Paz et al. 2023).

In contrast to phylogenetic reconstruction using a
limited number of markers, whole-genome sequencing
provides an abundance of informative data, which in-
creases resolution and helps disentangle major causes of
incongruence (e.g., Liu et al. 2015; Guo et al. 2023). In
comparison to plastid markers, nuclear loci are charac-
terized by higher evolutionary rates and biparental in-
heritance, rendering them suitable for resolving rapid
radiations and investigating hybridization (Léveillé-
Bourret et al. 2018; McKain et al. 2018; Debray et al.
2019). Due to costs and the complexities of molecular
lab work, a trade-off often exists between the scale of
sequencing and the extent of taxon sampling in phy-
logenomic studies (Mitchell et al. 2000). Regardless,
genome-wide molecular phylogenetics generally results
in greater phylogenetic informativeness and improved
accuracy (Rokas et al. 2003; Wortley et al. 2005; Léveillé-
Bourret et al. 2018; Debray et al. 2019).

Using whole-genome resequencing data and a re-
cently published reference genome, our study aims to
shed light on the evolutionary history of *Tillandsia* sub-
genus *Tillandsia*, focusing on Central American groups
within the radiated “clade K” and representatives of
South American species. Specifically, we seek to in-
fer relationships between representative species and in-
vestigate signals of phylogenetic conflict, with the ul-
timate aim to examine patterns of gene flow, which
may have contributed to adaptive trait shifts within this
Neotropical rapid radiation. Our detailed analyses re-
veal a reticulated history within the subgenus and
at least 2 colonization events from South into Cen-
tral America, surprisingly indicating that the previ-
ously described “clade K” (Barfuss et al. 2016; Granados
Mendoza et al. 2017; Vera-Paz et al. 2023) may not con-
stitute a monophyletic group.

MATERIALS AND METHODS

Plant Material

Leaf material was collected from 69 individuals be-
longing to 36 *Tillandsia* species. Material from field col-
lections was cut lengthwise and immediately dried in
powdered silica gel. Samples from the botanical collec-
tion of the University of Vienna were extracted fresh
(Supplementary Tables S1 and S3). Sixty-seven acces-

sions correspond to 34 species of *Tillandsia* subgenus *Tillandsia*, whereas the remaining 2 belong to subgenus *Allardtia* sensu Smith and Downs (1977) and were used as outgroups. Of the 34 ingroup species, 26 represent the Central-North American “clade K” radiation and 8 species represent South American *Tillandsia* (henceforth SA clades; Fig. 2; Supplementary Tables S1). Albeit not exhaustive, our sampling was performed with the intention to represent the variety of morphological and physiological syndromes within clades SA and “K” (Supplementary Tables S2). We sampled species from all subclades of clade K previously introduced: 6, 5, 2, and 13 from subclades *T. punctulata*, *T. utriculata*, *T. gutemalensis*, and *T. fasciculata*, respectively.

DNA Library Preparation

DNA extractions were performed as previously de-
scribed (Yardeni et al. 2022) using a modified CTAB pro-
tocol (Doyle and Doyle 1987) or the QIAGEN DNeasy®
Plant Mini Kit (Qiagen, USA). They were further pu-
rified using the NucleoSpin® gDNA clean-up kit from
Macherey-Nagel (Hudlow et al. 2011). The extracts were
subsequently diluted in water and quantified with a
Qubit® 3.0 Fluorometer (Life Technologies, Carlsbad,
CA, USA). Prior to library preparation, a maximum of
600 ng of DNA per accession was sheared at 4 °C to a
target average length of 400 bp using a Bioruptor® Pico
sonication device (Diagenode, Denville, NJ, USA). Illum-
mina libraries were prepared using a modified KAPA
protocol with the KAPA LTP Library Preparation Kit
(Roche, Basel, Switzerland) with adaptor ligation and
adaptor fill-in reactions based on Meyer and Kircher
(2010). Some libraries have been instead prepared with
the NEBNext Ultra II DNA PCR-free Library Prep Kit
(New England Biolabs, Ipswich, MA, USA). Samples
were double indexed either with a set of 60 dual-
index primers, as recommended by Kircher et al. (2012)
and described in Loiseau et al. (2019), or with Il-
lumina TrueSeq PCR-free dual indexes. The libraries
underwent size selection steps using AMPure Beads
(Agencourt). Finally, the libraries were pooled and se-
quenced at the Vienna BioCenter Core Facilities on Illum-
mina HiSeqV4 PE125 or on NovaSeqS1 PE150. Protocol
details for all accessions are listed in Supplementary
Tables S3.

Data Processing

The raw sequence data were demultiplexed using
deML v.1.1.3 (Renaud et al. 2015) and bamtools v.2.5.1
(Barnett et al. 2011), and converted from BAM to FASTQ
format using bedtools v.2.29.2 (Quinlan and Hall 2010).
Reads were trimmed for adapter content and quality
using trimgalore v.0.6.5 (Krueger 2019), a wrapper tool
around fastqc and cutadapt, using FastQC retain un-
paired. Sequence quality and adapter removal were con-
firmed with FastQC v.0.11.9 (Andrews 2010).
<figure><img src="image_2.png" /><figcaption>FIGURE 2. A coalescent-based species tree generated on 15,791 genomic windows with ASTRAL-III for 64 individuals representing 33 species of *Tillandsia* subgenus *Tillandsia*, plus 2 outgroups of subgenus *Allardtia* sensu Smith and Downs (1977). Branch lengths are given in coalescent units. Node values represent local posterior probabilities for the main topology and are equal to 1 unless noted otherwise. Pie charts at the nodes show levels of gene-tree discordance: the percentages of concordant gene trees (blue), the top alternative bipartition (green), other conflicting topologies (red), and uninformative gene trees (gray). Clade terminology follows previously proposed names (e.g., Barfuss et al. 2016; Pinzón et al. 2016).</figcaption></figure>

Quality- and adapter-trimmed reads were aligned to the *T. fasciculata* reference genome v.1.0 (GenBank Assembly accession GCA_029168755.1; Groot Crego et al. 2024) using bowtie2 v.2.3.5.1 (Langmead and Salzberg 2012) with the –very-sensitive-local option to maximize alignment rates. Low-quality mapped reads (mapQ < 10) were removed, alignments were sorted by position using samtools v.1.15.1 (Danecek et al. 2021),

and PCR duplicates were marked using MarkDuplicates from PicardTools v.2.25.2 (Picard Toolkit 2018). The average coverage was calculated with samtools v.1.15.1 (Li 2011).

To call variants, we used GATK HaplotypeCaller v.4.1.9.0 followed by joint calling with GenotypeGVCFs (Poplin et al. 2018). We proceeded to confirm technical replicates and to identify high relatedness between
accessions with kinship coefficients inferred through KING v.2.2.6 (Manichaikul et al. 2010). BAM files of technical replicates were merged before we called a final variant call file (VCF). Samples exhibiting a first-degree relationship with another sample (consistent with full-sib or direct parent-offspring relationships; i.e., kinship >0.177) or having an average genome-wide coverage below 3.5× were removed from subsequent analyses (Supplementary Tables S3). The resulting VCF was filtered using bcftools v.1.15 (Danecek et al. 2021) and GATK SelectVariants (DePristo et al. 2011) to exclude indels and any SNPs located within 3 bp of indels. Regions annotated as transposable elements were excluded using bedtools intersect (Quinlan and Hall 2010). We generally used the transition/transversion (Ts/Tv) ratio for guidance to define our filtering parameters, aiming for values leading to a set of SNPs exhibiting the highest ratio. Ts/Tv ratios were calculated with SnpSift (Cingolani et al. 2012). We finally used the following parameters for filtering: mapping quality <15, read depth coverage <4, quality by depth <4, Fisher strand bias >40, strand odds ratio >3, minor allele frequency (MAF) <0.045 (corresponding to the presence in at least 6 chromosomal complements), and missing rate <0.2. We prepared additional files with different MAF filtering thresholds, to examine the possible bias in branch lengths (see later and Supplementary Files 3). Summary statistics were generated using bcftools stats (Danecek et al. 2021). The full variant calling pipeline is available at https://github.com/giyany/Tillandsia_Phylo_hybridization/blob/main/notebooks/raw_to_vcf_GATK.md.

## Phylogenetic Tree Inference

We inferred phylogenetic relationships for all samples using both a maximum-likelihood (ML) and a coalescent-based method. We included a coalescent-based method in order to explicitly account for ILS, which may otherwise result in high support for an incorrect topology (Kubatko and Degnan 2007). Gene-tree incongruence further provides insight into molecular genome evolution, including the extent of ILS and genetic processes such as hybridization and introgression (Wendel and Doyle 1998; Galtier and Daubin 2008; Hibbins et al. 2020; Morales-Briones et al. 2021). Both analyses were repeated on data sets with different MAF filtering thresholds (Supplementary Files 3). The ML tree was inferred on a concatenated supermatrix, including both variant and invariant sites. First, a concatenated matrix was obtained by converting the VCF into a phylip file with vcf2phylip (Ortiz 2019). We then inferred a phylogeny with IQ-TREE v.2.1.3, using IQ-TREE's ModelFinder to select the best-fitting partitioning scheme and models for each genic region (Nguyen et al. 2015; Kalyaanamoorthy et al. 2017). Node support was inferred with 1000 nonparametric ultrafast bootstrap replicates (Chernomor et al. 2016). During the first

inference, the available *T. zoquensis* accession grouped within *T. fasciculata* (Supplementary Fig. S2), suggesting a need to revise its species status. We removed this sample from the data set and all analyses hereafter. We then repeated the analysis, constructing a whole-genome ML tree as described, as well as separate ML trees for each reference chromosome.

To reconstruct a species tree, ML trees were estimated on nonoverlapping genomic windows of 10 kb. Windows with fewer than 40 SNPs were excluded. Because loci with insufficient signal may reduce the accuracy of species tree estimation (Mirarab 2019), nodes with a bootstrap support below 10 were collapsed across all trees with newick utilities v.1.1.0 (Junier and Zdobnov 2010). ASTRAL-III v.5.7.8 (hereafter ASTRAL) was then used to infer a species tree, measuring branch support as posterior probabilities (Zhang et al. 2018). To explore gene-tree discordance, we used phyparts and calculated quartet support for the main, as well as the first and second alternative topologies. We calculated the number of concordant and discordant bipartitions on each node, imposing a cutoff of 70 for informative bootstrap values (-s 70; Smith et al. 2015). Gene-tree discordance was visualized with phypartspiecharts.py (available from https://github.com/mossmatters/MJPythonNotebooks).

## Investigating Monophyly and Deviations from a Bifurcating Tree Structure

We used several variations of the ABBA-BABA test (D-statistic) to further explore the relationships within and between subclades and to assess deviations from tree structure (Green et al. 2010; Durand et al. 2011; Martin et al. 2015). The D-statistic implemented in Dsuite v.0.5r45 (Malinsky et al. 2021) specifically uses allele frequency estimates, which allows to include several individuals per taxon and does not require the implicit assumption that the outgroup is fixed for the ancestral allele. For all D-statistic-related analyses, we used the full VCF with *T. complanata* as an outgroup. Significance was assessed using a jackknife procedure with a window size of 200 kb and family-wise error rate was calculated and corrected for multiple comparisons following the Holm-Bonferroni method with the p-adjust function in R (R Core Team 2020). We first examined the consistency of assigning individuals to a subclade, following the approach of Malinsky et al. (2018), by testing whether individuals assigned to the same subclade always share more derived alleles with each other than with individuals from other subclades. For each individual A, we calculated $D(A, G_1; G_2, O)$, where A is assigned to one subclade and $G_1$ and $G_2$ to another. Values greater than 0 indicate that allele sharing is greater within groups compared with across groups.

To characterize the relationships among subclades, we examined the consistency of subclade monophyly in genomic windows. We used the set of window-based trees previously generated for the ASTRAL analysis and
checked whether a subclade was monophyletic in each window with the check_monophyly function in the ETE Toolkit v.3 (Huerta-Cepas et al. 2016). Next, a cloudogram of ML trees was visualized. Window-based trees were modified using R packages phytool v.1.0-3 (Revell 2012) and treeio v.1.20.0 (Wang et al. 2020) to force branch lengths into an ultrametric structure for visualization purposes. A cloudogram was generated with ggdensitree from the R package ggtree v.3.4.0 (Yu 2020). We next used $D_{min}$, another variation of the *D*-statistic (Malinsky et al. 2018), to assess whether allele sharing among subclades is consistent with a tree-like structure. $D_{min}$ yields a conservative inference of gene flow within a trio by providing the minimal absolute value of *D* across all possible arrangements within it. A significantly positive score signifies allele sharing that is inconsistent with a single species tree. To identify and visualize the genetic structure between subclades, we performed a principal component analysis (PCA) on a reduced data set. This data set was pruned for linkage disequilibrium by randomly selecting SNPs that are at least 10 kb apart among all the biallelic SNPs, with <10% missing data. The PCA was performed with SNPRelate v.1.20.1 (Zheng et al. 2012).

Finally, to quantify and visualize the relationship among subclades along the genome, we used topology weighting by iterative sampling of subtrees with *Twisst* v.0.2 (Martin and Van Belleghem 2017). Using ML topologies inferred on genomic windows, this method considers all possible topologies and quantifies the contribution of each to the full tree, enabling to locate genomic regions that are associated with certain topologies. We focused on the 3 largest subclades and reduced the number of possible topologies by including only the subclades *T. extensa*, *T. utriculata*, and *T. fasciculata*. For this analysis, subtrees were constructed on windows of 50 SNPs following the strategy of Martin and Van Belleghem (2017). Individual VCFs for each window were produced by partitioning the original VCF with the biostar497922 script from Jvarkit (Lindenbaum 2015) and ML trees were inferred for each window as described previously. A summary and visualization of all topologies along genomic coordinates were produced with a slightly modified version of the plot_twisst.R script, which excludes windows with high rates of missing data in the computation of the average regional weightings (Martin and Van Belleghem 2017; script is available at https://github.com/giyany/Tillandsia_Phylo_hybridization/blob/main/notebooks/plot_twisst_mod.R).

## Characterizing Hybridization Events

To quantify the rates of hybridization between all species in our data set, we obtained a genome-wide signal of hybridization using the original implementation of the *D*-statistic and estimates of admixture fraction *f* (henceforth $f_4$-ratio) in *Dsuite* v.0.5r45 (Malinsky

et al. 2021). Given a certain level of uncertainty regarding the true relationship between species, we set no a priori knowledge of taxon relationships. Instead, *Dsuite* ordered each trio so that the BBAA pattern is more common, principally to focus on topologies with minimal discordant patterns. For all *D*-statistics-related analyses next, we again used the full VCF with *T. complanata* as an outgroup. We repeated the analysis separately for each reference chromosome and computed the *D*-statistic for each, obtaining *P*-values using a jackknife procedure with windows of 150 SNPs.

Because groups that are involved in hybridization may share branches on a phylogenetic tree, a single hybridization event can present multiple correlated instances of a significantly elevated *D*-statistic. This is especially expected when gene flow involves ancestral lineages, affecting internal branches of a phylogenetic tree. We used *Dsuite* to calculate the *f*-branch metric [henceforth $f_b(C)$], an estimator developed to create a summary of gene flow events with minimized correlation (Malinsky et al. 2021). The $f_b(C)$ results invited several hypotheses regarding hybridization events, which we examined by inferring a species network under a maximum pseudo-likelihood approach using PhyloNet v.3.8.0 (Than et al. 2008; Yu and Nakhleh 2015; Cao et al. 2019). For PhyloNet, we reduced our sampling to 1 outgroup and 18 ingroup taxa to include representatives from each highly supported subclade, and then inferred ML trees on nonoverlapping windows of 10 kb as previously described. We inferred networks that specified between 0 and 5 reticulations and repeated the search 5 times for each network, finally picking the network with the highest pseudo-likelihood score in each search.

## Characterizing Hybridization on Chromosome 18

Our previous calculation of the *D*-statistics for each reference chromosome (see previously) pointed at significantly elevated values on chromosome 18. Gene flow particularly involved the species *T. achyrostachys* (*T. punctulata* subclade) and most species in the CAM subclade *T. fasciculata* (see the “Results” section). Species in the *T. punctulata* subclade putatively possess an intermediate C3-CAM photosynthetic syndrome, expressing intermediary values of stable carbon isotope ratios ($\delta^{13}$C; Supplementary Tables S2; but see also Messerschmid et al. 2021). Within our sampling of this subclade, *T. achyrostachys* was documented to express the strongest CAM phenotype ($\delta^{13}$C of −14.7; Crayn et al. 2015). Additionally, all known species in the *T. fasciculata* subclade are putatively strong, constitutively CAM plants. We wished to investigate the possibility that gene flow introduced advantageous genes related to shifts in metabolic syndromes. Hence, we inferred the signature of hybridization along chromosome 18 to locate highly admixed loci. We focused on possible hybridization between *T. punctulata*, *T. butzii*, and *T. achyrostachys* (cor-
responding to the *T. punctulata* subclade, the *T. fasciculata* subclade, and the *T. punctulata* subclade, respectively). This allowed us to account for phylogenomic relatedness between *T. achyrostachys* and *T. punctulata* and for similar photosynthetic syndromes between *T. achyrostachys* and *T. butzii*. We used the Dinvestigate function from *Dsuite* on this species trio, performing analysis on windows of 50 SNPs with a step size of 10, obtained as previously described. The *D*-statistic itself shows large variance when applied to genomic windows (Martin et al. 2015); hence, we used *f<sub>dM</sub>* a statistic designed to investigate hybridization in small windows, which also accounts for allele sharing across all possible taxon arrangements in a trio (Malinsky et al. 2015; Martin et al. 2015). We identified regions exhibiting *D*-statistic values exceeding the 95% quantile of the distribution and inspected the genes annotated in regions of high admixture, specifically their functional annotation.

## RESULTS

### Read Mapping and Variant Calling

After removing samples with low coverage, those sharing high kinship coefficients, and the *T. zoquensis* accession (see the "Materials and Methods" section), we retained 64 *Tillandsia* accessions, corresponding to 35 recognized species. The average number of reads retained per accession was $5.1 \times 10^7$ (range $4.7 \times 10^6$–$1.2 \times 10^8$, standard deviation [SD] = $2.5 \times 10^7$) and average mapping rates were 89.6% (range 69.4–97.5%, SD = 6.3), with slightly higher rates for members of the *T. fasciculata* subclade and slightly lower for the SA subclades (97.5% and 89.6%, respectively). Differences in mapping rates between the 6 main subclades (Fig. 2) were, however, not significant (Kruskal–Wallis test, *P* = 0.22), suggesting no or limited biases toward the reference genome. An average coverage of $13.8 \times$ (range $4.1 \times$–$35.0 \times$, SD = $6.5 \times$) was obtained for the samples retained further. After variant calling and filtering, we retained 2,162,143 high-quality SNPs.

### Phylogenomic Inference Contradicts Clade Monophyly

The concatenated matrix was partitioned into 14,392 regions, whereas the coalescent-based analyses were performed on 15,791 genomic windows. The species tree (Fig. 2) and the concatenated ML tree (Supplementary Fig. S1) yielded inconsistent results regarding the main clades and relationships. In both trees, the Central American “K” subclades did not form a monophyletic group, contrary to among-subclade relationships previously reported for this subgenus (Fig. 1; Barfuss et al. 2016; Granados Mendoza et al. 2017; Rivera 2019; Vera-Paz et al. 2023). Specifically, our results are incongruent with a monophyly of the previous inferred “subclade K.2.” Instead, the *T. utriculata* subclade (K.2.1) was recovered as a sister to one or both of the South

American subclades (in the ML tree or the species tree, respectively). The assignment of species to subclades and all other relationships between subclades overall remained congruent with previous phylogenies (Granados Mendoza et al. 2017; Vera-Paz et al. 2023). Several, relatively minor within-clade topological differences appeared among trees inferred with different methods, especially within the *T. fasciculata* subclade. Relationships between species in this subclade were highly supported yet coupled with high levels of gene-tree incongruence, suggesting substantial allele sharing and high rates of gene flow within this group. Although the South American taxa in our sampling formed a monophyletic group in the species tree (Fig. 2), in the ML tree (Supplementary Fig. S1) these were separated into 2 subclades: one consisted of endemic Peruvian species and the widespread species *T. adpressiflora*, whereas a second subclade contained *T. marnierlapostollei* and *T. mima*. The latter was retrieved as paragraphetic to the *T. utriculata* subclade and all remaining South American taxa. As the subclades we retrieved do not correspond to the previously inferred *T. paniculata* and *T. secunda* clades (Vera-Paz et al. 2023), we renamed them *T. extensa* and *T. mima* subclades, respectively (Supplementary Table S1). We additionally renamed the Mexican “K” subclades to reflect the new assignment (Fig. 2; see the “Materials and Methods” section). Gene-tree discordance was widespread within the data set, affecting both deep and shallow nodes (Fig. 2; Supplementary Fig. S3). The relationships between the South American species and the *T. utriculata* subclade were particularly characterized by short internode distances and many alternative topologies. Frequently, the majority of inferred gene-tree topologies were discordant with a single main topology: For example, in the node preceding the separation of the SA subclades, only 3148 (19.9%) of all gene windows supported the main topology (Fig. 2; Supplementary Fig. S3). Similar levels of discordance characterized the internal nodes within the *T. fasciculata* subclade and high levels of discordance were also found within the *T. punctulata* subclade (Fig. 2; Supplementary Fig. S3).

ML trees constructed separately for concatenated matrices of each reference chromosome retrieved many different topologies: Solely considering relationships between subclades, we recognized 10 different topologies among the different 25 trees (see Supplementary File 1). The placement of the *T. mima* subclade showed the greatest incongruence between trees, as well as the placement of *T. adpressiflora*. In 3 trees, the relationships between subclades were similar to the topology retrieved from the whole-genome concatenated matrix. The South American species formed a monophyly in 6 chromosome trees and “clade K” was recovered as monophyletic in 2. Although each chromosome tree contains high amounts of gene-tree discordance, the abundance of different topologies implies that several evolutionary histories can be traced along the genome of
*Tillandsia*, as disparate genomic processes contribute to the differences between gene trees and the true species relationships. Trees inferred on data sets with different MAF filtering thresholds did not substantially differ in topology, whereas gene-tree discordance was minimally affected (Supplementary File 3). Branch lengths were, however, affected, as minor alleles contribute to long branches of specific species. We again retrieved inconsistency in the placement of the *T. mima* subclade when using different filtering parameters, indicating that the subclade cannot be confidently placed in our current data set.

## Lack of Monophyly and Deviations from Tree-Like Structure

The multitude of tree topologies along the genome led us to hypothesize that a single bifurcating tree misrepresents the true relationships among species in the subgenus *Tillandsia*, in particular affecting “clade K” (Supplementary Fig. S4). To investigate deviations from a tree-like structure, we analyzed the patterns of allele sharing between species and between subclades. Assuming no interspecific gene flow, allele sharing is expected to be consistent with the main tree topology, whereas asymmetrical allele sharing indicates deviations from that assumption.

We first used the *D*-statistic framework to test the robustness of the assignment of species to subclades. We calculated allele sharing between species, utilizing the *D*-statistic’s allele frequency estimates. We found high support for subclade monophyly, as species assigned to the same subclade shared more alleles with other species within the subclade than with species from other subclades. However, we observed low support for subclade monophyly based on the analyses in genomic windows performed with ETE: 52.4%, 34.4%, and 14.11% of the windows supported monophyly for subclades *T. guatemalensis*, *T. fasciculata*, and *T. punctulata*, respectively. A notable exception was observed for the subclade *T. utriculata*, for which 86.6% of the windows were consistent with a monophyly.

We next characterized allele sharing between subclades by computing the *D*<sub>min</sub>-statistic on a total of 7141 trios (Supplementary Fig. S5a). We found widespread deviations from a tree-like structure, mostly driven by allele sharing between Central American taxa (previous “clade K”): *D*<sub>min</sub> values were significantly elevated in 3650 comparisons (*P* < 0.05) and more than 95% of those were highly significant (*P* < 0.01). The rate of significant *D* values was highest for comparisons involving accessions assigned to the *T. utriculata* subclade (58.7% of the trios), followed by *T. punctulata*, *T. guatemalensis*, and *T. fasciculata* (57.9%, 56.1%, and 51.4%, respectively). We also employed a multivariate technique with PCA to investigate the interspecific genetic structure (Supplementary Fig. S5b). After distance pruning and removing SNPs with <10% missing data, we retained 16,204 SNPs. We found a consistent interspe-

cific genetic structure, separating the South American subclades and the *T. utriculata* subclade, whereas the remainder of the Central American subclades clustered densely—reflecting the distinct subclades found in tree-based and *D*-statistics approaches.

Finally, to obtain a better understanding of how subclade relationships vary along the genome, we considered 3 distinct tree topologies for topology weighting with *Twisst* (Supplementary Fig. S6): The most frequent topology was congruent with the one recovered in the species tree, which placed the *T. extensa* subclade as sister to the *T. utriculata* subclade in 42% of the genomic windows. A second topology, with *T. fasciculata* and *T. utriculata* as sister subclades, was recovered from 35% of the genomic windows, whereas a third topology appeared in 23%. The 3 topologies were broadly equally distributed along the genome—however, in some chromosomes (e.g., chromosome 4; Supplementary Fig. S6), the first topology was more frequent within centromeric regions. Given that topology weighing remains a descriptive method, it does not allow explicit testing for hybridization or ILS. Regardless, the prevalence of the main topology in regions of low recombination rates and reduced genic density suggests that it represents the backbone phylogeny, whereas other topologies are likely the result of gene flow or deep coalescence (Carneiro et al. 2009; Christmas et al. 2021).

## Correlated and Widespread Gene Flow Events

*D*-statistics results for all possible trios indicated that all subgenus *Tillandsia* species in our data set were involved in potential hybridization (Supplementary Fig. S7). Out of 7141 tests in total, 4331 returned significantly elevated values, with *D* values ranging between 0.021 and 0.581. The signal was not localized to a specific chromosome according to separate calculations on each chromosome (Supplementary File 2). A prominent signal revealed gene flow between species in the *T. utriculata* subclade and all other Central American subclades. Notable gene flow signals were also found within the *T. fasciculata* subclade. Overall, *f*<sub>4</sub>-ratio scores ranged between 0.0017 and 0.357, but for most hybridizations, the proportion of the genome involved was estimated smaller than 10% (Supplementary Fig. S8). Larger parts of the genome were admixed within the main subclades: For example, *T. caput-medusae* was involved in hybridization with *T. butzii* and *T. fasciculata* involving ca. 31.7% of the genome.

Further analysis confirmed that past hybridization involved the ancestor of the *T. urticulata* subclade and an ancestor of the other 3 Central American clades. In analysis based on the *f*-branch metric, elevated *f*<sub>b</sub>(C) scores were assigned to events involving these subclades with an average of 6.3% (Fig. 3). Among all between-subclades *f*<sub>b</sub>(C) values, 246 (19.6%) were significantly elevated (*P* < 0.05), although most signifi-
<figure><img src="image_3.png" /><figcaption>FIGURE 3. Heatmap summarizing the statistic $f_b(C)$, where excess sharing of derived alleles is inferred between the branch of the tree on the y-axis and the species C on the x-axis. The ASTRAL species tree was used as input topology for the branch statistic. The matrix is colored according to $f_b(C)$ values and gray squares correspond to tests that are inconsistent with the ASTRAL phylogeny. Dots within the matrix denote a significant P-value, estimated using a block jackknife procedure and corrected for family-wise error rate. Colors correspond to the clades in Fig. 1.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 4. Pseudo-likelihood species networks inferred with PhyloNet for 0–3 reticulations (network presented was scored with the highest log probability). Curved branches indicate reticulation events. Numbers next to curved branches indicate inheritance probabilities for each event. Colors correspond to the clades in Fig. 1.</figcaption></figure>

cant values (412, 62.6%) occurred within subclades. In a network analysis on a total of 14,535 genomic windows, allowing for 4 or more reticulations resulted in events involving the outgroup or fewer than 4 reticulations reported, so we present here results for up to 3 reticulations. Otherwise, the results repeatedly indicated the involvement of the *T. utriculata* subclade in hybridization with the *T. extensa* and other Central American subclades (Fig. 4). Notably, higher weights were assigned to the former hybridization. Additional gene flow occurred within subclades. Overall, these findings confirm the extensive violations of a tree-like structure revealed in the previous part of the analyses and suggest that interspecific gene flow occurred frequently throughout the evolutionary history of the radiation.

### Local Footprint of Hybridization between CAM-Expressing Species

Elevated *D*-statistics, signifying hybridization, generally affected all chromosomes (File 2). In contrast, a significantly elevated, but localized *D*-statistic was found on chromosome 18 involving *T. achyrostachys* (*T. punctulata* subclade)

and most of the species in the CAM subclade *T. fasciculata* (Supplementary File 2). To further investigate this highly localized signal, *f*<sub>dM</sub> statistics were calculated on a total of 9185 windows between *T. punctulata*, *T. butzii*, and *T. achyrostachys* (P1, P2, and P3). The first and last species are assigned to the *T. punctulata* subclade, whereas *T. butzii* is a member of the *T. fasciculata* subclade. Across all chromosomes, mean values were negative at an average of −0.043, fitting the expectation for higher rates of gene flow between closely related species (Fig. 5a). Chromosome 18 was, however, characterized by 2 regions of high positive values, which were also evident with topology weighting (Fig. 5b and c). Considering centromere localization, we found that the loci exhibiting high *f*<sub>dM</sub> values coincide with low genic density and contain 194 genes (Fig. 5d). In a detailed survey of gene annotations, we found that this set of genes was potentially associated with a variety of stress and metabolic functions. For example, serine/threonine-protein kinase prpf4B is known to have a role in pre-mRNA splicing in yeast and humans (Eckert et al. 2016), and was found to be associated with stress response in millet (Parvathi et al. 2019). Another example is the S-adenosyl-L-
<figure><img src="image_5.png" /></figure>

FIGURE 5. Signature of introgression and topology weighting on chromosome 18 between *T. punctulata*, *T. butzii*, and *T. achyrostachys* (P1, P2, and P3; *T. complanata* was used as the outgroup). (a) Values of *f*<sub>dM</sub>-statistic and gene content for each chromosome calculated between *T. punctulata*, *T. butzii*, and *T. achyrostachys*. (b) Topology weighting by iterative sampling of subtrees in genomic windows of 50 SNPs using *Twisst*. Colors represent the frequency of each topology in (c) along the chromosomal position in chromosome 18: White gaps indicate regions excluded due to high levels of missing data. (d) *f*<sub>dM</sub>-statistic (green dots and smoothed green line, scale on the right) calculated in genomic windows. Analysis was performed on windows of 50 SNPs with a step size of 10. Shared variation is quantified in positive values when shared between P2 and P3 and as negative values when shared between P3 and P1. Gene content as the number of predicted genes per 1 MB window is in lavender with scale on the left.

methionine-dependent methyltransferase superfamily protein (SAM-Mtase), a key enzyme in plant metabolic pathways such as the phenylpropanoid and flavonoid pathway (Joshi and Chiang 1998; Sistla and Rao 2004).

DISCUSSION

Hybridization is extensively studied for its various roles in plant speciation. Using a multispecies genomic data set and phylogenomic approaches, we present the subgenus *Tillandsia* as a striking example of a highly reticulated radiation. The radiation likely proceeded in the presence of rampant interspecific gene flow, rather than by refined reproductive barriers (Loiseau et al. 2019; Koch et al. 2022; Vera-Paz 2023). An abundance of sequencing data tends to generate high branch support in both concatenated alignments and a coalescent-based tree-building method (Salichos and Rokas 2013;

Pease et al. 2016), yet a detailed investigation in the subgenus *Tillandsia* revealed that a bifurcating tree offers an incomplete picture of the true relationships between species. This work thus provides a compelling example to the mounting evidence of the presence, and possibly contribution, of gene flow in species diversification (Nosil 2008; Seehausen et al. 2014; Keller et al. 2013; Arnold et al. 2016; Filiault et al. 2018; Stankowski et al. 2019; Zhang et al. 2021; Slovak et al. 2023).

In stark contrast to previously inferred phylogenies, we retrieved the Central American "clade K" as polyphyletic. First inferred based on morphological characters and later on a limited number of genomic markers, the generic and subgeneric classification of the subfamily Tillandsioideae in general and of the genus *Tillandsia* in particular shifted throughout decades of phylogenetic research. Overall, phylogenies based on plastid sequences and on relatively few nuclear genomic regions offered little resolution for shallow phyloge-
netics in this young genus, particularly within Central American taxa (Barfuss et al. 2016; Pinzón et al. 2016; Granados Mendoza et al. 2017). Recently, a full plastome phylogeny greatly improved the phylogenomic resolution, whereas statistical support remained low in shallow nodes (Vera-Paz et al. 2023). We suggest that incongruence between *Tillandsia* phylogenies reflects both the local scarcity of genomic divergence within this young radiation and the complicating role of hybridization. Moreover, considering the prevalence of recent interspecific gene flow, the monophyly of the Central American subclades in plastid phylogenies seems particularly driven by chloroplast capture (Pham et al. 2017). Beyond insight on subclade relationship, our analyses produced several novel insights regarding cryptic species: For example, our inference suggests that *T. mima* and *T. mima* var. *chiletensis* are distinct species, despite morphological similarities. Similarly, *T. fuchsii* and *T. fuchsii* var. *stephanii* did not form a monophyletic group, raising a need to revise their taxonomic status (Fig. 2; Supplementary Fig. S1).

Hybridization played a central role in the evolution of Central American *Tillandsia*, reflected in departures from tree structure encompassing all clades, subclades, and species. *D*-statistics revealed both recent and ancient gene flows, whereas allele sharing did not compromise species or subclade boundaries. Notably, the obtained *f*<sub>b</sub>(*C*) ranges were in general higher than those previously inferred for ancestral introgressions between species of Malawi cichlids (Malinsky et al. 2018) or of hares (Ferreira et al. 2021), but comparable to other rapid radiations (De-Kayne et al. 2022; Slovak et al. 2023). The interplay of genomic signals and our limited sampling complicate inference on the timing of gene flow events, yet their high prevalence suggests occurrence during both deeper and recent history.

High rates of gene flow are compatible with *Tillandsia* ecology and evolution: *Tillandsia* is known to produce natural hybrids (Luther 1985; Till 2000; Koch et al. 2022) and bears copious seed adapted to wind dispersal, which may facilitate high rates of inter- and intraspecific gene flow (Mondragon and Calvo-Irabien 2006; Victoriano-Romero et al. 2017). Previous studies proposed that the South American ancestor of *Tillandsia* subgenus *Tillandsia* colonized North and Central America in a long-distance dispersal event 4.86 Ma (Barfuss et al. 2016; Granados Mendoza et al. 2017; Vera-Paz et al. 2023), and our results offer intriguing hypotheses on multiple, likely polytopic origins of Central American *Tillandsia*. Instead of a single dispersal event, South American ancestors may have dispersed in several migration events into Central America. Strong population bottlenecks, producing phylogenies with relatively long internal branches, were followed by concomitant episodes of gene admixture and increased isolation. It is further possible that gene flow increased allelic diversity in founder populations, thus fueling adaptation (see next).

Apart from hybridization, other molecular processes can contribute to violations of a strictly bifurcating species tree, such as ILS, paralogy, and gene duplication and loss (Galtier and Daubin 2008; Edwards 2009; Smith et al. 2015). Previous studies on *Tillandsia* subgenus *Tillandsia* indeed found evidence for changes in population sizes and elevated rates of gene duplication and loss, specifically associated with photosynthetic syndrome shifts (de La Harpe et al. 2020; Yardeni et al. 2022; Groot Crego et al. 2024). We suggest that there is strong evidence for ancestral hybridization as the source of discordance between *T. utriculata* and other Central American subclades: *f*-statistic and *D*-statistic tests are robust to the presence of ILS, as is PhyloNet (Martin et al. 2015; Wen et al. 2018; Malinsky et al. 2021), and ancestral population structure is unlikely to segregate through the demographic events that accompanied *Tillandsia*'s dispersal into Central America. Regardless, the generality of the obtained signal could be influenced by our limited sampling; a study employing a wider sampling could further infer whether hybridization involved related or ancient lineages, while considering disparate molecular processes.

Genetic variation introduced through hybridization can be maladaptive, neutral, or adaptive (Zeberg and Pääbo 2020; Moran et al. 2021; Wong et al. 2022). Although its consequences remain difficult to predict and negative consequences may not be observable for long periods, an adaptive role of introgression has been demonstrated in numerous animal and plant taxa (Dasmahapatra et al. 2012; Lamichhaney et al. 2015; Arnold et al. 2016; Taylor and Larson 2019; Leroy et al. 2020). In this study, we did not directly test for the hypothesis that hybridization facilitated adaptation. However, we have identified 2 candidate regions on chromosome 18, with genomic signatures of hybridization occurring between species that share similar photosynthetic syndromes. Interestingly, the identified regions coincide with low gene density and one of them occurs at a chromosome edge. These findings echo the general expectation regarding the theory, posing that higher rates of hybridization correlate with regions of low gene density, because hybridized fragments in gene-rich regions are more likely to be detrimental (Barton and Bengtsson 1986; Sankararaman et al. 2014; Martin and Jiggins 2017). The incomplete taxon sampling in our current study limits our ability to draw generalities about the extent of adaptive hybridization and introgression in *Tillandsia*: However, considering the rapid accumulation of morphological and physiological disparity, we may hypothesize that hybridization in *Tillandsia* contributed to its success in a wide range of habitats, potentially introducing novel alleles through porous species boundaries and driving this evolutionary radiation. Further efforts will be needed to identify additional potential adaptive regions and uncover the key genes within.
We used whole-genome sequencing to deeply investigate the phylogenomics of a remarkable *Tillandsia* radiation. Recent research addressed questions on the genus' evolutionary history, simultaneously expanding bromeliad and *Tillandsia* genomic resources and facilitating a range of evolutionary analyses. Ultimately, the complex history of *Tillandsia* remains elusive, calling for further investigations to uncover the interplay of processes that drove this rapid diversification. Future studies employing both a wider sampling and a deeper genomic coverage are warranted to characterize the genomic properties associated with diversification and explore the prevalence and consequences of hybridization across different evolutionary scales.

## SUPPLEMENTARY MATERIAL

Data available from the Dryad Digital Repository: https://datadryad.org/dataset/doi:10.5061/dryad.x0k6djhrs

## FUNDING

This work was supported with funding from the Christian lexer professorship start-up BE772002 at the University of Vienna and in part by the Austrian Science Fund (FWF) [https://doi.org/10.55776/P35275] to O.P.

## ACKNOWLEDGMENTS

This paper is dedicated to the memory of Prof. Christian lexer, who conceived and guided this study, but passed away before its completion. The analyses took advantage of the Vienna Scientific Cluster. We thank Magnus Nordborg, James Pease, and Hannes Svardal for insightful discussions and advice. We thank the members of Swiss SNSF Sinergia project CRSII3_147630 and Juan P. Pinzón for accession sampling. We thank 3 anonymous reviewers for their insightful and detailed comments that greatly assisted to improve the manuscript.

## DATA AVAILABILITY

Alignment additional data can be found in the Dryad data repository: https://datadryad.org/dataset/doi:10.5061/dryad.x0k6djhrs

Additional information, scripts, and protocols to replicate the analyses performed in this study are available from: https://github.com/giyany/TillandsiaPhyl o/

All raw sequence reads are available at NCBI-SRA under BioProject PRJNA1052467.

## CONFLICT OF INTEREST

None declared.

## REFERENCES

* Abbott R., Albach D., Ansell S., Arntzen J.W., Baird S.J.E., Bierce N., Boughman J., Brelsford A., Buerkle C.A., Buggs R., Butlin R.K., Dieckmann U., Eroukhmanoff F., Grill A., Cahan S.H., Hermansen J.S., Hewitt G., Hudson A.G., Jiggins C., Jones J., Keller B., Marczewski T., Mallet J., Martinez-Rodriguez P., Möst M., Mullen S., Nichols R., Nolte A.W., Parisod C., Pfennig K., Rice A.M., Ritchie M.G., Seifert B., Smadja C.M., Stelkens R., Szymura J.M., Väinölä R., Wolf J.B.W., Zinner D. 2013. Hybridization and speciation. *J. Evol. Biol.* 26(2):229–246.

* Andrews S. 2010. FastQC: a quality control tool for high throughput sequence data. Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

* Arnold B.J., Lahner B., DaCosta J.M., Weisman C.M., Hollister J.D., Salt D.E., Bomblies K., Yant L. 2016. Borrowed alleles and convergence in serpentine adaptation. *Proc. Natl. Acad. Sci. USA* 113(29):8320–8325.

* Barfuss M.H.J., Samuel R., Till W., Stuessy T.F. 2005. Phylogenetic relationships in subfamily Tillandsioideae (Bromeliaceae) based on DNA sequence data from seven plastid regions. *Am. J. Bot.* 92(2):337–351.

* Barfuss M.H.J., Till W., Leme E.M.C., Pinzón J.P., Manzanares J.M., Halbritter H., Samuel R., Brown G.K. 2016. Taxonomic revision of Bromeliaceae subfam. Tillandsioideae based on a multi-locus DNA sequence phylogeny and morphology. *Phytotaxa* 279(1):1–97.

* Barnett D.W., Garrison E.K., Quinlan A.R., Strömberg M.P., Marth G.T. 2011. BamTools: a C++ API and toolkit for analyzing and managing BAM files. *Bioinformatics* 27(12):1691–1692.

* Barton N., Bengtsson B.O. 1986. The barrier to genetic exchange between hybridising populations. *Heredity.* 57(3):357–376.

* Benzing D.H. 2000. Bromeliaceae: profile of an adaptive radiation. *New York, USA: Cambridge University Press, .*

* Boschman L.M., Condamine F.L. 2022. Mountain radiations are not only rapid and recent: Ancient diversification of South American frog and lizard families related to Paleogene Andean orogenic and Cenozoic climate variations. *Global Planet. Change* 208:103704.

* Cao Z., Zhu J., Nakhleh L. 2019. Empirical performance of tree-based inference of phylogenetic networks 19th International Workshop Algorithms Bioinformatics (WABI 2019). 143:21:1 – 21:13 Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik.

* Carneiro M., Ferrand N., Nachman M.W. 2009. Recombination and speciation: loci near centromeres are more differentiated than loci near telomeres between subspecies of the European rabbit (*Oryctolagus cuniculus*). *Genetics* 181(2):593–606.

* Chernomor O., von Haeseler A., Minh B.Q. 2016. Terrace aware data structure for phylogenomic inference from supermatrices. *Syst. Biol.* 65(6):997–1008.

* Chew T., De Luna E., González D. 2010. Phylogenetic relationships of the pseudobulbous *Tillandsia* species (Bromeliaceae) inferred from cladistic analyses of ITS 2, 5.8 S ribosomal RNA gene, and ETS sequences. *Syst. Bot.* 35(1):86–95.

* Christmas M.J., Jones J.C., Olsson A., Wallerman O., Bunikis I., Kierczak M., Peona V., Whitley K.M., Larva T., Suh A., Miller-Struttmann N.E., Geib J.C., Webster M.T. 2021. Genetic barriers to historical gene flow between cryptic species of alpine bumblebees revealed by comparative population genomics. *Mol. Biol. Evol.* 38(8):3126–3143.

* Cingolani P., Platts A., Wang L.L., Coon M., Nguyen T., Wang L., Land S.J., Lu X., Ruden D.M. 2012. A program for annotating and predicting the effects of single nucleotide polymorphisms, *SnpEff. Fly* 6(2):80–92.

* Crayn D.M., Winter K., Schulte K., Smith J.A.C. 2015. Photosynthetic pathways in Bromeliaceae: phylogenetic and ecological significance of CAM and C3 based on carbon isotope ratios for 1893 species. *Bot. J. Linn. Soc.* 178(2):169–221.

* Danecek P., Bonfield J.K., Liddle J., Marshall J., Ohan V., Pollard M.O., Whitwham A., Keane T., McCarthy S.A., Davies R.M., Li H. 2021. Twelve years of SAMtools and BCFtools. *Gigascience* 10(2):giab008.

* Dasmahapatra K.K., Walters J.R., Briscoe A.D., Davey J.W., Whibley A., Nadeau N.J., Zimin A.V., Hughes D.S.T., Ferguson L.C., Martin S.H., Salazar C., Lewis J.J., Adler S., Ahn S.-J., Baker D.A., Baxter S.W., Chamberlain N.L., Chauhan R., Counterman B.A., Dal-
may T., Gilbert L.E., Gordon K., Heckel D.G., Hines H.M., Hoff K.J., Holland P.W.H., Jacquin-Joly E., Jiggins F.M., Jones R.T., Kapan D.D., Kersey P., Lamas G., Lawson D., Mapleson D., Maroja L.S., Martin A., Moxon S., Palmer W.J., Papa R., Papanicolaou A., Pauchet Y., Ray D.A., Rosser N., Salzberg S.L., Supple M.A., Surridge A., Tenger-Trolander A., Vogel H., Wilkinson P.A., Wilson D., Yorke J.A., Yuan F., Balmuth A.L., Eland C., Gharbi K., Thomson M., Gibbs R.A., Han Y., Jayaseelan J.C., Kovar C., Mathew T., Muzny D.M., Ongeri F., Pu L.-L., Qu J., Thornton R.L., Wormaly K.C., Wu Y.-Q., Linares M., Blaxter M.L., ffrench-Constant R.H., Joron M., Kronforst M.R., Mullen S.P., Reed R.D., Scherer S.E., Richards S., Mallet J., Owen McMillan W., Jiggins C.D.; The Heliconius Genome Consortium. 2012. Butterfly genome reveals promiscuous exchange of mimicry adaptations among species. Nature 487:94-98.

de La, Harpe M., Paris M., Hess J., Barfuss M.H.J., Serrano-Serrano M.L., Ghatak A., Chaturvedi P., Weckwerth W., Till W., Salamin N., Wai C.M., Ming R., Larer C. 2020. Genomic footprints of repeated evolution of CAM photosynthesis in a Neotropical species radiation. Plant Cell Environ. 43:2987-3001.

Debray K., Marie-Magdelaine J., Ruttink T., Clotault J., Foucher F., Malécot V. 2019. Identification and assessment of variable single-copy orthologous (SCO) nuclear loci for low-level phylogenomics: A case study in the genus Rosa (Rosaceae). BMC Evol. Biol. 19(1):1-19.

De-Kayne R., Selz O.M., Marques D.A., Frei D., Seehausen O., Feulner P.G.D. 2022. Genomic architecture of adaptive radiation and hybridization in Alpine whitefish. Nat. Commun. 13(1):4479.

DePristo M.A., Banks E., Poplin R., Garimella K.V., Maguire J.R., Hartl C., Philippakis A.A., del Angel G., Rivas M.A., Hanna M., McKenna A., Fennell T.J., Kernytsky A.M., Sivachenko A.Y., Cibulskis K., Gabriel S.B., Altshuler D., Daly M.J. 2011. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat. Genet. 43(5):491-498.

Doyle J.J., Doyle J.L. 1987. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem. Bull. 19:11-15.

Drummond C.S., Eastwood R.J., Miotto S.T.S., Hughes C.E. 2012. Multiple continental radiations and correlates of diversification in *Lupinus* (Leguminosae): testing for key innovation with incomplete taxon sampling. Syst. Biol. 61(3):443-460.

Durand E.Y., Patterson N., Reich D., Slatkin M. 2011. Testing for ancient admixture between closely related populations. Mol. Biol. Evol. 28(8):2239-2252.

Eckert D., Andrée N., Razanau A., Zock-Emmenthal S., Lützelberger M., Plath S., Schmidt H., Guerra-Moreno A., Cozzuto L., Ayté J., Käufer N.F. 2016. Prp4 kinase grants the license to splice: control of weak splice sites during spliceosome activation. PLoS Genet. 12(1):e1005768.

Edelman N.B., Frandsen P.B., Miyagi M., Clavijo B., Davey J., Dikow R.B., García-Accinelli G., Van Belleghem S.M., Patterson N., Neafsey D.E., Challis R., Kumar S., Moreira G.R.P., Salazar C., Chouteau M., Counterman B.A., Papa R., Blaxter M., Reed R.D., Dasmahapatra K.K., Kronforst M., Joron M., Jiggins C.D., McMillan W.O., Di Palma F., Blumberg A.J., Wakeley J., Jaffe D., Mallet J. 2019. Genomic architecture and introgression shape a butterfly radiation. Science 366(6465):594-599.

Edwards S.V. 2009. Is a new and general theory of molecular systematics merging? Evolution 63(1):1-19.

Espejo-Serna A., López-Ferrari A.R., Ramírez-Morillo I., Holst B.K., Luther H.E., Till W. 2004. Checklist of Mexican Bromeliaceae with notes on species distribution and levels of endemism. Selbyana 25:33-86.

Ferreira M.S., Jones M.R., Callahan C.M., Farelo L., Tolesa Z., Suchentrunk F., Boursot P., Mills L.S., Alves P.C., Good J.M., Melo-Ferreira J. 2021. The legacy of recurrent introgression during the radiation of hares. Syst. Biol. 70(3):593-607.

Filiault D.L., Ballerini E.S., Mandáková T., Aköz G., Derieg N.J., Schmutz J., Jenkins J., Grimwood J., Shu S., Hayes R.D., Hellsten U., Barry K., Yan J., Mihaltcheva S., Karafiátová M., Nizhynska V., Kramer E.M., Lysak M.A., Hodges S.A., Nordborg M. 2018. The *Aquilegia* genome provides insight into adaptive radiation and re-

veals an extraordinarily polymorphic chromosome with a unique history. eLife 7:e36426.

Galtier N., Daubin V. 2008. Dealing with incongruence in phylogenomic analyses. Philos. Trans. R. Soc. B Biol. Sci. 363(1512):4023-4029.

Giarla T.C., Esselstyn J.A. 2015. The challenges of resolving a rapid, recent radiation: empirical and simulated phylogenomics of Philippine shrews. Syst. Biol. 64(5):727-740.

Givnish T.J. 2015. Adaptive radiation versus "radiation" and "explosive diversification": why conceptual distinctions are fundamental to understanding evolution. New Phytol. 207(2):297-303.

Givnish T.J., Barfuss M.H.J., Ee B.V., Riina R., Schulte K., Horres R., Gonsiska P.A., Jabaily R.S., Crayn D.M., Smith J.A.C., Winter K., Brown G.K., Evans T.M., Holst B.K., Luther H., Till W., Zizka G., Berry P.E., Sytsma K.J. 2011. Phylogeny, adaptive radiation, and historical biogeography in Bromeliaceae: insights from an eight-locus plastid phylogeny. Am. J. Bot. 98(5):872-895.

Givnish T.J., Barfuss M.H.J., Ee B.V., Riina R., Schulte K., Horres R., Gonsiska P.A., Jabaily R.S., Crayn D.M., Smith J.A.C., Winter K., Brown G.K., Evans T.M., Holst B.K., Luther H., Till W., Zizka G., Berry P.E., Sytsma K.J. 2014. Adaptive radiation, correlated and contingent evolution, and net species diversification in Bromeliaceae. Mol. Phylogenet. Evol. 71:55-78.

Granados Mendoza C., Granados-Aguilar X., Donadío S., Salazar G.A., Flores-Cruz M., Hágsater E., Starr J.R., Ibarra-Manríquez G., Golovano-Martínez I., Magallón S. 2017. Geographic structure in two highly diverse lineages of *Tillandsia* (Bromeliaceae). Botany 95(7):641-651.

Green R.E., Krause J., Briggs A.W., Maricic T., Stenzel U., Kircher M., Patterson N., Li H., Zhai W., Fritz M.H.-Y., Hansen N.F., Durand E.Y., Malaspinas A.-S., Jensen J.D., Marques-Bonet T., Alkan C., Prüfer K., Meyer M., Burbano H.A., Good J.M., Schultz R., Aiximu-Petri A., Butthof A., Höber B., Höffner B., Siegemund M., Weihmann A., Nusbaum C., Lander E.S., Russ C., Novod N., Affourtit J., Egholm M., Verna C., Rudan P., Brajkovic D., Kucan Ž., Gušic I., Doronichev VB., Golovanova L.V., Lalueza-Fox C., de la Rasilla M., Fortea J., Rosas A., Schmitz R.W., Johnson P.L.F., Eichler E.E., Falush D., Birney E., Mulkalin J.C., Slatkin M., Nielsen R., Kelso J., Lachmann M., Reich D., Pääbo S. 2010. A draft sequence of the Neandertal genome. Science 328(5979):710-722.

Groot Crego C., Hess J., Yardeni G., de La Harpe M., Priemer C., Beclin F., Saadain S., Cauz-Santos L.A., Temsch E.M., Weiss-Schneeweiss H., Barfuss M.H.J., Till W., Weckwerth W., Heyduk K., Larer C., Paun O., Leroy T. 2024. CAM evolution is associated with gene family expansion in an explosive bromeliad radiation. Plant Cell 36(10):4109-4131.

Guo C., Luo Y., Gao L.M., Yi T.S., Li H.T., Yang J.B., Li D.Z. 2023. Phylogenomics and the flowering plant tree of life. J. Integr. Plant Biol. 65(2):299-323.

Harrison R.G., Larson E.L. 2014. Hybridization, introgression, and the nature of species boundaries. J. Hered. 105:795-809.

He Z., Li X., Yang M., Wang X., Zhong C., Duke N.C., Wu C., Shi S. 2019. Speciation with gene flow via cycles of isolation and migration: insights from multiple mangrove taxa. Natl. Sci. Rev. 6(2): 275-288.

Hibbins M.S., Gibson M.J., Hahn M.W. 2020. Determining the probability of hemiplasy in the presence of incomplete lineage sorting and introgression. eLife 9:e63753.

Hudlow W.R., Krieger R., Meusel M., Sehhat J.C., Timken M.D., Buoncristiani M.R. 2011. The NucleoSpin® DNA Clean-up XS kit for the concentration and purification of genomic DNA extracts: an alternative to microdialysis filtration. Forensic Sci. Int. Genet. 5(3):226-230.

Huerta-Cepas J., Serra F., Bork P. 2016. ETE 3: reconstruction, analysis, and visualization of phylogenomic data. Mol. Biol. Evol. 33(6):1635-1638.

Hughes C.E., Atchison G.W. 2015. The ubiquity of alpine plant radiations: from the Andes to the Hengduan Mountains. New Phytol. 207(2):275-282.

Hughes C.E., Nyffeler R., Linder H.P. 2015. Evolutionary plant radiations: where, when, why and how? New Phytol. 207(2):249-253.
Joshi C.P., Chiang V.L. 1998. Conserved sequence motifs in plant S-
adenosyl-L-methionine-dependent methyltransferases. Plant Mol.
Biol. 37(4):663-674.

Junier T., Zdobnov E.M. 2010. The Newick utilities: high-throughput phylogenetic tree processing in the UNIX shell. Bioinformatics 26(13):1669-1670.

Kalyaanamoorthy S., Minh B.Q., Wong T.K., Von Haeseler A., Jermiin L.S. 2017. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat. Methods 14(6):587-589.

Keller I., Wagner C.E., Greuter L., Mwaiko S., Selz O.M., Sivasundar A., Wittwer S., Seehausen O. 2013. Population genomic signatures of divergent adaptation, gene flow and hybrid speciation in the rapid radiation of Lake Victoria cichlid fishes. Mol. Ecol. 22(11):2848-2863.

Kircher M., Sawyer S., Meyer M. 2012. Double indexing overcomes in-accuracies in multiplex sequencing on the Illumina platform. Nucleic Acids Res. 40(1):e3.

Koch M.A., Kiefer C., Möbus J., Quandt D., Merklinger F., Harpke D., Benavides F.V. 2022. Range expansion and contraction of *Tillandsia landbeckii* lomas in the hyperarid Chilean Atacama Desert indicates ancient introgression and geneflow. Perspect. Plant Ecol. Evol. Syst. 56:125689.

Krueger F. 2019. Trim Galore: a wrapper tool around Cutadapt and FastQC to consistently apply quality and adapter trimming to FastQ files. Cambridge: Babraham Bioinformatics, Babraham Institute.

Kubatko L.S., Degnan J.H. 2007. Inconsistency of phylogenetic estimates from concatenated data under coalescence. Syst. Biol. 56(1):17-24.

Lagomarsino L.P., Condamine F.L., Antonelli A., Mulch A., Davis C.C. 2016. The abiotic and biotic drivers of rapid diversification in Andean bellflowers (Campanulaceae). New Phytol. 210(4): 1430-1442.

Lamichhaney S., Berglund J., Almén M.S., Maqbool K., Grabherr M., Martinez-Barrio A., Promerová M., Rubin C.-J., Wang C., Zamani N., Grant B.R., Grant P.R., Webster M.T., Andersson L. 2015. Evolution of Darwin's finches and their beaks revealed by genome sequencing. Nature. 518(7539):371-375.

Langmead B., Salzberg S.L. 2012. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9(4):357-359.

Leroy T., Louvet J.M., Lalanne C., Le Provost G., Labadie K., Aury J.M., Delzon S., Plomion C., Kremer A. 2020. Adaptive introgression as a driver of local adaptation to climate in European white oaks. New Phytol.. 226(4):1171-1182.

Léveillé-Bourret É., Starr J.R., Ford B.A., Moriarty Lemmon E., Lemmon A.R. 2018. Resolving rapid radiations within angiosperm families using anchored phylogenomics. Syst. Biol. 67(1): 94-112.

Li H. 2011. A statistical framework for SNP calling, mutation discovery, association mapping and population genetic parameter estimation from sequencing data. Bioinformatics 27(21): 2987-2993.

Linck E., Battey C.J. 2019. On the relative ease of speciation with periodic gene flow. bioRxiv. doi:10.1101/758664.

Lindenbaum P. 2015. JVarkit: Java-based utilities for bioinformatics [computer software]. doi:10.6084/m9 figshare.1425030.

Linder H.P. 2008. Plant species radiations: where, when, why? Philos. Trans. R. Soc. B Biol. Sci. 363(1506):3097-3105.

Liu L., Wu S., Yu L. 2015. Coalescent methods for estimating species trees from phylogenomic data. J. Syst. Evol. 53(5):380-390.

Loiseau O., Olivares I., Paris M., de La, Harpe M., Weigand A., Koubínová D., Rolland J., Bacon C.D., Balslev H., Borchsenius F., Cano A., Couvreur T.L.P., Delnatte C., Fardin F., Gayot M., Mejía F., Mota-Machado T., Perret M., Roncal J., Sanin M.J., Stauffer F., LEXER C., Kessler M., Salamin N. 2019. Targeted capture of hundreds of nuclear genes unravels phylogenetic relationships of the diverse Neotropical palm tribe Geonomateae. Front. Plant Sci. 10:864.

Luther H.E. 1985. Notes on hybrid tillandsias in Florida. Phytologia 53:175-176.

McKain M.R., Johnson M.G., Uribe-Convers S., Eaton D., Yang Y. 2018. Practical considerations for plant phylogenomics. Appl. Plant Sci. 6(3):e1038.

Malinsky M., Challis R.J., Tyers A.M., Schiffels S., Terai Y., Ngatunga B.P., Miska E.A., Durbin R., Genner M.J., Turner G.F. 2015. Genomic islands of speciation separate cichlid ecomorphs in an East African crater lake. Science 350(6267):1493-1498.

Malinsky M., Matschiner M., Svardal H. 2021. Dsuite—fast D-statistics and related admixture evidence from VCF files. Mol. Ecol. Resour. 21(2):584-595.

Malinsky M., Svardal H., Tyers A.M., Miska E.A., Genner M.J., Turner G.F., Durbin R. 2018. Whole-genome sequences of Malawi cichlids reveal multiple radiations interconnected by gene flow. Nat. Ecol. Evol. 2(12):1940-1955.

Mallet J., Besansky N., Hahn M.W. 2016. How reticulated are species? Bioessays 38(2):140-149.

Manichaikul A., Mychaleckyj J.C., Rich S.S., Daly K., Sale M., Chen W.-M. 2010. Robust relationship inference in genome-wide association studies. Bioinformatics 26(22):2867-2873.

Martin S.H., Davey J.W., Jiggins C.D. 2015. Evaluating the use of ABBA-BABA statistics to locate introgressed loci. Mol. Biol. Evol. 32(1):244-257.

Martin S.H., Jiggins C.D. 2017. Interpreting the genomic landscape of introgression. Curr. Opin. Genet. Dev. 47:69-74.

Martin S.H., Van Belleghem S.M. 2017. Exploring evolutionary relationships across the genome using topology weighting. Genetics 206(1):429-438.

Meier J.I., Stelkens R.B., Joyce D.A., Mwaiko S., Phiri N., Schliewen U.K., Selz O.M., Wagner C.E., Katongo C., Seehausen O. 2019. The coincidence of ecological opportunity with hybridization explains rapid adaptive radiation in Lake Mweru cichlid fishes. Nat. Commun. 10(1):5391.

Messerschmid T.F.E., Wehling J., Bobon N., Kahmen A., Klak C., Los J.A., Nelson D.B., dos Santos P., de Vos J.M., Kadereit G. 2021. Carbon isotope composition of plant photosynthetic tissues reflects a Crassulacean acid metabolism (CAM) continuum in the majority of CAM lineages. Perspect. Plant Ecol. Evol. Syst. 51:125619.

Meyer M., Kircher M. 2010. Illumina sequencing library preparation for highly multiplexed target capture and sequencing. Cold Spring Harb. Protoc. 2010(6):pdb.prot5448.

Mirarab S. 2019. Species tree estimation using ASTRAL: practical considerations. arXiv. doi:10.48550/arXiv:1904.03826.

Mitchell A., Mitter C., Regier J.C. 2000. More taxa or more characters revisited: combining data from nuclear protein-encoding genes for phylogenetic analyses of Noctuoidea (Insecta: Lepidoptera). Syst. Biol. 49(2):202-224.

Mondragon D., Calvo-Irabien L.M. 2006. Seed dispersal and germination of the epiphyte *Tillandsia brachycaulos* (Bromeliaceae) in a tropical dry forest, Mexico. Southwestern Nat. 51(4):462-470.

Morales-Briones D.F., Kadereit G., Tefarikis D.T., Moore M.J., Smith S.A., Brockington S.F., Timoneda A., Yim W.C., Cushman J.C., Yang Y. 2021. Disentangling sources of gene tree discordance in phylogenomic data sets: testing ancient hybridizations in Amaranthaceae s.l. Syst. Biol. 70(2):219-235.

Moran B.M., Payne C., Langdon Q., Powell D.L., Brandvain Y., Schumer M. 2021. The genomic consequences of hybridization. eLife 10:e69016.

Naciri Y., Linder H.P. 2020. The genetics of evolutionary radiations. Biol. Rev. 95(4):1055-1072.

Nguyen L.-T., Schmidt H.A., von Haeseler A., Minh B.Q. 2015. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32(1): 268-274.

Nosil P. 2008. Speciation with gene flow could be common. Mol. Ecol. 17(9):2103-2106.

Oliver J.C. 2013. Microevolutionary processes generate phylogenomic discordance at ancient divergences. Evolution 67(6):1823-1830.

Ortiz E.M. 2019. vcf2phylip v2.0: convert a VCF matrix into several matrix formats for phylogenetic analysis. Zenodo. doi:10.5281/zenodo.2540861.

Parins-Fukuchi C., Stull G.W., Smith S.A. 2021. Phylogenomic conflict coincides with rapid morphological innovation. Proc. Natl. Acad. Sci. USA 118(19):e2023058118.
Parvathi M.S., Nataraja K.N., Reddy Y.A., Naika M.B., Gowda M.V. 2019. Transcriptome analysis of finger millet (Eleusine coracana (L.) Gaertn.) reveals unique drought responsive genes. J. Genet. 98(2):1-12.

Paun O., Turner B., Trucchi E., Munzinger J., Chase M.W., Samuel R. 2016. Processes driving the adaptive radiation of a tropical tree (Diospyros, Ebenaceae) in New Caledonia, a biodiversity hotspot. Syst. Biol. 65(2):212-227.

Pease J.B., Haak D.C., Hahn M.W., Moyle L.C. 2016. Phylogenomics reveals three sources of adaptive variation during a rapid radiation. PLoS Biol. 14(2):e1002379.

Pérez-Escobar O.A., Chomicki G., Condamine F.L., Karremans A.P., Bogarín D., Matzke N.J., Silvestro D., Antonelli A. 2017. Recent origin and rapid speciation of Neotropical orchids in the world's richest plant biodiversity hotspot. New Phytol. 215(2):891-905.

Pham K.K., Hipp A.L., Manos P.S., Cronn R.C. 2017. A time and a place for everything: phylogenetic history and geography as joint predictors of oak plastome phylogeny. Genome 60(9):720-732.

Picard Toolkit. 2018. Broad Institute, GitHub Repository. Available from: http://broadinstitute.github.io/picard/

Pinzón J.P., Ramírez-Morillo I.M., Carnevali G., Barfuss M.H., Till W., Tun J., Ortiz-Díaz J.J. 2016. Phylogenetics and evolution of the *Tillandsia utriculata* complex (Bromeliaceae, Tillandsioideae) inferred from three plastid DNA markers and the ETS of the nuclear ribosomal DNA. Bot. J. Linn. Soc. 181(3):362-390.

Poplin R., Ruano-Rubio V., DePristo M.A., Fennell T.J., Carneiro M.O., der Auwera G.A.V., Kling D.E., Gauthier L.D., Levy-Moonshine A., Roazen D., Shakir K., Thibault J., Chandran S., Whelan C., Lek M., Gabriel S., Daly M.J., Neale B., MacArthur D.G., Banks E. 2018. Scaling accurate genetic variant discovery to tens of thousands of samples. BioRxiv. doi:10.1101/201178.

Quinlan A.R., Hall I.M. 2010. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26(6):841-842.

R Core Team. 2020. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing.

Renaud G., Stenzel U., Maricic T., Wiebe V., Kelso J. 2015. deML: robust demultiplexing of Illumina sequences using a likelihood-based approach. Bioinformatics 31(5):770-772.

Revell L.J. 2012. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3(2):217-223.

Richardson J.E., Pennington R.T., Pennington T.D., Hollingsworth P.M. 2001. Rapid diversification of a species-rich genus of Neotropical rain forest trees. Science 293(5538):2242-2245.

Rivera N.L. 2019. Exploring genomic relationships in *Tillandsia* subgenus *Tillandsia* [master's thesis]. University of Vienna. doi:10.25365/thesis.57906.

Rokas A., Williams B.L., King N., Carroll S.B. 2003. Genome-scale approaches to resolving incongruence in molecular phylogenies. Nature 425(6960):798-804.

Rundell R.J., Price T.D. 2009. Adaptive radiation, nonadaptive radiation, ecological speciation and nonecological speciation. Trends Ecol. Evol. 24(7):394-399.

Salichos L., Rokas A. 2013. Inferring ancient divergences requires genes with strong phylogenetic signals. Nature 497(7449):327-331.

Sankararaman S., Mallick S., Dannemann M., Prüfer K., Kelso J., Pääbo S., Patterson N., Reich D. 2014. The genomic landscape of Neanderthal ancestry in present-day humans. Nature 507(7492): 354-357.

Schluter D. 2000. The ecology of adaptive radiation. Oxford: OUP.

Seehausen O., Butlin R.K., Keller I., Wagner C.E., Boughman J.W., Hohenlohe P.A., Peichel C.L., Saetre G.-P., Bank C., Brännström Å. 2014. Genomics and the origin of species. Nat. Rev. Genet. 15(3): 176-192.

Sistla S., Rao D.N. 2004. S-Adenosyl-L-methionine-dependent restriction enzymes. Crit. Rev. Biochem. Mol. Biol. 39(1):1-19.

Slovak M., Melichárková A., Štubňová E.G., Kučera J., Mandáková T., Smyčka J., Lavergne S., Passalacqua N.P., Vďačný P., Paun O. 2023. Pervasive introgression during rapid diversification of the European mountain genus *Soldanella* (L.) (Primulaceae). Syst. Biol. 72(3):491-504.

Smith L.B., Downs R.J. 1977. Tillandsioideae (Bromeliaceae). Flora Neotrop. 14:798-800.

Smith S.A., Moore M.J., Brown J.W., Yang Y. 2015. Analysis of phylogenomic datasets reveals conflict, concordance, and gene duplications with examples from animals and plants. BMC Evol. Biol. 15(1): 150.

Soltis P.S., Folk R.A., Soltis D.E. 2019. Darwin review: angiosperm phylogeny and evolutionary radiations. Proc. R. Soc. B Biol. Sci. 286(1899):20190099.

Stankowski S., Chase M.A., Fuiten A.M., Rodrigues M.F., Ralph P.L., Streisfeld M.A. 2019. Widespread selection and gene flow shape the genomic landscape during a radiation of monkeyflowers. PLoS Biol. 17(7):e3000391.

Straub S.C.K., Moore M.J., Soltis P.S., Soltis D.E., Liston A., Livshultz T. 2014. Phylogenetic signal detection from an ancient rapid radiation: effects of noise reduction, long-branch attraction, and model selection in crown clade Apocynaceae. Mol. Phylogenet. Evol. 80:169-185.

Stroud J.T., Losos J.B. 2016. Ecological opportunity and adaptive radiation. Annu. Rev. Ecol. Evol. Syst. 47(1):507-532.

Suarez-Gonzalez A., Lawyer C., Cronk Q.C.B. 2018. Adaptive introgression: a plant perspective. Biol. Lett. 14(3):20170688.

Taylor S.A., Larson E.L. 2019. Insights from genomes into the evolutionary importance and prevalence of hybridization in nature. Nat. Ecol. Evol. 3(2):170-177.

Than C., Ruths D., Nakhleh L. 2008. PhyloNet: a software package for analyzing and reconstructing reticulate evolutionary relationships. BMC Bioinf. 9(1):1-16.

Till W. 2000. Tillandsioideae. In: Benzing D.H., editor. Bromeliaceae: profile of an adaptive radiation. New York: Cambridge University Press. p. 555-571.

Vera-Paz S.I., Granados Mendoza C., Díaz Contreras Díaz D.D., Jost M., Salazar G.A., Rossado A.J., Montes-Azcué C.A., Hernández-Gutiérrez R., Magallon S., Sánchez-González L.A., Gouda E.J., Cabrera L.I., Ramírez-Morillo I.M., Flores-Cruz M., Granados-Aguilar X., Martínez-García A.L., Hornung-Leoni C.T., Barfuss M.H.J., Wanke S. 2023. Plastome phylogenomics reveals an early Pliocene North- and Central America colonization by long-distance dispersal from South America of a highly diverse bromeliad lineage. Front. Plant Sci. 14:1205511.

Victoriano-Romero E., Valencia-Díaz S., Toledo-Hernández V. H., Flores-Palacio A. 2017. Dispersal limitation of *Tillandsia* species correlates with rain and host structure in a central Mexican tropical dry forest. PLoS One 12(2):e0171614.

Wang L.-G., Lam T.T.-Y., Xu S., Dai Z., Zhou L., Feng T., Guo P., Dunn C.W., Jones B.R., Bradley T., Zhu H., Guan Y., Jiang Y., Yu G. 2020. Treeio: an R package for phylogenetic tree input and output with richly annotated and associated data. Mol. Biol. Evol. 37(2): 599-603.

Wen D., Yu Y., Zhu J., Nakhleh L. 2018. Inferring phylogenetic networks using PhyloNet. Syst. Biol. 67(4):735-740.

Wendel J.F., Doyle J.J. 1998. Phylogenetic incongruence: window into genome history and molecular evolution. In: Soltis D.E., Soltis P.S., Doyle J.J., editors. Molecular systematics of plants II: DNA sequencing. Boston, MA: Springer US. p. 265-296.

Whitfield J.B., Lockhart P.J. 2007. Deciphering ancient rapid radiations. Trends Ecol. Evol. 22(5):258-265.

Winkler S. 1990. Zur Evolution der Gattung *Tillandsia* L. Bot. Jahrb. Syst. 112:43-77.

Wogan G.O.U., Yuan M.L., Mahler D.L., Wang I.J. 2023. Hybridization and transgressive evolution generate diversity in an adaptive radiation of Anolis lizards. Syst. Biol. 72(4): 874-884.

Wong E.L., Hiscock S.J., Filatov D.A. 2022. The role of interspecific hybridisation in adaptation and speciation: insights from studies in Senecio. Front. Plant Sci. 13:907363.

Wortley A.H., Rudall P.J., Harris D.J., Scotland R.W. 2005. How much data are needed to resolve a difficult phylogeny? Case study in Lamiales. Syst. Biol. 54(5):697-709.

Yardeni G., Viruel J., Paris M., Hess J., Groot Crego C., de La Harpe M., Rivera N., Barfuss M.H.J., Till W., Guzmán-Jacob V., Krömer T.,
Lexer C., Paun O., Leroy T. 2022. Taxon-specific or universal? Using target capture to study the evolutionary history of rapid radiations. Mol. Ecol. Resour. 22(3):927–945.

Yu G. 2020. Using ggtree to visualize data on tree-like structures. Curr. Protoc. Bioinformatics 69(1):e96.

Yu Y., Nakhleh L. 2015. A maximum pseudo-likelihood approach for phylogenetic networks. BMC Genomics 16(S10):S10.

Zeberg H., Pääbo S. 2020. The major genetic risk factor for severe COVID-19 is inherited from Neanderthals. Nature 587(7835):610–612.

Zhang C., Rabiee M., Sayyari E., Mirarab S. 2018. ASTRAL-III: polynomial time species tree reconstruction from partially resolved gene trees. BMC Bioinf. 19(S6):153.

Zhang D., Rheindt F.E., She H., Cheng Y., Song G., Jia C., Qu Y., Alström P., Lei F. 2021. Most genomic loci misrepresent the phylogeny of an avian radiation because of ancient gene flow. Syst. Biol. 70(5):961–975.

Zheng X., Levine D., Shen J., Gogarten S.M., Laurie C., Weir B.S. 2012. A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics 28(24):3326–3328.