Yardeni et al. 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 The recent rapid radiation of Tillandsia subgenus Tillandsia (Bromeliaceae) provides an attractive system to study the drivers and limits 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 MY while expanding through one major event from South into Central America within the last 5 MY. However, disagreements between phylogenies and lack of resolution at shallow nodes suggested that hybridization occurred throughout the radiation, together with frequent incomplete lineage sorting and/or considerable gene family evolution. We used whole-genome resequencing data and a newly available reference genome to explore the evolutionary history of 34 representative ingroup species employing both a tree-based and a network approach. 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, while topology weighting uncovered high heterogeneity of the phylogenetic signal along the genome. High rates of hybridization within and among clades suggest that, in contrast to previous hypotheses, the expansion of subgenus Tillandsia into Central America proceeded in several dispersal events, punctuated by episodes of diversification and gene flow. Network analysis revealed reticulation as a prominent propeller during radiation and establishment in 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. Neotropical diversity, species network 34 35 36 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Rapid evolutionary radiations are characterised 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 (Hughes et al. 2015; Linder 2008; Naciri and Linder 2020; Rundell and Price 2009; Stroud and Losos 2016). Rapid radiations provide ample opportunities to investigate the mechanisms driving diversification. While not all rapid radiations are adaptive, i.e., not necessarily associated with an increase in ecological occupancy (Givnish 2015; Schluter 2000; 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 particular, plant radiations remain to date understudied and enigmatic, especially when compared to animal systems (but see e.g., Drummond et al. 2012; Lagomarsino et al. 2016; Linder 2008; Paun et al. 2016; Pérez-Escobar et al. 2017;

as low sequence divergence and incomplete reproductive barriers increase the probability of short internal branches and impede phylogenomic resolution (Giarla and Esselstyn 2015; Straub et al. 2014; Whitfield and Lockhart 2007). Frequent shifts at key traits can be distinctly associated with phylogenetic conflict, stemming from population-level processes like changes in population size, incomplete lineage sorting (ILS) and reticulation (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 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83

research indicates that episodes of phylogenomic conflict represent typical signatures of microevolutionary processes during rapid species diversification (Filiault et al. 2018; ParinsFukuchi et al. 2021). Hence, gene-tree conflict and non-bifurcating relationships should be viewed and studied as defining features of diversification processes during rapid radiations (Malinsky et al. 2018; Slovak et al. 2023; Wogan et al. 2023).

Abbott et al. 2013; Seehausen 2004) with gene flow regarded as an important process during speciation (Barth et al. 2020; Linck and Battey 2019; Mallet et al. 2016). Gene flow between recently diverged species is thought to increase their genetic variability and enrich the genomic substrate for natural selection to act on. Interspecific hybridization can prompt ecological diversification through heterosis and adaptive introgression, a process in which advantageous alleles are incorporated in new genomic backgrounds from one gene pool to another. Furthermore, introgression can introduce advantageous alleles upon which selection has already acted, thus catalysing evolutionary radiations (Abbott et al. 2013; Edelman et al. 2019; Harrison and Larson 2014; Meier et al. 2019; Suarez-Gonzalez et al. 2018).

provides a particularly relevant study system to investigate the drivers of radiations. Subgenus Tillandsia is the most diverse of seven subgenera within genus Tillandsia (family Bromeliaceae). It consists of more than 250 species of predominantly epiphytic plants, distributed from 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 disparate habitats, from tropical rainforests to deserts, and from low- to highlands (Benzing 2000; Givnish et al. 2011). Adaptations such as adventitious 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 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 recovered largely inconsistent phylogenetic relationships within the subgenus (Barfuss et al. 2005; 2016; Chew et al. 2010; Granados Mendoza et al. 2017; Pinzón et al. 2016). Using plastome phylogenomics, Vera-Paz et al. (2023) confirmed the monophyly of the subgenus and the presence of three main clades within the subgenus, including a monophyletic, radiated ‘clade K’ with members distinctly distributed in North and Central America (Barfuss et al. 2016; Granados Mendoza et al. 2017). They further identified clades K.1 and K.2 as subclades of clade K which we have further renamed to reflect phylogenetic groups (see Methods; Supporting Table S1). Several studies (Givnish et al. 2014; Till 2000; Vera-Paz et al. 2023; Winkler, 1990) discussed that the ancestor of ‘clade K’ was most likely South American and expanded to Central America approximately 4.9 Mya, and subsequently dispersed multiple times to North America. However, the relationships within the clade remained elusive and disagreements between gene trees are common, calling to consider the role of different evolutionary processes in the Tillandsia radiation. As a result, hybridization was put forward as especially prevalent in the subgenus (Till 2000; Vera-Paz et al. 2023).

study aims to shed light on the evolutionary history of the radiation of Tillandsia subgenus Tillandsia, focusing on Central American groups within the radiated ‘clade K’ and representatives of South American species. Specifically, our objectives are to infer relationships between representative species and investigate the signals of phylogenetic

Pervasive hybridization in radiated Tillandsia 107 108 109 110 conflict, with a ultimate aim to characterise and quantify hybridization and gene-flow. Additionally, we highlight patterns of gene-flow which may have contributed to adaptive trait shift in this Neotropical rapid radiation. Our detailed analyses provide evidence of a deeply reticulated history within the subgenus and of at least two colonisation events from South into Central America, surprisingly refuting the monophyly of the so-called ‘clade K’. 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133

Materials and Methods

Material from field collections was cut lengthwise and immediately dried in powdered silica gel. Samples from the botanical collection of the University of Vienna were extracted fresh (Supporting Tables S1, S3). 69 accessions correspond to 34 species of Tillandsia subg. Tillandsia, while the remaining two belong to subg. Allardtia sensu Smith & Downs (1977) and were used as outgroups. Of the 34 ingroup species, 26 represent the Central-North American ‘clade K’ radiation and eight species represent South American Tillandsia (henceforth: SA clades; Supporting Fig. 1, Supporting Table 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’ (Supporting Table S2). Species in ‘clade K’ are mostly epiphytic and ornithophilous, whereas clade SA includes epiphytic, saxicolous and tank forming xerophitic species, predominantly endemic to Peru (with the exception of Tillandsia adpressiflora; Supporting Tables S1-2). The T. fasciculata clade (previously K.2.3) includes species largely exhibiting CAM photosynthetic syndrome with typical waterabsorbing trichomes. In addition, species within this clade present varying distributions, from endemic to widespread, some extending into the Antilles. The T. guatemalensis clade (previously K.2.2) is a small group of widespread and predominantly-C3 plants, exhibiting smooth leaves and a prominent tank habit. The T. punctulata (K.1) and T. utriculata (K.2.1) clades include widespread or endemic species with intermediate CAM-C3 syndromes or various CAM-like photosynthetic syndromes, respectively.

Pervasive hybridization in radiated Tillandsia 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154

a modified CTAB protocol (Doyle and Doyle, 1987) or the QIAGEN DNeasy® Plant Mini Kit (Qiagen, USA). They have been further purified using the Nucleospin® gDNA cleanup 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). Illumina 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 either double indexed with a set of 60 dual-index primers, as recommended by Kircher et al. (2012) and described in Loiseau et al. (2019), or with Illumina TrueSeq PCR-free dual indexes. Finally, the libraries underwent size selection steps using AMPure Beads (Agencourt). Libraries were then pooled and sequenced at Vienna BioCenter Core Facilities on Illumina HiSeqV4 PE125 or on NovaSeqS1 PE150. Protocol details for all accessions are listed in Supporting Table S3.

genome v.1.0 (GenBank Assembly accession GCA_029168755.1, Groot Crego et al., 2023) using bowtie2 v.2.3.5.1 (Langmead and Salzberg 2012) with the --very-sensitive-local option to increase accuracy. Low-quality mapped reads were removed, alignments were sorted by position using samtools v.1.15.1 (Li et al. 2009), and PCR duplicates were marked using MarkDuplicates from PicardTools v.2.25.2 (Picard Toolkit 2018). Relatedness between accessions was inferred with kinship coefficients in KING v.2.2.6 (Manichaikul et al. 2010), both for confirming technical replicates and to infer unexpected relatedness. Samples with an average coverage below 2x or a first degree relationship (consistent with full-sib or direct parent-offspring relationships; i.e., kinship > 0.177) were removed from subsequent analysis (Supporting Table S3).

Next, variants were called using a dedicated pipeline available at https://github.com/giyany/TillandsiaPhylo. In brief, we used

HaplotypeCaller v.4.1.9.0 followed by joint calling with GenotypeGVCFs (Poplin et al. 2018). The resulting variant call file (VCF) was filtered to exclude indels and any SNPs 3bp or closer from indels using bcftools v.1.15 (Li 2011) and GATK SelectVariants (DePristo et al. 2011). Regions annotated as transportable elements were excluded using bedtools intersect (Quinlan and Hall 2010). We used the transition/transversion ratio as a guideline to define the filtering parameter values leading to a set of SNPs exhibiting the highest ratio, calculated with SnpSift

(Cingolani et al. 2012). The following filtering parameters were finally used: MQ < 15, DP < 4, QD < 4, FS > 40, SOR > 3, MAF < 0.045 and missing rate > 0.2. Summary statistics were generated using bcftools stats (Li 2011).

the relationships within and between clades and to assess deviations from tree structure (Durand et al. 2011; Green et al. 2010; Martin et al. 2015). The D-statistic we 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 (Malinsky et al. 2021). Tillandsia complanata was used as the outgroup for all comparisons. We first examined the consistency of assigning individuals to a clade. Following Malinsky et al. (2018), we tested whether individuals assigned to the same clade always share more derived alleles with each other than with individuals from another clade. For this comparison and all related D-statistics hereafter, 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244

significance was assessed using a jackknife procedure with 200kb window size, and familywise error rate (FWER) was calculated and corrected for multiple comparisons following the

of clade monophyly within our set of window-based trees using the check_monophyly function in the ETE Toolkit v.3 (Huerta-Cepas et al. 2016). A cloudogram of ML trees was visualised using ggdensitree from the R package ggtree v.3.4.0 (Yu 2020). Trees were modified using R packages phytool v.1.0-3 (Revell 2012) and treeio v.1.20.0 (Wang et al. 2020), with branch lengths forced into ultrametric structure for visualisation purposes only. We next used Dmin, another variation of the D-statistic (Malinsky et al. 2018), to assess if allele-sharing among clades is consistent with a tree-like structure. Dmin measures the minimal amount 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, even in the presence of incomplete lineage sorting. In addition, a principal components analysis (PCA) was performed with SNPRelate v.1.20.1 (Zheng et al. 2012) on a set of biallelic SNPs, filtered to allow a maximum of 10% missing data and pruned to be 10kb or more apart.

used topology weighting by iterative sampling of subtrees with Twisst v.0.2 (Martin and Van Belleghem 2017). This method considers all the possible topologies relating several taxonomic groups and quantifies the contribution of each individual taxon topology to the full tree, enabling to locate genomic regions that are associated with certain topologies. We focused on the three largest clades and reduced the number of possible topologies by 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 including only clades T. extensa, T. utriculata and T. fasciculata, then chose windows of 50 SNPs following the strategy of Martin and Van Belleghem (2017). Subtrees were produced by partitioning the VCF with the biostar497922 script from Jvarkit (Lindenbaum 2015) and ML trees were inferred for each window as described above. A summary of all topologies and visualisation of topologies along genomic coordinates were produced using the plot_twisst.R script, modified to avoid averaging the signal over regions with high rates of missing data (Martin and Van Belleghem 2017).

Characterising Hybridization and Introgression

To quantify the rates of hybridization between all species in our dataset, we obtained a genome-wide signal of introgression using the original implementation of the D-statistic and estimates of admixture fraction f (henceforth: ƒ4-ratio) in Dsuite v.0.5r45 (Malinsky et al. 2021). We additionally computed D-statistic for each reference chromosome, obtaining pvalues using a jackknife procedure with windows of 150 SNPs. 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.

Since groups that are involved in hybridization may share branches on a phylogenetic tree, a single hybridization event can present multiple correlated instances of significantly elevated D-statistic. This is especially expected when introgression involves ancestral lineages, hence affecting internal branches of a phylogenetic tree. We used Dsuite to calculate the f-branch metric (henceforth: fb(C)), an estimator developed to create a summary of gene flow events with minimized correlation by utilising multiple ƒ4-ratio calculations 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289

(Malinsky et al. 2021). The fb(C) results invited several hypotheses regarding hybridization events, which we examined by inferring a species network under a maximum pseudolikelihood approach using PhyloNet v.3.8.0 (Cao et al. 2019; Than et al. 2008; Y. Yu and Nakhleh 2015), allowing different numbers of reticulation events. We reduced our sampling to one outgroup and 18 ingroup taxa to include representatives from each highly-supported clade and inferred ML trees on non-overlapping windows as previously described.

To investigate the signature of introgression on specific loci and to locate highly admixed loci in regions of interest (see Results), we used the Dinvestigate function implemented in Dsuite which calculates statistics suitable for analysis in sliding genomic windows. Our analysis was performed on windows of 50 SNPs with a step size of ten. The Dstatistic itself shows large variance when applied to genomic windows (Martin et al. 2015), hence we used fdM, a statistic designed to investigate introgression in small genomic windows. fdM accounts for allele sharing across all possible taxon arrangements in a species trio (Malinsky et al. 2015; Martin et al. 2015). After finding a significantly elevated D-statistic on chromosome 18 involving T. achyrostachys (from T. punctulata clade) and most species in the CAM clade T. fasciculata (see Results), we investigated patterns of hybridization between T. punctulata, T. butzii and T. achyrostachys (P1, P2 and P3). 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). Species in the T. punctulata clade putatively possess an intermediate C3-CAM photosynthetic syndrom, as they tend to express intermediary values of stable carbon isotope ratios (δ13C: see Supporting Table S2, but see also Messerschmid et al. 2021) and exhibit small to intermediate tank structures. T. achyrostachys was documented to express the strongest CAM phenotype and 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 δ13C values (−14.7; Crayn et al. 2015) within our sampling. We finally identified regions exhibiting D-statistic values exceeding the 95% quantile of the distribution to gain insight into the possible functional basis of genes that were involved in ancient hybridization.

shallow nodes (Fig. 1; Supporting Fig. S3). The relationships between the South American species and the T. utriculata clade were characterised 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 clades only 587 (15.5%) of all gene windows supported the main topology (Fig. 1; Supporting Fig. S3). Similar levels of discordance characterised many of the internal nodes within the T. fasciculata clade and high levels of discordance were also found within the intermediate CAM clade T. punctulata (Fig. 1; Supporting Fig. S3).

reference chromosome retrieved many different topologies: solely considering relationships between main clades we recognized ten different topologies among the 25 different trees (see Supporting file 1). The placement of the T. mima clade showed the greatest incongruence between trees, as well as the placement of T. adpressiflora. Within-clade relationships, like the whole-genome topology, were recovered from three chromosome-by-chromosome trees. The South American species formed a monophyly in six chromosome trees and ‘clade K’ was 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 recovered as monophyletic in two. While each chromosome tree contains high amounts of gene tree discordance, the abundance of different topologies imply that several evolutionary histories can be traced along the genome of Tillandsia, as disparate genomic processes muddy the relationship between gene trees and the true species tree.

Lack of Monophyly and Deviations from Tree-like Structure The multitude of tree topologies along the genome led us to hypothesise that a single bifurcating tree misrepresents the true relationships between Central American Tillandsia (i.e., previous ‘clade K’). To investigate potential deviations from a tree-like structure, we analysed the patterns of allele sharing between species and between clades. Assuming no inter-specific 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 calculated allele sharing between all species to test the robustness of the assignment of species to clades, utilising allele frequency estimates as performed with Dstatistic. We found high support for clade monophyly, as species assigned to the same clade always shared more alleles with other species within the clade than with species from other clades. Except for the T. punctulata clade, consistency of clade monophyly was further supported in analyses performed in genomic windows: monophyly was confirmed in 92.0%, 79.4%, 67.0% and 39.7% of the windows for clades T. utriculata, T. guatemalensis, T. fasciculata and T. punctulata, respectively.

We next characterised allele sharing between clades by computing the Dmin statistic for a total of 7,141 trios. We found widespread deviations from a tree-like structure, driven mostly by allele sharing between Central American taxa (previous ‘clade K’): Dmin values 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803

Mallet J., Besansky N., Hahn M.W. 2016. How reticulated are species? BioEssays. 38:140–

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: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: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: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: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:pdb.prot5448. Mirarab S. 2019. Species tree estimation using ASTRAL: practical considerations. ArXiv:1904.03826. 51:462-470.

Mondragon D., Calvo-Irabien, L. M. 2006. Seed dispersal and germination of the epiphyte Tillandsia brachycaulos (Bromeliaceae) in a tropical dry forest, Mexico. Southwest. Nat. 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: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: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:268–274. 824

Nosil P. 2008. Speciation with gene flow could be common. Mol Ecol. 17:2103-2106. 821 822 823 825 826 827 828 829 830 831 832 833 834 835 836

Oliver J.C. 2013. Microevolutionary processes generate phylogenomic discordance at ancient divergences. Evolution. 67:1823–1830.

Ortiz, E.M. 2019. vcf2phylip v2.0: convert a VCF matrix into several matrix formats for phylogenetic analysis. DOI:10.5281/zenodo.2540861 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:720-732.

Parins-Fukuchi C., Stull G.W., Smith S.A. 2021. Phylogenomic conflict coincides with rapid morphological innovation. Proc. Natl. Acad. Sci. 118: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:1–12.

Pervasive hybridization in radiated Tillandsia 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 Processes driving the adaptive radiation of a tropical tree (Diospyros, Ebenaceae) in New

http://broadinstitute.github.io/picard/. Broad Institute.

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, 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869

Tillandsioideae) inferred from three plastid DNA markers and the ETS of the nuclear ribosomal DNA. Bot. J. Linn. Soc. 181:362–390.

