Full plastomes of a total of 200 species were analyzed, of which 162 were newly sequenced and assembled for this study and 38, generated by Guisinger et al. (1 sp.; 2010) and Vera-Paz et al. (37 spp.; 2022), were downloaded from GenBank ( Supplementary Material 1 ). The sampled species were selected to construct a phylogenetic context with key nodes represented for 1) performing secondary calibrations, i.e., stem and crown nodes of Bromeliaceae and the crown node of Tillandsioideae, and 2) reconstructing the geographic origin and dispersal history of the expanded clade K. We aimed at achieving a dense geographic and taxonomic sampling of the expanded clade K. For this, we included 56 species previously recovered by Barfuss et al. (2016) and Granados Mendoza et al. (2017) within the expanded clade K, plus 37 additional Tillandsia subg. Tillandsia species mainly distributed in Central America and Mexico, which we hypothesized might belong in the expanded clade K given the geographic distribution pattern reported by Granados Mendoza et al. (2017). Sampling within T. subg. Tillandsia was complemented with three and five species previously recovered by Barfuss et al. (2016) within the T. secunda and T. paniculata clades, respectively, plus two additional species distributed in the Caribbean and South America and included here for the first time in a phylogenetic context, resulting in a total of 103 sampled species for T. subg. Tillandsia. Other subgenera and species complexes within Tillandsia, recognized by Barfuss et al. (2016), were sampled as follows: T. subg. Aerobia (8 spp.), T. subg. Anoplophytum s. str. (5 spp.), T. subg. Diaphoranthema (5 spp.), T. subg. Phytarrhiza s. str. (3 spp.), T. subg. Pseudovriesea (2 spp.), T. subg. Viridantha (4 spp.), T. australis complex (1 sp.), T. biflora complex (21 spp.), T. disticha complex (1 sp.), T. gardneri complex (3 spp.), T. purpurea complex (1 sp.), T. rauhii complex (2 spp.), and T. sphaerocephala complex (5 spp.), as well as the unclassified species T. albertiana and T. esseriana. Tribe Tillandsieae was further represented by species of the genera Barfussia (2 spp.), Guzmania (2 spp.), Lemeltonia (1 sp.), Pseudalcantarea (2 spp.), Racinaea (3 spp.), and Wallisia (2 spp.). Within tribe Vrieseeae, subtribe Cipuropsidinae was represented by the genera Goudaea (2 spp.), Lutheria (2 spp.), Mezobromelia (1 sp.), Werauhia (3 spp.), and Zizkaea (1 sp.), whereas for subtribe Vrieseinae the genera Alcantarea (2 spp.), Stigmatodon (2 spp.), and Vriesea (4 spp.) were sampled. We included three species of the genus Catopsis as representatives of tribe Catopsideae. Brocchinia micrantha (subfamily Brocchinioideae) and Typha latifolia (Typhaceae) were sampled to represent the crown and stem nodes of Bromeliaceae, respectively. Typha latifolia was used to root the phylogenetic tree following the Poales familial phylogenetic relationships from the Kew Tree of Life Explorer (https://treeoflife.kew.org/). Lineages of Tillandsioideae not available to us include tribe Glomeropitcairnieae, and genera Gregbrownia (tribe Tillandsieae), Josemania, Cipuropsis, Jagrantia, and Waltillia (tribe Vrieseeae; Barfuss et al., 2016; Leme et al., 2017). Classification within Tillandsioideae follows Barfuss et al. (2016) and Gouda et al. (continuously updated), whereas the subfamilial classification of Bromeliaceae follows Givnish et al. (2007).

New plastomes were assembled from raw data derived from a Hyb-Seq approach that applies a modified version of the universal probe kit for angiosperms of Buddenhagen et al. (2016), which includes additional Bromeliaceae target sequences. A previous publication (Vera-Paz et al., 2022) demonstrated the feasibility of assembling complete plastomes from data derived from this Hyb-Seq project, and herein-employed molecular methods were according to these authors. DNA isolation was performed at the “Laboratorio de Biología Molecular, Laboratorio Nacional de Biodiversidad (LANABIO)” of the Institute of Biology of the National Autonomous University of Mexico (IB-UNAM). The DNeasy Plant Pro Kit (Qiagen) was used to isolate genomic DNA from silica-gel-dried leaf tissue. A minimum DNA concentration of 1.0 µg/µl and a DNA purity 260/280 ratio ≥1 was targeted. DNA concentration and purity were quantified with a Qubit™ fluorometer v. 3.0 (Invitrogen™ Thermo Fisher Scientific) using the Qubit™ dsDNA BR Assay Kit (Invitrogen™ Thermo Fisher Scientific) and a NanoDrop 2000 spectrophotometer (Thermo Scientific™), respectively. Genomic DNA molecular weight was assessed in 1% agarose test gels, ran at 85 V and 500 mA for 50 min. Library preparation, enrichment, and sequencing were carried out by Daicel Arbor Biosciences (https://arborbiosci.com/). There, DNA concentration was additionally assessed by spectrofluorimetric assays with a PicoGreen assay kit (Thermo Fisher Scientific). A Qsonica Q800 ultrasonicator was used to fragment genomic DNA to a target insert length of 300–800 nt. A proprietary modification by Daicel Arbor Biosciences of the KAPA HyperPrep kit (Roche) protocol was applied for library preparation, using custom unique dual-index combinations for the samples. Spectrofluorimetric assays and quantitative PCR assays using a KAPA Library Quantification Kit (Roche) were carried out to quantify the indexed libraries. Capture was carried out in pools of 12 libraries per reaction following the myBaits v. 5 protocol (https://arborbiosci.com/mybaits-v5-chemistry/). After capture, enrichment reactions were amplified for 12 cycles and quantified with spectrofluorimetric and quantitative PCR assays. Based on the number of libraries in each capture, the captures were pooled in approximately equimolar ratios and combined with equimolar pools of non-captured libraries at a 70:30 ratio. Sequencing was performed on the Illumina NovaSeq 6000 platform on S4 PE150 lanes targeting 2 Gbp of sequencing effort per library.

In general, we followed the plastome assembly and annotation methods described in Vera-Paz et al. (2022). In short, raw read quality was visualized in FastQC v. 0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Then, the reads’ leading and trailing low-quality or N bases as well as adapters were trimmed with Trimmomatic v. 39 (Bolger et al., 2014). Success of the trimming process was verified in FastQC. Trimmed reads were used to perform de novo assemblies in GetOrganelle v. 1.5.0 (Jin et al., 2020), applying default settings for 150 bp long paired data. In cases where GetOrganelle assembled partial plastome sequences as separated scaffolds, the SPAdes v. 3.15.5 (Prjibelski et al., 2020) assembler was used in an effort to try and expand the scaffolds further. Using Tillandsia utriculata (ON398158; Vera-Paz et al., 2022) as reference, annotations were automatically transferred to the newly assembled plastomes in Geneious Prime^® 2021.2.2 (https://www.geneious.com/) and manually curated. Inverted repeat regions were identified with the repeat finder plugin of Geneious Prime^®. Newly sequenced and assembled full plastomes can be found in GenBank under the accession numbers OQ935587–OQ935748.

Full plastomes were aligned using the MAFFT v. 7.490 (Katoh and Standley, 2013) plugin of Geneious Prime^©. The resulting alignment was manually curated and inspected to locate inversions. Gblocks v. 0.91b (Castresana, 2000) was then used to identify ambiguously aligned regions using default settings, except for the allowed gap position parameter which was set to the option “all.” Then, the -E command of RAxML v. 8.2.10 (Stamatakis, 2014) was applied to exclude: 1) small (<200 pb long) inversions, 2) ambiguously aligned regions detected by Gblocks, and 3) one of the inverted repeat regions. Long inversions (>200 pb long) detected in genus Vriesea, three species from the Tillandsia biflora complex clade N (i.e., T. deppeana, T. heterophylla, and T. imperialis), and Werauhia gladioliflora, were replaced by “?” in the species where they were present and considered as missing data. This modified alignment was partitioned by protein-coding genes, tRNA genes, rRNA genes, introns, and intergenic spacers (IGS) and analyzed in PartitionFinder v. 2.1.1 (Lanfear et al., 2016) to select the best-fit partitioning scheme and corresponding models of molecular evolution using the Bayesian Information Criterion (BIC).

