Introgression across evolutionary scales suggests reticulation contributes to Amazonian tree diversity

Hybridization has the potential to generate or homogenize biodiversity and is a particularly common phenomenon in plants, with an estimated 25% of plant species undergoing interspecific gene flow. However, hybridization in Amazonia's megadiverse tree flora was assumed to be extremely rare despite extensive sympatry between closely related species, and its role in diversification remains enigmatic because it has not yet been examined empirically. Using members of a dominant Amazonian tree family (Brownea, Fabaceae) as a model to address this knowledge gap, our study recovered extensive evidence of hybridization among multiple lineages across phylogenetic scales. More specifically, using targeted sequence capture our results uncovered several historical introgression events between Brownea lineages and indicated that gene tree incongruence in Brownea is best explained by reticulation, rather than solely by incomplete lineage sorting. Furthermore, investigation of recent hybridization using ~19,000 ddRAD loci recovered a high degree of shared variation between two Brownea species that co‐occur in the Ecuadorian Amazon. Our analyses also showed that these sympatric lineages exhibit homogeneous rates of introgression among loci relative to the genome‐wide average, implying a lack of selection against hybrid genotypes and persistent hybridization. Our results demonstrate that gene flow between multiple Amazonian tree species has occurred across temporal scales, and contrasts with the prevailing view of hybridization's rarity in Amazonia. Overall, our results provide novel evidence that reticulate evolution influenced diversification in part of the Amazonian tree flora, which is the most diverse on Earth.