Poplin R., Ruano-Rubio V., DePristo M.A., Fennell T.J., Carneiro M.O., Auwera G.A.V. der, 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:201178. Quinlan A.R., Hall I.M. 2010. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 26:841–842.

R Core Team. 2020. R: A language and environment for statistical computing. Vienna,

Rivera N.L. 2019. Exploring genomic relationships in Tillandsia subgenus Tillandsia. [Master‘s Thesis]. University of Vienna. htpps://doi.org/10.25365/thesis.57906 Rundell R.J., Price T.D. 2009. Adaptive radiation, nonadaptive radiation, ecological speciation and nonecological speciation. Trends Ecol. Evol. 24:394–399.

Salichos L., Rokas A. 2013. Inferring ancient divergences requires genes with strong phylogenetic signals. Nature. 497: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:354–357. 879

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. Smith, L.B. & Downs, R.J. 1977. Tillandsioideae (Bromeliaceae). Flora Neotrop. 14:798Soltis P.S., Folk R.A., Soltis D.E. 2019. Darwin review: angiosperm phylogeny and evolutionary radiations. Proc. R. Soc. B Biol. Sci. 286: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:e3000391.

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:150. flow. Syst. Biol. 65:843–851.

Solís-Lemus C., Yang M., Ané C. 2016. Inconsistency of species tree methods under gene 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916