Estimation of phylogenetic relationships and divergence times was carried out simultaneously in BEAST2 v. 2.6.1 (Bouckaert et al., 2014). In BEAUti v. 2.5, we set the analysis priors for two data partitions retrieved as the best-fit partitioning scheme by PartitionFinder, applying the GTR and GTR+G models of molecular evolution, respectively, and empirical base frequencies for both site models. An uncorrelated relaxed clock model was applied, with rates sampled from a log-normal prior (Drummond et al., 2006) and a birth-death model tree prior. Based on the phylogenetic relationships recovered by Vera-Paz et al. (2022) and to reduce the computational cost, the topology was constrained for the major clades Bromeliaceae, Tillandsioideae, Catopsideae, Vrieseeae, and Tillandsieae, leaving unconstrained relationships among and within these clades. Given the lack of fossils that can be unambiguously assigned to Bromeliaceae and dated with confidence (e.g., Baresch et al., 2011), three secondary calibrations were applied using uniform prior distributions obtained from the relaxed (conservative fossil set) calibration performed by Ramírez-Barahona et al. (2020), as follows: 1) the root of the tree, i.e., stem node of Bromeliaceae (lower: 90.8908 Ma, upper: 123.4522 Ma), 2) the crown node of Bromeliaceae (lower: 15.3615 Ma, upper: 37.7002 Ma), and 3) the crown node of Tillandsioideae (lower: 7.3867, upper: 19.3439). Seven independent Markov Monte Carlo (MCMC) chains were run for 271–300 million generations each, sampling every 5,000 steps, until convergence was achieved, which was verified in Tracer v1.7.1 (Suchard et al., 2018). MCMC chain results were summarized with LogCombiner v2.6.1, setting the burn-in to 10% and sampling every 50,000 trees. The maximum clade credibility (MCC) tree was generated in TreeAnnotator v. 2.6.0 and visualized in FigTree v1.4.4. A second phylogenetic analysis without topological constraints was carried out under maximum likelihood in IQ-TREE (Nguyen et al., 2015) with an automatic selection of the best-fit substitution model for the individual partitions and assessing node support with 1,000 ultrafast bootstrap replicates.