| INTRODUC TI ON
Reproductive isolation is often seen as a prerequisite for speciation and as a defining feature of species (Barraclough, 2019;Mayr, 1942).
Despite this, hybridization between species is known to occur and has several different outcomes, from the erosion of evolutionary divergence (Kearns et al., 2018;Vonlanthen et al., 2012) to the formation of entirely new "hybrid species" (Mallet, 2007). Neotropical rainforests, and in particular Amazonia, harbour the highest levels of plant diversity on Earth (Cardoso et al., 2017;Ulloa Ulloa et al., 2017), the evolution of which was influenced by many drivers (Antonelli & Sanmartín, 2011). However, to date there has been little convincing evidence that hybridization consistently occurs between tree species therein. Indeed, the prevailing view has been that hybridization between tropical tree species is an exceptionally rare event (Ashton, 1969;Ehrendorfer, 1970;Gentry, 1982). Although a few observations of reproductive biology support this, such as the intersterility between Costa Rican lineages of Inga, a species-rich legume genus (Koptur, 1984), evidence for hybridization has been poorly tested empirically in Amazonian trees (in contrast to other Neotropical regions and taxa (e.g. Morales-Briones, Liston, & Tank, 2018)). Indeed, within Amazonia the only potential molecular evidence of hybridization among tree species involved two species of Carapa (Scotti-Saintagne et al., 2013), and a recent review of hybrid zones (which are one potential outcome of hybridization) identified no studies focussed on Amazonian taxa (Abbott, 2017). Since Amazonia is the largest expanse of rainforest in the world, containing at least 6,800 tree species (Cardoso et al., 2017) and given that hybridization is widely acknowledged as a powerful creative force in plant evolution, being known to generate morphological and genetic novelty (Rieseberg et al., 2007), the paucity of studies examining it as an aspect of evolution in the hyperdiverse rainforests of Amazonia is surprising.
One factor which may allow more widespread hybridization in Amazonian tree species than was previously assumed, if their reproductive isolation is not absolute, is the remarkable level of sympatry found for closely related species. In several Amazonian lineages, including Inga, Guatteria and Swartzia, many recently diverged species co-occur (Dexter et al., 2017), often to a remarkable degree. One such example of this is the co-occurrence of 19 Inga species in a single hectare of the Ecuadorian Amazon (Valencia, Balslev, & Miño, 1994). As such, for many rainforest taxa the opportunity for hybridization is constantly present.
Hybridization can have a range of evolutionary consequences.
In many cases, hybridization simply results in the formation of sterile, maladapted offspring with poor reproductive fitness, with subsequent selection against them maintaining isolation between species (such as in reinforcement, where selection against hybrids drives the evolution of prezygotic isolation (Hopkins, 2013)). In other cases, there may be a significant movement of genetic material from one lineage to another facilitated by backcrossing (Rieseberg & Wendel, 1993), which is known as "introgression." This transfer of genetic material through hybridization may confer a selective advantage to resultant offspring (Taylor & Larson, 2019), in which case it is referred to as "adaptive introgression." Adaptive introgression is often observed between closely related taxa during the invasion of new habitats (e.g. Suarez-Gonzalez, Lexer, & Cronk, 2018; e.g. Whitney, Randell, & Rieseberg, 2010). Furthermore, hybridization can lead to rapid evolutionary radiations. This occurs through the re-assembly of standing genetic variation which has accumulated between diverging lineages, and which has already been subject to selection. This "combinatorial" process is much more rapid than the gradual accumulation of variation through mutation, and the passage of these variants through hybridization often fuels rapid diversification events (Marques, Meier, & Seehausen, 2019). This has so far been demonstrated largely in animals, most notably in African cichlid fish (Meier et al., 2017).
Introgression can occur at different rates in different regions of the genome (Payseur, 2010). Regions under divergent selection may remain distinct due to reduced fitness of hybrid genotypes at such loci, resulting in a low rate of introgression, as demonstrated in temperate-zone tree species, where divergence between hybridizing lineages is maintained through selection (Hamilton & Aitken, 2013;Sullivan, Owusu, Weber, Hipp, & Gailing, 2016). This explains why species that hybridize can remain as biologically distinct (and taxonomically identifiable) groups despite undergoing genetic exchange with other species (Abbott et al., 2013;Seehausen et al., 2014). Conversely, regions under little or no selection tend to introgress more freely, becoming homogenized between species.
Moreover, if there is selection for hybrid genotypes (as in adaptive introgression), the rate of introgression for a region may be further increased relative to the rest of the genome (Gompert, Parchman, & Buerkle, 2012).
Brownea, a member of the legume family (Fabaceae), is a characteristic tree genus of lowland Neotropical rainforests and contains around 27 species distributed across northern South America, particularly within Amazonia. Previous work indicates that there is a broad degree of phylogenetic incongruence evident in this genus (Schley et al., 2018) which might indicate hybridization, although this could also result from incomplete lineage sorting. However, there are numerous Brownea hybrids in cultivation (e.g. B. x crawfordii (Crawford & Nelson, 1979) and B. hybrida (Backer, 1911)), as well as several instances of putative hybridization among Brownea lineages in the wild. The most notable of these instances is that proposed between the range-restricted Brownea jaramilloi (Pérez, Klitgård, Saslis-Lagoudakis, & Valencia, 2013) and the wide-ranging B. grandiceps, which co-occur in the Ecuadorian Amazon ( Figure 1). There are multiple morphological distinctions between these two species, including differences in inflorescence colour and structure, growth habit, tree height and leaf morphology (Pérez et al., 2013). Although they co-occur, these species favour different habitats: B. jaramilloi grows on ridge tops and hillsides, whereas B. grandiceps shows a slight preference for swamps and valleys but is more evenly distributed (pers. obs. & Klitgaard, 1991;Pérez et al., 2013). Despite this, hybridization appears to occur, as evidenced by the existence of a putative hybrid between these two species known as B. "rosada" (Figure 1). Brownea "rosada" displays an intermediate morphology between its two parental species, producing pink flowers. The hypothesis of a B. jaramilloi x B. grandiceps hybrid has not yet been tested, however, using molecular data.
As a member of the legume family, which dominates Amazonian forests (ter Steege et al., 2013), and with its apparent propensity for hybridization, Brownea is an excellent system with which to study the phylogenetic patterns and genomic architecture of introgression in Amazonian trees. Systematically documenting hybridization at a range of timescales and taxonomic levels within this group could help understand how admixture has contributed to the assembly of one of the world's richest floras. Accordingly, this study investigates the role of hybridization in the diversification of this rainforest tree genus at two levels. First, we test whether introgression played a role during the long-term diversification of the genus, by investigating reticulation at deep phylogenetic levels. Second, we test whether recent gene flow between Brownea species occurs at the landscape scale and if so, whether it occurs evenly across the genome or differentially for subsets of loci.