Phylogenetic signal detection from an ancient rapid radiation: Effects of noise reduction, long-branch attraction, and model selection in crown clade Apocynaceae. Mol. Phylogenet. Suarez-Gonzalez A., Lexer C., Cronk Q.C.B. 2018. Adaptive introgression: a plant perspective. Biol. Lett. 14:20170688.

Suh A., Smeds L., Ellegren H. 2015. The dynamics of incomplete lineage sorting across the ancient adaptive radiation of neoavian birds. PLOS Biol. 13:e1002224.

Suvorov A., Kim B.Y., Wang J., Armstrong E.E., Peede D., D’Agostino E.R.R., Price D.K., Waddell P.J., Lang M., Courtier-Orgogozo V., David J.R., Petrov D., Matute D.R., Schrider D.R., Comeault A.A. 2022. Widespread introgression across a phylogeny of 155 Drosophila genomes. Curr. Biol. 32:111-123.e5.

Taylor S.A., Larson E.L. 2019. Insights from genomes into the evolutionary importance and prevalence of hybridization in nature. Nat. Ecol. Evol. 3:170–177. 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933

Than C., Ruths D., Nakhleh L. 2008. PhyloNet: a software package for analyzing and reconstructing reticulate evolutionary relationships. BMC Bioinformatics. 9:1–16. Till W. 2000. Tillandsioideae. D.H. Benzing, ed. Bromeliaceae: profile of an adaptive radiation. New York, New York: Cambridge University Press. p. Pp. 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ánchezGonzá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.,Michael H.J. 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

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:599–603.

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. Springer US: pp. 265– 296.