The alignment of the complete plastomes resulted in a 180,766 bp matrix. A total of 48,040 bp corresponding to one inverted repeat region, small inversions, and ambiguously aligned regions were excluded resulting in a 132,726 bp long matrix that was used for phylogenetic inference. Nearly all phylogenetic relationships of the BEAST2 analysis received strong statistical support (PP ≥0.85), excepting as noted below ( Figures 2 – 5 ). Genera, subgenera, and species complexes from which we included more than one accession were recovered as monophyletic, except for Goudaea, Lutheria, Tillandsia, and the T. biflora complex. Tribe Catopsideae was recovered as the sister lineage of tribes Vrieseeae plus Tillandsieae. In subtribe Vrieseinae, Alcantarea is sister of a clade including Stigmatodon and Vriesea, whereas subtribe Cipuropsidinae is further divided in two lineages, one including Lutheria splendens as sister of Werauhia and the other including Goudaea ospinae as successive sister of Goudaea aff. chrysostachys, Zizkaea, Mezobromelia, and Lutheria bibeatricis ( Figure 2 ). Within tribe Tillandsieae, Guzmania is the sister lineage of a clade containing the remaining sampled species, which is divided into two main lineages. The first of these lineages (PP = 0.60) consists of a clade including the Tillandsia disticha complex and Pseudalcantarea as successive sisters (PP = 0.74) of T. subg. Pseudovriesea, Barfussia, Lemeltonia, and Wallisia plus Racinaea. The second lineage (PP = 0.34) is further divided into two clades, the first of them including a clade of T. purpurea complex and T. subg. Viridantha as the successive sister of 1) a clade including T. australis and T. sphaerocephala complexes, 2) the T. biflora complex clade N, 3) a clade including the T. rauhii complex and the T. biflora complex clade P, 4) the T. gardneri complex, 5) T. albertiana, 6) T. esseriana, 7) T. subg. Anoplophytum s. str., 8) T. subg. Diaphoranthema, and 9) a clade of T. subg. Aerobia and T. subg. Phytarrhiza s. str. ( Figure 3 ). The second clade corresponds to T. subg. Tillandsia, where all but a few recent and shortly spaced apart nodes within the later mentioned clade K.2 were strongly supported ( Figures 4 , 5 ). Within T. subg. Tillandsia, three main clades were recovered as successive sister lineages ( Figure 4 ). In the first of them, T. extensa is the sister species of a clade composed of T. marnier-lapostollei, T. spiraliflora, and T. propagulifera (corresponding to T. secunda clade). The second main clade, T. paniculata clade, includes two lineages, the first one including T. hildae as sister of T. funckiana–T. flexuosa, whereas in the second T. elongata is sister to T. malzinei–T. heliconioides. The third main clade is highly diverse, herein including 93 spp., equivalent to the expanded clade K, which consists of two main clades herein identified as K.1 (31 spp.) and K.2 (62 spp.) following the clade nomenclature proposed by Granados Mendoza et al. (2017). Given the substantial number of species integrating these two clades, internal phylogenetic relationships will not be described in detail, but they are shown in full in Figures 4 , 5 . When considering only strongly supported relationships (BS ≥85 and PP ≥0.85), the ML analysis resulted in a tree topology almost entirely congruent with the BEAST2 analysis, except for the alternative sister position of T. somnians and T. ovatispicata relative to the clade of T. huarazensis, T. roezlii, and T. reversa, as well as the position of Goudaea aff. chrysostachys as sister of Mezobromelia capituligera instead of sister to the clade of Mezobromelia capituligera, Lutheria bibeatricis, and Zizkaea tuerckheimii.