| Phylogenetic analyses
In order to test for ancient reticulation between Brownea species, sequences from 225 nuclear genes were used to elucidate evolutionary relationships. This was done using a targeted bait capture sequencing approach along with phylogenomic methods, for which 23 of 27 lineages were sampled within Brownea, including the three subspecies of B. coccinea and B. "rosada," a putative hybrid lineage of B. grandiceps and B. jaramilloi. In total, 59 accessions were sampled within Brownea, as well as an additional 13 outgroup taxa from the genera Macrolobium, Heterostemon, Paloue and Browneopsis, which form part of the "Brownea clade" (Fabaceae, subfamily Detarioideae (LPWG, 2017)). The list of accessions and their associated information can be found in Table S1.
DNA sequence data were generated using leaf material collected from herbarium specimens and silica-gel-dried accessions. Genomic DNA was extracted using the CTAB method (Doyle & Doyle, 1987), and sequencing libraries were prepared using the NEBNext® Ultra™ II DNA Kit (New England Biolabs, Massachusetts, USA) with a modified protocol to account for fragmented DNA, and including a ~600 bp size-selection step. Targeted bait capture was performed using the MyBaits protocol (Arbor Biosciences, Michigan, USA) to target 289 phylogenetically informative nuclear genes, using a bait kit designed for the legume subfamily Detarioideae within which Brownea is nested (Ojeda et al., 2019 DNA sequencing reads were quality-checked with the program FASTQC v0.11.3 (Andrews, 2010) and were subsequently trimmed using Trimmomatic v.0.3.6 (Bolger, Lohse, & Usadel, 2014) in order to remove adapter sequences and to quality filter reads.
Trimmomatic settings permitted <4 mismatches, a palindrome clip threshold of 30 and a simple clip threshold of 6. Bases with a quality score <28 and reads shorter than 36 bases long were removed from the data set. Following quality filtering, loci were assembled using SPAdes v3.11.1 (Bankevich et al., 2012) by the HybPiper pipeline v1.2 (Johnson et al., 2016) and potentially paralogous loci were removed using the Python (Python Software Foundation, 2010) script "paralog_investigator.py," distributed with the HybPiper pipeline. All sequences were aligned by gene region using MAFFT v7.215 (Katoh & Standley, 2013) and were cleaned using the -automated1 flag in trimal (Capella-Gutiérrez, Silla-Martínez, & Gabaldón, 2009), which calculates the optimal method to remove sequence gaps from alignments and removes them. In order to infer gene trees for phylogenetic network analysis, the 225 recovered gene regions were further refined to include only 20 taxa, representing a single accession per lineage in Brownea, F I G U R E 1 Two co-occurring Brownea lineages (Brownea grandiceps (photograph © Rowan Schley) and Brownea jaramilloi (photograph © Xavier Cornejo)) as well as their putative hybrid Brownea "rosada" (photograph © J. L. Clark) [Colour figure can be viewed at wileyonlinelibrary.com] using Macrolobium colombianum as the outgroup taxon. Where applicable, accessions were chosen by comparing the sequence recovery of conspecific samples and choosing the individual with the best gene recovery. This resulted in 220 single-accession-per-lineage gene alignments. Further details of phylogenomic sampling, labwork, sequencing and data processing are detailed in Methods S1.
Inference was performed using 1,000 rapid bootstrap replicates and the GTRCAT model of nucleotide substitution (which is a parameter-rich model and hence is most suitable for large data sets) on the CIPRES web portal ( (Miller, Pfeiffer, & Schwartz, 2010), https:// www.phylo.org/). In addition, gene trees were generated for each of the full and single-accession gene alignments using RAxML v.8.0.26 (Stamatakis, 2014) with the same settings as above. Macrolobium colombianum was used to root both the full and single-accession data set analyses since it was the Macrolobium accession that was present in the most gene alignments.
Gene trees from both data sets were used to generate species trees under the multispecies coalescent model in the heuristic version of ASTRAL v.5.6.1 (Zhang, Rabiee, Sayyari, & Mirarab, 2018) using the default parameters. Monophyly was not enforced for individuals belonging to the same species (the "-a" option in ASTRAL).

| Inferring ancient introgression
Phylogenetic networks were inferred for 220 gene trees from the single-accession-per-lineage data set to understand whether introgression occurred over Brownea's evolutionary history. Networks Networks were estimated by calculating quartet concordance factors (CF) for each node from the single-accession-per-lineage gene trees, since SNaQ! requires that each tip of the phylogenetic trees represents a single lineage. The network with the number of hybridization events (h) best describing the data was chosen using negative log-pseudolikelihood comparison. Finally, the observed gene tree CFs were compared with those of the best-fit phylogenetic network (i.e. the model including hybridization) and those expected under a coalescent model (i.e. a "tree-like" model which accounts for incomplete lineage sorting but not hybridization). The fit of the observed gene tree topologies to either the "network-like" or the "tree-like" model was assessed using the Tree Incongruence Checking in R (TICR) test (Stenz, Larget, Baum, & Ané, 2015). This was done with the function "test.one.species.tree" in the R v3.4.4 (R Development Core Team, 2013) package PHYLOLM (Ho & Ané, 2014). Proportions of genes contributed between lineages by hybridization events were taken from the "Gamma value" output of PhyloNetworks. Further details of phylogenetic network analysis are contained in Methods S1.
In order to further investigate patterns of ancient reticulation, we used the Dtrios function in Dsuite (Malinsky, 2019) to perform Patterson's D-statistic tests (Durand, Patterson, Reich, & Slatkin, 2011;Green et al., 2010). Patterson's D-statistic uses asymmetry in gene tree topologies to quantify introgression between either of two lineages which share a common ancestor (P1 and P2) and one other lineage (P3) that diverged from the common ancestor of P1 and P2 or earlier. We generated an input VCF file for Dtrios based on the single-accession-per-lineage data set using SNPsites (Page et al., 2016) and ran Dtrios with Macrolobium colombianum as an outgroup for all 969 tests. In addition, we used the single-accession-per-lineage ASTRAL tree to inform taxon relationships for each test. We assessed significance of each test using 100 jackknife resampling runs and visualized the most conservative D-statistic estimates (the "D_min" output of Dtrios) with the Ruby script "plot_d.rb," available from https://github.com/mmats chiner.

| Population-level introgression between B. grandiceps and B. jaramilloi
Having examined the degree of historical reticulation among Brownea lineages using targeted bait capture, population genomic data were generated with ddRAD sequencing and used to investigate recent introgression at a finer taxonomic scale. More specifically, the degree of shared genetic variation was visualized for individuals of the ecologically divergent species B. grandiceps and B.
jaramilloi focussing on individuals which co-occur in Yasuní National Park in the Ecuadorian Amazon, as this area appears to be a hybrid zone for these two taxa. Following this, to make inferences about the potential evolutionary significance of recent introgression the rates at which different loci introgress relative to the rest of the genome were estimated for these taxa using genomic clines.
One hundred and seventy-one specimens in total were genotyped using ddRADseq (Peterson, Weber, Kay, Fisher, & Hoekstra, 2012). jaramilloi trees was collected from the Yasuní National Park 50-ha forest plot in Napo, Ecuador, and dried in a herbarium press. Samples were identified based on the characters listed in Table S2, with these morphological differences between B. grandiceps and B. jaramilloi displayed in Figure S1. The sample list is shown in Table S3.
Genomic DNA was extracted from dried leaf material with the CTAB method and purified using a Qiagen Plant Minikit (Qiagen, Hilden, Germany) column cleaning stage. Sequencing libraries were prepared by digesting the DNA template using the restriction enzymes EcoRI and mspI (New England BioLabs, Massachusetts, USA), in accordance with the ddRADseq protocol. Samples were then ligated to universal Illumina P2 adapters and barcoded using 48 unique Illumina P1 adapters. Samples were pooled and size-selected to between 375 and 550 bp using a Pippin Prep Electrophoresis machine (Sage Science, Massachusetts, USA). Following amplification and multiplexing, sequencing was performed using a single-lane, paired-end 150bp run on the HiSeq 3/4000 platform.
DNA sequencing reads from the ddRADseq genotyping were processed de novo (i.e. without the use of a reference genome) using the Stacks pipeline v2.1 (Catchen, Amores, Hohenlohe, Cresko, & Postlethwait, 2011). Reads were quality-filtered using Stacks by removing reads which were of poor quality (i.e. had a Phred score <10), following which "stacks" of reads were created, and SNPs were identified among all de novo "loci" and across individuals. This was done using a minimum coverage depth (the "-m" flag in Stacks) of three and a within-individual mismatch distance (-M) of seven nucleotides.
Individuals with sequencing coverage under 7.5x were removed. Loci found in fewer than 40% of individuals and sites with a minor allele frequency threshold of 5% were filtered out using the "populations" In order to understand the patterns of introgression at the population level, the single-SNP ddRAD data set containing 19,130 RAD loci was used to visualize the degree of shared variation between B. grandiceps and B. jaramilloi. This was performed using principal component analysis implemented in R v3.4.4 followed by plotting with gg-plot2 (Wickham, 2016). To further understand the patterns of shared variation, we visualized relationships between individuals by inferring a neighbour net plot in SPLITSTREE v4.14.6 (Huson & Bryant, 2005) and inferred population structure with the program fastSTRUCTURE v1.0 (Raj, Stephens, & Pritchard, 2014). The number of lineages "K" (hereafter called "groups," as is standard in FastSTRUCTURE) was chosen using the value which provided the largest improvement in marginal likelihood. In addition, a fastSTRUCTURE run incorporating 40 random individuals from each species was performed to account for any bias which may have been incurred by differences in sample size. Finally, NEWHYBRIDS v1.1 (Anderson & Thompson, 2002)(https://github. com/eriqa nde/newhy brids) was used to categorize individuals into different hybrid classes, using three runs on a subset of 500 randomly selected ddRAD loci due to the computational limits of the program.
Runs were performed with 50,000 MCMC sweeps following 50,000 burn-in sweeps under the default parameters of the program.
For each locus, bgc compares the probability of ancestry at the locus relative to an individual's genome-wide ancestry, thereby allowing it to estimate two parameters for each locus. These parameters are α, which approximately equates to the "direction" of introgression, and β, which may be summarized as the "rate" of introgression for a locus (Gompert & Buerkle, 2011;Gompert, Parchman, et al., 2012).
In order to estimate these parameters from the single-SNP data set consisting of 19,130 RAD loci, 50,000 MCMC generations with a 50% burn-in were performed in bgc. Admixture proportions (i.e. mean Q values) generated by fastSTRUCTURE were used to assign each individual to one of three groups (Brownea grandiceps, B. jaramilloi and admixed), resulting in an admixed population containing 27 individuals. In addition, hybrid index and interspecific heterozygosity were estimated with bgc using the -p and -i flags, respectively, and were plotted using the triangle.plot function in the R package Introgress (Gompert & Buerkle, 2010).
Convergence was checked for the MCMC output from bgc in Tracer v1.6 (Rambaut, Suchard, Xie, & Drummond, 2015) and with the R package coda (Plummer, Best, Cowles, & Vines, 2006) using Geweke's diagnostic (Geweke, 1991). Loci with significant "excess ancestry" were identified by ascertaining whether the 99% posterior probability estimates of the α and β parameters included zero (i.e. by identifying positive or negative nonzero estimates of the parameters). In addition, loci which were extreme "introgression outliers" were identified for both parameters by extracting loci whose median estimates were not included in the 99% posterior probability credible intervals (Gompert & Buerkle, 2011).
Finally, to assess whether our bgc analysis was limited in its ability to detect introgression outliers due to sample size, we ran an extra set of analyses based on a randomly subsampled data set consisting of 9,430 RADseq loci. These data were taken from a similar plant speciation study (Royer, Streisfeld, & Smith, 2017) and were refined to contain the same number of admixed individuals as our study (n = 27). We then ran the bgc analyses with the parameters and methods described above. Further details of visualizing shared variation, hybrid category assignment and bgc analysis are shown in Methods S2.

| Phylogenetic analyses
Using DNA sequence data from 225 nuclear gene regions, we produced a concatenated phylogenetic tree using RAxML (Figure 2). This resulted in well-supported relationships between most bipartitions (>90% bootstrap support (BS)), with lower support for some interspecific relationships, especially within the "Grandiceps" subclade.
This tree also indicated a very strong signal of geographical structure among species, with several wide-ranging species (such as Brownea grandiceps) being polyphyletic and being more closely related to accessions from the same region than to conspecifics. The ASTRAL tree run using the same data set resulted in a nearly identical topology F I G U R E 2 RAxML phylogenetic tree inferred from 220 concatenated genes. Numbers at nodes of the tree signify bootstrap support (BS). Taxa within the top box belong to the "Grandiceps" subclade, taxa within the next box down box belong to the "Coccinea" subclade, taxa within third box down represent the "Stenantha" subclade, and those within the bottom box belong to the "Chocóan" subclade. The province and country from which each accession was collected are noted in parentheses next to each tip label, and colours correspond to the inset map (top left). Green corresponds to Western Amazonia, blue is the Chocó-Darién region and Andean cordilleras, while orange is Eastern Amazonia and the Guiana shield (modified from the phytogeographical regions identified by Gentry (1982)  Brownea grandiceps S100 (Napo, Ecuador)

Browneopsis cauliflora S88 (San Martín, Peru)
Brownea grandiceps S34 (Ecuador) with similar levels of geographical structure and inferred a high degree of discordance among gene tree topologies at many nodes (quartet score = 0.52), with multiple alternative bipartitions reconstructed at most nodes. This is evident from the presence of many more conflicting gene trees (numbers below branches) than congruent gene trees (numbers above branches) in Figur S2.

| Ancient introgression in Brownea species
Phylogenetic networks estimated to ascertain whether diversification in Brownea was tree-like or reticulated pinpointed two hybridization events within the genus (Figure 3a). This is indicated by the -log-pseudolikelihood values in Figure S3, since the number of hybridization events (h) which best describe the data gives the largest improvement in -log pseudolikelihood. In Figure S3, -log pseudolikelihoods increased steadily between h = 0 and h = 2, after which the increasing number of hybridization events only made minimal improvements to -log pseudolikelihood.

| Population-level introgression between B. grandiceps and B. jaramilloi
Our population genomic analyses revealed a broad degree of shared variation between B. grandiceps and B. jaramilloi, along with a homogeneous distribution of introgression rates across loci. In total, ddRAD sequencing resulted in 640,898,950 reads between 350 and 550bp in length for 171 individuals. Among these reads, 28,503,796 (4.45%) were discarded due to poor recovery. There was a mean  Figure S5a), which was the value that resulted in the largest increase in marginal likelihood. Since the "best" value of K is an estimate, fastSTRUCTURE plots generated with other K values are shown in Figure S5b). Figure 4b shows that most individuals of B. jaramilloi have varying fractions of B. grandiceps ancestry in their genome. In addition, the individuals identified as B. "rosada" appear to have inherited most of their ancestry from B. grandiceps, with only a minimal contribution from B. jaramilloi.
The same pattern was recovered from a fastSTRUCTURE run incorporating 40 random individuals from each species, performed to account for any bias which may have been incurred by differences in sample size ( Figure S5c).

Brownea ariza S90 (Magdalena, Colombia)
Brownea rosa-de-monte S66 (Antioquia, Colombia) suggesting backcrossing. In the latter, the probability (Q) of samples belonging to any of the two species was below 90%; otherwise, they were identified as either B. grandiceps or B. jaramilloi.

Most patterns of shared variation inferred by Splitstree,
FastSTRUCTURE and NEWHYBRIDS were congruent with the F I G U R E 3 a) Phylogenetic network with two hybridization events (h = 2), estimated using SNaQ! in the Julia package PhyloNetworks. Light blue horizontal branches indicate inferred hybridization events, and numbers next to the branches show the estimated proportion of genes contributed by each lineage in the hybridization event. Red branches signify the ancestral lineage and what proportion of the modern lineage's genes were contributed from it. Taxa were chosen in order to represent one accession per lineage, as per the assumptions of PhyloNetworks, and the individual with the highest gene recovery was used for each lineage. Taxa within the top box belong to the "Grandiceps" subclade, taxa within the next box down box belong to the "Coccinea" subclade, taxa within third box down represent the "Stenantha" subclade, and those within the bottom box belong to the "Chocóan" subclade. The province and country from which each accession was collected are noted in parentheses next to each tip label, and colours correspond to the inset map (top left). Green corresponds to Western Amazonia, blue is the Chocó-Darién region and Andean cordilleras, while orange is Eastern Amazonia and the Guyana shield (modified from the phytogeographical regions identified by Gentry (1982)). b) Heatmap summarizing the most conservative four-taxon D-statistic estimates and their P-values from 969 tests. Taxa P2 and P3 are displayed on the x-and y-axes in the same order as in Figure 3a, and each square represents the highest estimate of each combination of P2 and P3 taxa. The colour of each square signifies the D-statistic estimate (blue = low estimate; red = high estimate), and the colour saturation of these colours represents the P-value for that test (see inset box, bottom right). D-statistic tests for which p < 1x10 -6 (i.e. were below a Bonferroni-adjusted P-value threshold of 0.001) are marked with an asterisk (*). As in Figure 3a, the province and country from which each accession was collected are noted in parentheses next to each tip label, and colours correspond to the inset map ( (Table S5). Splitstree was the most congruent (92.98% of accessions with the same classification as identified morphologically), followed by NEWHYBRIDS (92.31%) and FastSTRUCTURE (85.96%). The same hybrid class assignment was recovered for 84.62% of accessions across methods.
The Bayesian estimation of genomic clines (bgc) analysis, which was used to determine whether loci showed outlying rates of admixture relative to the genome-wide average, recovered a signal of asymmetric introgression, but did not detect any differential rates of introgression between loci. Of the 19,130 loci under study, bgc recovered 251 loci (1.3%) with positive excess ancestry for the α parameter and 1,089 loci (5.69%) with negative excess ancestry for α.

| D ISCUSS I ON
Our study represents the first clear case of reticulate evolution among Amazonian trees documented at multiple phylogenetic levels. This contrasts with previous work, which suggested that hybridization is an extremely rare phenomenon in Amazonian trees (Ashton, 1969;Ehrendorfer, 1970;Gentry, 1982). We demonstrate that within Brownea, a characteristic Amazonian tree genus, reticulate evolution has occurred over the course of its evolutionary history, with evidence of hybridization deep in the history of the genus and more recently, between the closely related B. grandiceps and B. jaramilloi.
Our study recovered a low degree of support for bifurcating interspecific relationships ( Figure 2) and a high degree of incongruence among gene trees ( Figure S2). These results can be largely explained by the strong signal of geographically structured reticulation that we recovered, rather than being caused purely incomplete lineage sorting, which may be common in rainforest trees (Pennington & Lavin, 2016). Indeed, the low quartet scores and minimal gene tree concordance found at the species level in this study mirror those observed in Lachemilla, a montane Neotropical plant genus, in which gene tree discordance was shown to be explained by both historical and recent hybridization (Morales-Briones et al., 2018), both of which we found evidence for.

| Reticulation has occurred at deep phylogenetic levels in Brownea
Phylogenetic network analysis suggested that reticulation has taken place between Brownea lineages in the past, with two separate hybridization events inferred (Figure 3a). The concordance factors obtained from gene trees were not adequately explained by a tree-like model (i.e. one accounting only for incomplete lineage sorting), suggesting that a network-like model (i.e. one including reticulation) best describes the diversification patterns within Brownea (P = 1.240 × 10 -19 , Χ 2 = 91.149). Since the inferred events of ancient introgression both involved the common ancestors of two sets of present-day species, they are likely to have occurred several million years ago, before the divergence of the descendent species. These introgression events are likely to have occurred since the Miocene period (~23Ma), as date estimates from previous work suggest that Brownea originated around this time (Schley et al., 2018). Thus, the detected signatures of ancient reticulation suggest that their effects in genomes have persisted over evolutionary time.
A distinct pattern of geographical structure was evident from our phylogenetic tree and network (Figure 2; Figure 3a). Wide-ranging species such as B. grandiceps and B. coccinea were polyphyletic, and co-occurring individuals from different species were more closely related to each other than to morphologically identified conspecifics.
While this polyphyly could be due to stabilizing selection on certain aspects of the genome constraining morphological divergence while the rest of the genome diverges (Struck et al., 2018), this would not explain the high degree of reticulation between co-occurring species that we inferred. We found that both historical introgression events reconstructed by PhyloNetworks were between taxa that co-occurred (e.g. B. grandiceps, B. leucantha and B. birschelli in Northern Venezuela, and B. coccinea, B. santanderensis and B. enricii within the Colombian cordilleras), which was in agreement with our D-statistic tests (Figure 3b, discussed below). This geographically structured reticulation agrees with previous work (Schley et al., 2018) which indicated that the "stem" lineages within Brownea had a shared distribution in northern South America, perhaps facilitating hybridization.
Geographical structure is common in species complexes of tropical trees, where geography often correlates more closely with genetic affinity between related species than does morphology (Dexter, Pennington, & Cunningham, 2010). As alluded to above, this pattern can result from gene flow between sympatric, interfertile congeners which together form a syngameon ). Combined with the poor dispersal likely in Brownea given its explosively dehiscent fruit pods and tendency to form single-species stands (Klitgaard, 1991), we believe that the syngameon model best describes the reticulate evolutionary history of this group. It is also possible that our results reflect a similar model where widespread species (such as B. grandiceps) are progenitor lineages from which new, localized species diverge in peripatry through ecological speciation (Misiewicz and Fine, 2014;Pennington & Lavin, 2016) and subsequently hybridize in secondary contact. However, the latter seems unlikely given that relatively distantly related, co-occurring lineages show evidence of hybridization events.
The syngameon model seems also to describe the reticulate patterns of diversification in Brownea when considering the results of our D-statistic tests (Figure 3b). These tests also recovered a geographically structured pattern of reticulation with three main centres: Ecuador, the Colombian Cordilleras and Venezuela. Geography has been shown to be a factor structuring syngameons, based simply on the increased opportunity of interbreeding when in sympatry (Boecklen, 2017). Two other significant introgression events were also inferred: one between B. "rosada" PhyloNetworks can only detect single hybridization events at each node rather than those with multiple participant lineages. In addition, hybridization events between sister lineages can only be detected if both lineages give rise to another two descendent species (Nevado, Contreras-Ortiz, Hughes, & Filatov, 2018)  hybrids.