Whitfield J.B., Lockhart P.J. 2007. Deciphering ancient rapid radiations. Trends Ecol. Evol. 22:258–265. 939

Winkler S. 1990. Zur Evolution der Gattung Tillandsia L. Bot. Jahrb. Syst. 112:43–77. 934 935 936 937 938 940 941 942 943 944 945 946 947 948 949 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. syad026. 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. Yan H., Jiang J. 2007. Rice as a model for centromere and heterochromatin research.

Yu G. 2020. Using ggtree to visualize data on tree-like structures. Curr. Protoc. Bioinforma. Zeberg H., Pääbo S. 2020. The major genetic risk factor for severe COVID-19 is inherited from Neanderthals. Nature. 597:610–612.

Zhang C., Rabiee M., Sayyari E., Mirarab S. 2018. ASTRAL-III: polynomial time species tree reconstruction from partially resolved gene trees. BMC Bioinformatics. 19: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:961–975.

Zheng X., Levine D., Shen J., Gogarten S.M., Laurie C., Weir B.S. 2012. A highperformance computing toolset for relatedness and principal component analysis of SNP data.

Bioinformatics. 28:3326–3328.

Yardeni et al. 964 965 966 967 968 969 970 971 ASTRAL-III for 64 individuals representing 34 species of Tillandsia subg. Tillandsia, plus two outgroups. Branch lengths are given in coalescent units. Node values represent local posterior probabilities for the main topology and are equal to one 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 (grey).

Yardeni et al. 972 973 974 975 976 was used as the outgroup in all tests. The four taxa in each test have been rearranged to always obtain positive D values and P2 and P3 are shown on the axes. Colour indicates the value of D and the log value of p-value, the latter estimated using a block jackknife procedure with a 200kb window size and corrected for family-wise error rate. Colours on the bars correspond to the clades in Figure 1.

Pervasive hybridization in radiated Tillandsia 978 979 980 981 982 983 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 the legend of the fb(C) values and grey squares correspond to tests that are inconsistent with the ASTRAL phylogeny. Dots within the matrix denote a significant pvalue, estimated using a block jackknife procedure and corrected for family wise error rate.

to three reticulation events. Curved branches indicate reticulation events. Numbers next to curved branches indicate inheritance probabilities for each event. Colours correspond to the clades in Figure 1.