Divergence between Bromeliaceae and Typhaceae was estimated to have occurred 96.71 Mya (95% HPD = 90.89–114.14 Mya), whereas estimated crown ages of Bromeliaceae and Tillandsioideae were 28.69 Mya (95% HPD = 21.06–37.58 Mya) and 17.7 Mya (95% HPD = 13.72–19.34 Mya), respectively. Tribe Catopsideae, with a crown age of 5.38 Mya (95% HPD = 2.73–8.69 Mya), diverged from tribes Vrieseeae and Tillandsieae 17.7 Mya (95% HPD = 13.72–19.34 Mya), whereas divergence of the latter two tribes was 11.42 Mya (95% HPD = 7.66–15.01 Mya). The crown age of tribe Vrieseeae was 6.25 Mya (95% HPD = 3.81–8.98 Mya), whereas that of its subtribes Vrieseinae and Cipuropsidinae were 5.39 Mya (95% HPD = 3.16–7.92 Mya) and 5.01 Mya (95% HPD = 2.97–7.25 Mya), respectively ( Figure 2 ). The crown age of tribe Tillandsieae was 8.83 Mya (95% HPD = 5.94–12 Mya), and within it, Tillandsia subg. Tillandsia had stem and crown ages of 7.89 (95% HPD = 5.22–10.73 Mya) and 7.34 Mya (95% HPD = 4.8–9.98 Mya), respectively ( Figure 3 ). The stem and crown ages of the expanded clade K focal group were 6.94 (95% HPD = 4.49–9.42 Mya) and 4.86 Mya (95% HPD = 3.03–6.79 Mya), respectively, whereas the crown ages of its two subclades K.1 and K.2 were 2.82 (95% HPD = 1.72–4.1 Mya) and 3.9 Mya (95% HPD = 2.45–5.5 Mya), respectively ( Figures 4 , 5 ). The estimated stem and crown ages of other sampled genera and species complexes are shown in Supplementary Material 4 , and point-age estimations for all recovered nodes are shown in Figures 2 – 5 .