| Recent hybridisation also occurs between
The higher number of variant sites and the higher nucleotide diversity (π) recovered for B. jaramilloi (Table S4) could also reflect the hybrid ancestry of many individuals identified as this species. A similar introgression-driven increase in nucleotide diversity has been shown in closely related species of Mimulus which undergo asymmetric introgression (Sweigart & Willis, 2003).
The bgc analysis indicated that most introgression was asymmetric and predominated by B. grandiceps, with more loci exhibiting negative excess ancestry for the α parameter ( to "pollen swamping," whereby pollen transfer from another, more populous species is more likely than pollen transfer from less numerous conspecifics (Buggs & Pannell, 2006). Pollen swamping can serve as a mechanism of invasion (e.g. in Quercus (Petit, Bodénès, Ducousso, Roussel, & Kremer, 2004)), and the associated asymmetric introgression may lead to the extinction of the rarer species, especially when hybrids are inviable or sterile (Balao, Casimiro-Soriguer, García-Castaño, Terrab, & Talavera, 2015). However, due to the ongoing introgression evident between B. grandiceps and B.
jaramilloi it appears that hybrids are not always sterile, or at least that "foreign" alleles are not always subject to strong purifying selection at the landscape scale.
The best evidence of this lack of selection against hybrids is the absence of loci with extreme values of the β parameter recovered by our bgc analysis. Extreme deviations in β are mainly expected in the presence of gene flow when selection against hybrids is strong, resulting in underdominance and a paucity of heterozygous sites (Gompert, Parchman, et al., 2012). As such, our inferred lack of outlying β values suggests persistence of hybrid genotypes within the Yasuní 50-ha plot. Similarly, the 27 admixed individuals in this study recovered a range of hybrid indices and low interspecific heterozygosity (Fig. S6), which suggests that these individuals are the result of many generations of interspecific gene flow, resulting in "asymmetric advanced-generation introgression," as found in previous work (De La Torre, Ingvarsson, & Aitken, 2015), although this may also result from the low Fst values we recovered (as in Gompert, Lucas, et al., 2012).
Further to this, the observed asymmetry in introgression could be caused by adaptive introgression favouring the preferential passage of certain loci from one parental species (e.g. Yang et al., 2018).
However, the lack of outlying β estimates for any loci renders this very unlikely, and even so, it is difficult to ascertain whether adaptive introgression has occurred without measuring the impact of this introgression on variation in the phenotype and its fitness effects .
The three methods used to identify lineages and hybrid classes (Splitstree, FastSTRUCTURE and NEWHYBRIDS) largely recovered the same identification as was made using morphology (Table S5) (Atchison, 1951). Moreover, the bias for clustering by ploidy level uncovered by Stift et al. (2019) was most evident when population differentiation was weak, rather than between species as in the case of our study. More generally, estimating ancestry proportions using STRUCTURE analysis can be affected by demographic factors (e.g. recent bottlenecks or ancient population structure) which may be the case for our data, and so, interpretation of their results should take into account that these types of analyses try to parsimoniously explain patterns of shared variation among individuals rather than explicitly test between evolutionary models (Lawson, Van Dorp, & Falush, 2018). Similarly, the fact that a reduced data set of 500 (out of a possible 19,130) SNPs was used in the NEWHYBRIDS analysis means that we may have subsampled certain loci in such a way that they were not representative of some individual's genetic ancestry, giving different conservative hybrid class assignments.
The discrepancies may also have had a biological cause. For example, differential inheritance or selection on regions of the genome responsible for morphological features, when compared to the rest of the genome, may have resulted in individuals that were identified as one species morphologically, but which in actuality had hybrid ancestry (e.g. the many hybrid individuals identified as B. jaramilloi) (Abbott et al., 2013;Seehausen et al., 2014). Given the highly divergent phenotypes of B. grandiceps and B. jaramilloi ( Figure S1; Table S2.), it is unlikely that these discrepancies were purely due to poor identification.

| The contribution of hybridization to Neotropical rainforest tree diversity
Studies such as ours that document persistent hybridization across evolutionary time between tree species within Amazonia, and within tropical rainforests in general, are rare. Indeed, it was previously suggested that interspecific hybrids between rainforest tree species are poor competitors and that fertile hybrid populations are nearly nonexistent (Ashton, 1969). While there is some evidence of introgression in tropical trees from within Amazonia (Scotti-Saintagne et al., 2013), most available studies substantiating it are based on trees from other tropical regions or habitats (e.g. Shorea in Asia (Kamiya et al., 2011;Kenzo et al., 2019), or Rhizophora in Indo-Pacific mangroves (Lo, 2010)). Many of these other instances appear to occur only in degraded habitats or involve infertile first-generation hybrids with minimal backcrossing, which contrasts with the findings of our study. Within Brownea, we uncovered introgression across taxonomic, spatial and temporal scales, with evidence of backcrossing and the persistence of hybrid genotypes at the landscape scale.
It is possible that a lack of selection against hybrids allows the passage of variation between Brownea species, which may persist over evolutionary time and could explain why we found evidence of introgression at both macroevolutionary and microevolutionary scales. As mentioned above, this is congruent with a "syngameon" model, where large, disparate populations of individuals belonging to multiple interfertile species occasionally hybridize and exchange genetic variants. This can aid in the maintenance of genetic diversity in species which exist in very low population densities, as is typical of tropical rainforest trees, thereby preventing Allee effects (Cannon & Lerdau, 2015). This has been shown in temperate tree species , although it typically occurs at the edge of species ranges and its prevalence in hyperdiverse systems such as tropical rainforests is a topic of recent debate (e.g. Cannon & Lerdau, 2019).
While it is difficult to determine whether hybridization between Amazonian tree species occurs within many syngameons of closely related lineages, or whether it is a tendency unique to certain genera such as Brownea, it may be prudent to review how relationships among tropical tree lineages are inferred.
Accordingly, the use of phylogenetic networks in resolving the relationships between such groups may provide additional insight into whether reticulate evolution has contributed to diversification within Amazonian rainforest, which is among the most species-rich environments on Earth (Burbrink & Gehara 2018; Solís-Lemus et al., 2017).

ACK N OWLED G EM ENTS
Laboratory work, fieldwork and herbarium visits to NY and US were funded by the NERC SSCP DTP and a Genetics Society Heredity Fieldwork grant funded part of the fieldwork in Ecuador.
Phylogenomic sequencing was also part-funded by a grant from the

AUTH O R CO NTR I B UTI O N S
This study was conceived and all analyses were performed by R.J.S., and the study was supervised by T.G.B., F.F. and B.K.
Population genomic data were generated by R.J.S., and phylog-

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available from online repositories. All raw reads generated with the targeted bait capture and ddRADseq methods are available on the NCBI