Log likelihoods and estimated j, d, and e parameter values for each tested biogeographic model for both analyses are shown in Table 1 . The model that best fitted the data in both cases is BAYAREALIKE+j, explaining ca. 99% and 72% of the total predictive power of all assessed models in the analyses with the complete taxon sampling and Tillandsia subg. Tillandsia, respectively. Among all analyzed areas, the combined area of the Boreal and South Brazilian dominions (B), South American transition zone (C), and Pacific dominion (D) was recovered as the most probable ancestral area for the nodes corresponding to subfamily Tillandsioideae (8%), core Tillandsioideae (45%), tribe Tillandsieae (59%), and a node including all tribe Tillandsieae representatives except Guzmania (59.20%). Then, the three nodes preceding the expanded clade K, including that of T. subg. Tillandsia (53%), were reconstructed as having experienced a range contraction to the Boreal and South Brazilian dominions (B) and South American transition zone (C), followed by a non-contiguous area shift to the Mexican transition zone + Mesoamerican dominion + Nearctic region individual area (F; 98.60%). Within T. subg. Tillandsia, the node subtending the T. secunda clade is similarly reconstructed in the combined area of the Boreal and South Brazilian dominions (B) and South American transition zone (C; 78.60%), whereas the T. paniculata clade expanded its ancestral area back to the Pacific dominion (D; 48.50%; Figure 6 , see Supplementary Material 3C for the second most probable ancestral area).

Key nodes relevant for the biogeographic discussion of the expanded clade K are indicated in Figure 7 with the numbers 1–20 (see Supplementary Material 3D for the second most probable ancestral area). The most recent common ancestor of this focal group is estimated to have occupied the combined area of the Pacific Lowlands and Balsas Basin provinces + El Cabo de Baja California district (G), Mexican transition zone (H), and the Mosquito, Veracruzan, and Yucatán Peninsula provinces (I; 31%). From this ancestral area, four lineages within clade K.1 were reconstructed to have reduced their range to either the Mexican transition zone (H) and the Mosquito, Veracruzan, and Yucatán Peninsula provinces (I; node 1 and 5; 36% and 89%) or to the Pacific Lowlands and Balsas Basin provinces + El Cabo de Baja California district (G) and the Mexican transition zone (H; nodes 4 and 6; 50% and 38%), whereas one lineage experienced a range expansion to the Nearctic region (J; node 3; 53%). Further colonization of the Nearctic region (J) occurred via additional area range expansions from the Mexican transition zone (H) and the Mosquito, Veracruzan, and Yucatán Peninsula provinces (I; node 2, 47%), as well as the Pacific Lowlands and Balsas Basin provinces + El Cabo de Baja California district (G) and the Mexican transition zone (H; nodes 7, 12, 28% and 62%). Additional dispersals to the Mosquito, Veracruzan, and Yucatán Peninsula provinces (I) involved range expansions from the combined area of the Pacific Lowlands and Balsas Basin provinces + El Cabo de Baja California district (G), the Mexican transition zone (H), and the Nearctic region (J; node 8, 54%) or from the Pacific Lowlands and Balsas Basin provinces + El Cabo de Baja California district (G) and the Mexican transition zone (H), however, with a simultaneous area expansion to the Nearctic region (J, node 11, 54%). Successive area contraction and expansion were modeled at nodes 9 (47.9%) and 10 (80%) in combined areas involving the Pacific Lowlands and Balsas Basin provinces + El Cabo de Baja California district (G), the Mexican transition zone (H), and the Nearctic region (J). Within clade K.1, a single area range expansion back to the Pacific, Boreal and South Brazilian dominions + South American transition zone (Y) is reconstructed in a terminal branch for Tillandsia punctulata.

Ancestral area reconstruction within clade K.2 is only described for nodes that received strong statistical support (PP ≥0.85) in the phylogenetic analysis. Within this clade, several range contractions were reconstructed to either the Pacific Lowlands and Balsas Basin provinces + El Cabo de Baja California district (G) and the Mexican transition zone (H; nodes 14, 15, and 19; 74%, 72.50%, and 83.88%); the Mexican transition zone (H; node 13, 80.20%); or the Mosquito, Veracruzan, and Yucatán Peninsula provinces (I; node 18, 99%). Expansion to the Nearctic region (J) in clade K.2 occurred at nodes 17 (80%) and 20 (79%) from the ancestral area of the expanded clade K and from the Pacific Lowlands and Balsas Basin provinces + El Cabo de Baja California district (G) and the Mexican transition zone (H; node 16, 57%), however, with a simultaneous area range expansion to the Mosquito, Veracruzan, and Yucatán Peninsula provinces (I) and back to the Pacific, Boreal, and South Brazilian dominions + South American transition zone (Y). Eight additional area range expansions back to the Pacific, Boreal, and South Brazilian dominions + South American transition zone (Y) were estimated at terminal branches of clade K.2 corresponding to Tillandsia filifolia, T. leiboldiana, T. pseudobaileyi, T. ionantha, T. variabilis, T. schiedeana, T. brachycaulos, and T. tricolor. Colonization of the Antillean subregion + SE USA (E) occurred exclusively within clade K.2 in six independent events from either the combined area of the Pacific Lowlands and Balsas Basin provinces + El Cabo de Baja California district (G), Mexican transition zone (H), and the Mosquito, Veracruzan, and Yucatán Peninsula provinces (I, T. schiedeana, T. variabilis, and T. streptophylla), from the latter combined area plus the Nearctic region (J; T. setacea) or from an unknown area due to low statistical support of the phylogenetic relationships (T. bartramii and T. capitata).

Our study contributes 162 newly sequenced and assembled complete plastomes to Bromeliaceae genomic resources, more densely sampled within the expanded clade K of Tillandsia subg. Tillandsia but also representing many other Tillandsia subgenera and species complexes, as well as other genera of Tillandsioideae. Additionally, the complete plastome of Brocchinia micrantha (subfamily Brocchinioideae) was generated here as part of our outgroup sampling. Our analyses of 200 complete plastomes derived in an overall strongly supported phylogenetic hypothesis, with only 13 out of the 198 nodes (excluding the root) with PP <0.85. At tribe and subtribe levels, our results agree with previously proposed plastid phylogenetic contexts (Barfuss et al., 2016; Machado et al., 2020; Loiseau et al., 2021; Vera-Paz et al., 2022). As in the plastid-based tree of Loiseau et al. (2021), we here recovered the genus Goudaea (subtribe Cipuropsidinae) as non-monophyletic (see also Supplementary Material 5 ), albeit herein with strong statistical support. However, the nuclear-based inference of Loiseau et al. (2021) supports the monophyly of this genus. Interestingly, those authors propose a potential hybrid origin for Goudaea chrysostachys, which could potentially explain the discordance between nuclear and plastid data sources across studies. The genus Lutheria is also non-monophyletic in our phylogenetic context with strong statistical support. Remarkably, the alternative phylogenetic positions shown by the herein analyzed Lutheria species (L. splendens and L. bibeatricis) follow the same alternative positions proposed by Loiseau et al. (2021) for this genus upon comparing their plastid and nuclear trees. Our study is the first testing the phylogenetic position of L. bibeatricis; therefore, its affiliation to the genus Lutheria and the reasons (e.g., hybridization) for its unexpected phylogenetic position should be further investigated. Other non-monophyletic groups recovered in our analysis with strong statistical support are the genus Tillandsia and the T. biflora complex. The non-monophyly of Tillandsia results from the sister relationship of T. disticha and Pseudalcantarea, since the obtained sister relationship of T. subg. Pseudovriesea to the clade including Barfussia, Racinaea, Wallisia, and Lemeltonia received low statistical support. A comparatively long branch subtends T. disticha in a maximum likelihood tree ( Supplementary Material 5 ), a result shared with previous studies (e.g., Barfuss et al., 2016), suggesting that the phylogenetic position of this species could also be influenced by long branch attraction, leaving open the question if Tillandsia, as circumscribed by Barfuss et al. (2016), is or is not monophyletic. Additional statistical support is required to resolve recalcitrant nodes along the backbone of Core Tillandsieae (genus Tillandsia s.l., Smith and Downs, 1977) to validate or refute the segregation of various genera from Tillandsia s.l. (Grant and Zijlstra, 1998; Smith and Till, 1998; Barfuss et al., 2016). Since its recognition in Barfuss et al. (2016), the T. biflora complex (136 spp.) was defined as a polymorphic assemblage of species with contrasting geographic distributions, highlighting the presence of Caribbean and Mexican species as part of this complex. Although highly supported in the Bayesian analysis of Barfuss et al. (2016), the sister relationship of the T. biflora clades N and P received low statistical support in their maximum likelihood and parsimony analyses, which is in favor of the non-monophyly of this group. Also, Tillandsia biflora complex clades N and P show distinct biogeographical patterns, discussed below. The herein obtained phylogenetic relationships could help define which morphological attributes characterize these two independent lineages.

Three strongly supported lineages in successive sister relationship conform Tillandsia subg. Tillandsia, namely, the T. secunda clade, the T. paniculata clade, and the expanded clade K. These results are consistent with Barfuss et al. (2016), except that these authors recovered the T. secunda clade as the sister lineage of the expanded clade K instead of the T. paniculata clade as in our phylogeny. While we obtained strong statistical support for these relationships, Barfuss et al. (2016) obtained low statistical support in the ML and parsimony analyses, albeit a PP = 0.85 in the Bayesian analysis. The expanded clade K and its internal clades K.1 an K.2 have consistently been recovered as strongly supported monophyletic groups in the previous (Chew et al., 2010; Sidoti, 2015; Barfuss et al., 2016; Pinzón et al., 2016; Montes-Azcué, 2020; Vera-Paz et al., 2022; Yardeni et al., 2022) and present studies. Furthermore, species relationships and composition of clades K.1 and K.2 are congruent across previous studies (Sidoti, 2015; Barfuss et al., 2016; Pinzón et al., 2016; Montes-Azcué, 2020; Vera-Paz et al., 2022; Yardeni et al., 2022) and our study. While all internal relationships of the less diverse clade K.1 are strongly supported, several relationships within clade K.2, which contains the type species of the genus Tillandsia, T. utriculata, are subtended by short branches and lack strong statistical support. Resolution of these low supported phylogenetic relationships would require the use of more variable sequence data, such as low- or single copy nuclear loci. Although the matK-trnK region analysis of Granados Mendoza et al. (2017) included a dense taxon sampling of the expanded clade K, internal resolution of clades K.1 and K.2 was extremely limited. Other studies focused on particular species complexes or groups, such as the T. ionantha complex (Ancona et al., 2022), pseudobulbous species (Chew et al., 2010), the T. erubescens complex (Granados Mendoza, 2008; Martínez-García et al., 2022), the T. utriculata complex (2016), and the T. fasciculata complex (Sidoti, 2015), used combinations of Sanger-sequenced DNA markers and/or morphology. When considering strongly supported phylogenetic relationships, our study is in overall agreement with plastid-based phylogenetic relationships recovered in those previous studies. However, there are some disagreements with the nuclear ribosomal phylogenetic context of Chew et al. (2010), who recovered a pair of species from clade K.2 (T. pseudobaileyi) and K.1 (T. achyrostachys) as sister to one another, with strong support.

Many of the phylogenetic relationships within clades K.1 and K.2 are assessed here for the first time. Phylogenetic positions of species of proposed hybrid origin, such as Tillandsia maya (T. balbisiana × T. brachycaulos, Ramírez and Carnevali, 2003), T. may-patii (T. brachycaulos × T. dasyliriifolia, Ramírez and Carnevali, 2003), and T. rothii (T. jaliscomonticola × T. rolandgosselinii, Rauh, 1976), were consistently recovered as sister lineages of one of their proposed parental lineages, in agreement with the proposed hybrid origin of these species. Future studies using nuclear sequence data and phylogenetic network analyses could explicitly test these hybrid origin hypotheses and explore the prevalence of hybridization across the expanded clade K and its influence in the diversification of the group.

We present a time-calibrated phylogenetic framework for the expanded clade K based on a comprehensive amount of plastome sequence data. To include key nodes for secondary calibrations, a representative sampling (16 out of 22 genera) of other Tillandsioideae lineages was incorporated, allowing to date the age of the reconstructed ancestors of the expanded clade K all the way back to the stem node of Bromeliaceae. The following discussion will focus on the nodes preceding this focal group.

Regardless of the differences in phylogenetic dating method, taxonomic sampling, and calibration sources (fossils vs. secondary calibrations), most of our crown age estimates fall within the age ranges estimated in previous studies, including the crown ages of Bromeliaceae (28.69 vs. 96–19 Mya; Givnish et al., 2004; Janssen and Bremer, 2004; Givnish et al., 2007; Givnish et al., 2011; Bouchenak-Khelladi et al., 2014; Givnish et al., 2014; Zhou et al., 2018; Ramírez-Barahona et al., 2020), Core Tillandsioideae (11.42 vs. 12.9–8.7 Mya; Givnish et al., 2011; Givnish et al., 2014; Kessous et al., 2020; Loiseau et al., 2021; Möbus et al., 2021), and tribe Tillandsieae (8.8 vs. 8.8–6.5 Mya; Kessous et al., 2020; Loiseau et al., 2021; Möbus et al., 2021). In contrast, crown ages of subfamily Tillandsioideae (17.7 vs. 15.2–13.3 Mya; Givnish et al., 2011; Givnish et al., 2014; Ramírez-Barahona et al., 2020; Möbus et al., 2021) and Tillandsia subg. Tillandsia (7.34 vs. 6.6 Mya; Möbus et al., 2021) are estimated to be older than in previous studies. In both cases, differences in the employed taxon samplings and the use herein of a considerably greater amount of sequence data could explain these age differences. Our sampling, for instance, does not include representatives from tribe Glomeropitcairnieae, sampled by Givnish et al. (2011; 2014) and Ramírez-Barahona et al. (2020), which could have influenced our crown age estimates of subfamily Tillandsioideae. In the case of T. subg. Tillandsia, Möbus et al. (2021) only sampled two species, T. utriculata from the expanded clade K and T. malzinei from the T. paniculata clade, therefore representing the node succeeding T. subg. Tillandsia from our phylogeny for which we estimated an age of 6.94 Mya, which is closer to their age estimate for T. subg. Tillandsia. Our taxon sampling design, in combination with a comprehensive amount of sequence data, allowed us to calibrate for the first time several nodes, not only within our focal group, the expanded clade K, but also in other Tillandsioideae lineages. We expect that this dated phylogenetic framework will facilitate future macroevolutionary studies and provide reference age estimates for performing secondary calibrations for other Tillandsioideae lineages.