Chapter 3: Topography drives microgeographic adaptations of closely-related species in two tropical tree species complexes

Submitted the 16th of August 2020 and under review in New Phytologist

Sylvain Schmitt\(^1\) , Niklas Tysklind\(^2\), Bruno Hérault\(^{3,4,5}\) , Myriam Heuertz\(^1\) ,

\(^1\) Univ. Bordeaux, INRAE, BIOGECO, 69 route d’Arcachon, 33610 Cestas France; \(^2\) INRAE, UMR EcoFoG (Agroparistech, CNRS, Cirad, Université des Antilles, Université de la Guyane), Campus Agronomique, 97310 Kourou, French Guiana; \(^3\) CIRAD, UPR Forêts et Sociétés, Yamoussoukro, Côte d’Ivoire; \(^4\) Forêts et Sociétés, Univ Montpellier, CIRAD, Montpellier, France; \(^5\) Institut National Polytechnique Félix Houphouët-Boigny, INP-HB, Yamoussoukro, Côte d’Ivoire;

Summary

• Closely-related tree species that grow in sympatry are abundant in rainforests. However, little is known of the eco-evolutionary processes that shape their niches and allow their local coexistence. Here, we assessed genetic species delimitation in closely-related sympatric species belonging to two Neotropical tree species complexes and investigated their genomic adaptation to a fine-scale topographic gradient with associated edaphic and hydrologic features.
• Combining LiDAR-derived topographic data, comprehensive tree inventories, and single nucleotide polymorphisms (SNPs) from gene capture experiments, we explored genome-wide population genetic structure, covariation of environmental variables, and genotype-environment association to assess microgeographic adaptations to topography within the species complexes Symphonia (Clusiaceae), and Eschweilera (Lecythidaceae) with three species per complex and 385 and 257 individuals genotyped, respectively.
• Within species complexes, closely-related tree species had different realized optima for topographic niches defined through the topographic wetness index or the relative elevation, and species displayed genetic signatures of adaptations. Symphonia species were differentially adapted to water and nutrient distribution, whereas Eschweilera species avoided hydromorphic soils and were differentially adapted to soil chemistry.
• Our results suggest that varied topography represents a powerful driver of tropical forest biodiversity, driving differential adaptations and stabilizing local coexistence of closely-related tree species within tree species complexes.

Keywords

Ecological niche; French Guiana; relative elevation; species coexistence; syngameon; topographic wetness index; tropical forests

Introduction

Closely-related tree species are abundant at the local and regional scales in Neotropical forests (Gentry 1988, Steege et al. 2013, Pinheiro et al. 2018). The five thousand tree species botanically identified in Amazonia belong to only height hundred ten genera (Steege et al. 2013), many of which thus contain multiple related species. Even at the local scale, tropical forests shelter up to several hundred tree species per hectare (Gentry 1988), including genera with multiple closely-related species growing in sympatry (Caron et al. 2019). Closely-related species are expected to share similar niche and functional strategies due to phylogenetic constraints (Wiens et al. 2010). Such similar strategies can lead to increased competition between species, and a heightened risk of local competitive exclusion (Turcotte and Levine 2016). The abundance of sympatrically growing closely-related tropical tree species represents thus an intriguing paradox that warrants the investigation of the eco-evolutionary forces that shape their ecological niches and allow their local coexistence.

Genetic studies have shown that species in species-rich genera have a higher level of genetic polymorphism (Caron et al. 2019) and hybridize more often between congeners (Whitney et al. 2010) than those in species-poor genera. Species-rich genera thus frequently contain species complexes, defined as clades of morphologically similar species and/or species that share large amounts of genetic variation due to recent common ancestry and/or hybridization (Pernès and Lourd 1984). Species-rich genera also frequently contain syngameons, defined as species connected by limited, but recurrent, interspecific gene flow (Suarez-Gonzalez et al. 2018). Syngameons are evolving at two contrasting taxonomic levels, they maximize each species adaptation to its ecological niche, and they evolve at the syngameon level which benefits all the constituent species through adaptive introgressions. The syngameon is thus unlikely to just represent a transitional stage before complete speciation; instead it is probably a successful evolutionary strategy that maximizes adaptive evolution while decreasing the overall risk of extinction (Cannon and Lerdau 2015, Cannon and Petit 2019). Species complexes and syngameons frequently include sympatric species coexisting at local scales (Gonzalez et al. 2009, Caron et al. 2019, Cannon and Petit 2019), but how these species remain distinct and persist in sympatry despite phylogenetic constraints and inter-specific gene flow remains poorly known.

The local coexistence of ecologically and phenotypically similar, closely-related species is governed by ecological and evolutionary processes, contingent on the history of speciation and niche differentiation (Weber and Strauss 2016). Once in contact, ecologists consider foremost the competition for resources and the resulting niche differentiation (Chesson 2000a, Turcotte and Levine 2016), whereas evolutionary biologists are more interested in the level of reproductive isolation (Weber and Strauss 2016). For stable coexistence, niche theory evokes that species maximize ecological niche differences (Weiher and Keddy 1995, Lortie et al. 2004b), while the “emerging similarity theory” reconciles niche and neutral theories by suggesting coexistence of distinct groups of species that are functionally similar within groups (Scheffer and Van Nes 2006). Functionally similar species must have similar fitness as an additional condition for stable coexistence (Chesson 2000a), and recently diverged species must have evolved sufficient reproductive isolation to avoid the break-down of differences and their genetic homogenization upon secondary contact (Levin et al. 1996, Taylor et al. 2006, Abbott et al. 2013).

The heterogeneity of resource distribution in space and time defines fine-scale habitat structure where species can coexist. In tropical forests, topography explains the local spatial distribution of water and nutrient availability, showing a strong association with soil nitrogen, carbon, and phosphorus content (Ferry et al. 2010b). Topography has been shown to explain pervasive differentiation in habitat preference among species (Gunatilleke et al. 2006, Engelbrecht et al. 2007, Kraft et al. 2008, Allié et al. 2015). In particular, soil nutrients, also influenced by topography through hydromorphy, directly shape the spatial distribution of forest tree species (John et al. 2007). Topography has also been shown to drive functional responses among and within species (Schmitt et al. 2020). Topography, for example measured through the topographic wetness index (Schmitt et al. 2020) or relative elevation (Allié et al. 2015), is thus a proxy of numerous habitat features related to the distribution of nutrients and water. Therefore, fine-scale mapped topography is a good candidate factor that may explain adaptation and local coexistence of closely-related species in Neotropical tree species complexes.

Closely-related species growing in sympatry in differentiated ecological niches are frequently the product of an adaptive radiation (Seehausen 2004), such as Darwin’s finches in the Galapagos (Grant and Grant 2019). Evolutionary history behind sympatric species in adaptive radiations falls within a continuum from sympatric ecological speciation to secondary contacts of species ecologically specialised in allopatry or parapatry (Rundell and Price 2009, Weber and Strauss 2016). In particular, species complexes can result from adaptive radiations and species segregation along environmental gradients. Species complexes may combine and reshuffle genetic features among species in hybrid swarms (Seehausen 2004, Grant and Grant 2019), thus multiplying the number of potential ecological niches and helping species to achieve reproductive isolation (Runemark et al. 2019). A certain degree of reproductive isolation is necessary to prevent hybrid or derived species from becoming an evolutionary melting pot. Nevertheless, species-specific adaptations can be maintained or even maximised under gene flow, especially with selective pressures varying in space and or time (Tigano and Friesen 2016). For instance, the European white oaks form one of the best known syngameons (Cannon and Petit 2019) with extensive hybridisation (Petit et al. 2002), but each species has a unique ecological niche concerning tolerance to drought, cold, and alkaline soils (Leroy et al. 2019, Cannon and Petit 2019). The coexistence of the different species is partly due to the genes that allow them to survive in different, adjacent, ecological niches (Leroy et al. 2019). To our knowledge, only few studies have reported genetic evidence for differential adaptation to topography or abiotic habitats in tropical tree clades (Fine et al. 2004, Pillon et al. 2014, Paun et al. 2016). At the within-species level though, topography has been highlighted as a driver of genetic divergence (Brousseau et al. 2015).

In the present study, we assessed genetic species delimitation and genotypic diversity of closely-related sympatric tree species belonging to two Neotropical species complexes, the genus Symphonia (Clusiaceae) and the clade Parvifolia of the genus Eschweilera (Lecythidaceae). We addressed fine-scale spatial adaptation along a topographic gradient associated with nutrient and water distribution. Combining LiDAR-derived topographic data, comprehensive tree inventories, and single nucleotide polymorphisms (SNPs) from gene capture experiments, we explored genome-wide population genetic structure, covariation of environmental variables and genotype-environment association to address the following questions:

1. How is genetic diversity structured among and within species of Neotropical tree species complexes ?
2. How are soil water and nutrients distributed at fine spatial scale along the topographic gradient exploited by the species complexes?
3. Are tree species and individuals adapted to the fine-scale topographic gradient?

We hypothesized species to be delineated into gene pools corresponding to already described taxonomic morphotypes. We expected misidentifications due to overlapping intraspecific morphological variability and the difficulty to access reproductive material (Mori et al. 1990). In this context, gene pool delimitation based on molecular data can reveal cryptic species or lead to the merging of taxa previously thought to be distinct (Ewédjè et al. 2020). We hypothesized relative elevation and topographic wetness to explain the distributions of water and nutrients at fine scale (Ferry et al. 2010b, Allié et al. 2015). We expected genetic variation among and within species to be structured along the topographic gradient and to detect genetic signatures of adaptation, following preliminary evidence of topography determining species’ ecological niches [Schmitt et al., in prep]. However, countering this process, interspecific gene flow could reduce neutral and adaptive signatures of differentiation among species within the species complexes (Tigano and Friesen 2016).

Material and Methods

Study site

The study was conducted in the Paracou field station, in the coastal forests of French Guiana, South America. The site is characterized by an average of 3,041 mm annual rainfall and a mean air temperature of 25.71 °C (Aguilos et al. 2018). Old tropical forest with an exceptional richness (i.e. over 750 woody species) grows across the succession of small hills of this area, which rise to 10–40 m a.s.l. (Gourlet-Fleury et al. 2004). The site comprises 16 permanent plots (fifteen 6.25 ha plus one 25 ha) which have been censused (DBH>10) every 1-2 years for more than 35 years. Nine of the plots were logged and subjected to human-induced disturbance in 1986 (details on the experiment in Hérault and Piponiot 2018).

Plant material

Four hundred and two individuals of Symphonia globulifera (Clusiaceae) and 417 individuals belonging to the clade Parvifolia of the genus Eschweilera (Lecythidaceae, Huang et al. 2015, Mori et al. 2016) were sampled in 2017 and 2018 during the dry season (i.e. from September to December) in Paracou (Fig. 7). Both genera are locally abundant and some species within them are Amazonian hyperdominants (Steege et al. 2013). Symphonia globulifera L.f (Clusiaceae) was previously recognized as composed of two morphotypes in French Guiana (Sabatier et al. 1997, Molino and Sabatier 2001, Baraloto et al. 2007). Symphonia globulifera sensu stricto and Symphonia sp.1 occur in sympatry but in differentiated habitats, with S. globulifera preferentially growing in valley bottoms with an acquisitive functional strategy and S. sp1. preferentially exploiting a variety of drier habitats with a conservative functional strategy (Allié et al. 2015, Schmitt et al., in prep; Schmitt et al. 2020). Reciprocal transplantation experiments of Symphonia seedlings have shown that survival and growth performance of each morphotype is better in their home environment than in the opposite environment, showcasing how the two morphotypes are differently adapted to their respective environments (Tysklind et al., 2020). Similarly, Eschweilera sagotiana Miers, E. decolorans Sandwith, and E. coriacea (DC.) S.A.Mori exhibit niche differentiation, differential responses to water stress, and different functional traits along topography (Allié et al. 2015, Schmitt et al., in prep; Baraloto et al. 2007, Schmitt et al. 2020). The genera Symphonia and Eschweilera have been highlighted as species complexes with low (phylo-)genetic species resolution and high levels of plastid DNA sharing among closely related species (Gonzalez et al. 2009, Baraloto et al. 2012a, Huang et al. 2015, Torroba-Balmori et al. 2017, Caron et al. 2019, Heuertz et al. 2020). In addition, outgroups for genetic analysis in Symphonia were comprised of 13 individuals of Symphonia globulifera from Africa (Sao Tome, Gabon, Cameroun, Congo, Benin, Liberia, Ivory Coast, and Ghana), seven Symphonia globulifera from South America (Brazil, Costa Rica and Panama), two Symphonia nectarifera Jum. & H. Perrier from Madagascar, two Symphonia urophylla (Decne. ex Planch. & Triana) Benth. & Hook.f. ex Vesque from Madagascar, five Pentadesma butyracea Sabine from Benin and Cameroon and one Pentadesma grandifolia Baker f. from Cameroon. For Eschweilera, outgroups were selected from other Eschweilera and Lecythis clades in Paracou (Lecythis persistens Sagot, E. simiorum (Benoist) Eyma, and E. chartacea (O.Berg) Eyma). Leaves were collected from the 864 individuals (402 Symphonia + 30 outgroups; 417 Eschweilera + 15 outgroups) and dessicated using silica gel.

Figure 7: Map of sampled individuals for species complexes Symphonia and Eschweilera clade Parvifolia. Colours represent species determined based on the genetic structure of capture datasets, which was not necessarily in agreement with original botanical identification. Grey level represents hill shading to visualize topography.

Sequence capture

Design of probes sets

We designed in silico two sets of 20,000 80-mer probes for sequence capture, one for Symphonia globulifera, and a second one compatible with both Eschweilera sagotiana and Eschweilera coriacea (Fig. 36, Fig. 37).

For Symphonia globulifera, the genomic and transcriptomic resources used for the design were comprised of a published low-coverage draft genome obtained from an individual from Cameroon (Olsson et al. 2017), an unpublished draft genome from an individual from French Guiana [Scotti et al., in prep], an unpublished transcriptome from 20 juveniles from French Guiana [Tysklind et al., in prep], and reduced-representation genomic sequence reads of individuals from French Guiana [Torroba-Balmori et al., unpublished]. Details on the design of the probes set for Symphonia are given in Method S1.

For Eschweilera, we used transcriptomes from Eschweilera sagotiana and Eschweilera coriacea (Vargas et al. 2019), and unpublished reduced representation genomic reads of E. coriacea, E. sagotiana, and E. decolorans (M. Heuertz pers. com.). Details on the design of the probes set for Eschweilera are given in Method S2.

Genomic libraries and sequence capture

Genomic DNA was extracted from 5 mg of dried leaf tissue with a CTAB protocol (Doyle and Doyle 1987). DNA extracts were digested with ‘Ultra II FS Enzyme Mix’ (new England Biolabs Inc, MA, USA) for a target size of 150 bp, and libraries built with the ‘NEBNext Ultra II FS DNA Library Prep kit for Illumina’ (New England Biolabs Inc, MA, USA). We amplified and tagged libraries using 5 \(\mu L\) of adaptor-ligated DNA, 8.3 \(\mu L\) of ‘NEBNext Ultra II Q5 Master Mix’ (new England Biolabs Inc, MA, USA), 2x 1.6 \(\mu L\) of Index Primer i5 and i7 from ‘NEBNext Multiplex Oligos for Illumina (Dual Index Primers Set 1 and Set 2)’ (new England Biolabs Inc, MA, USA). Initial denaturation (98°C for 30 s) was followed by 8 cycles (98°C for 10 s and 65°C for 1 min 30 s) and a final extension (65°C for 5 min). We pooled libraries in four equimolar multiplexes for each genus. We used a custom made set of 20,000 80-mer probes for each genus (myBaits Custom 1-20K, Arbor Biosciences, MI, USA) and conducted the capture experiments using the corresponding myBaits V4 protocol with a hybridization time of 80 hours. We pooled the four multiplexes and sequenced them in two lanes of an Illumina HiSeq 4000 instrument obtaining 2x150bp pair-end reads for each genus.

SNP calling, filtering and annotation in Symphonia

We assessed the quality off raw reads using multiqc (Ewels et al. 2016) and trimmed them with trimmomatic (Bolger et al. 2014). We kept only pair-end reads without adaptors and a phred score above 15 in a sliding window of 4. Seventy percent of trimmed reads mapped off-targets using bwa (Li and Durbin 2009). We thus mapped trimmed reads on the hybrid reference built for the sequence capture experiment using bwa (Li and Durbin 2009), picard (Broad Institute 2018), samtools (Li et al. 2009) and bedtools (Quinlan and Hall 2010). We called variants for each individual using HaplotypeCaller, aggregated variants using GenomicsDBImport and jointly-genotyped individuals using GenotypeGVCFs all in GATK4 software (Auwera et al. 2013). We filtered biallelic SNPs with a quality above 30, a quality by depth above 2, a Fisher strand bias below 60 and a strand odds ratio above 3 using GATK4 (Auwera et al. 2013). Finally, we filtered individuals and SNPs for missing data with a maximum of 95% and 15% of missing data per individual and SNP permitted, respectively, using plink2 (Chen et al. 2019). We obtained 454,262 biallelic SNPs over 385 individuals, without outgroups, for population genetic analysis. Since low-frequency alleles and linkage disequilibrium (LD) will bias the number of fixed loci and increase the number of false-positives in genomic scans for outliers (Foll and Gaggiotti 2008), we built a second dataset for outlier and environmental association analysis, filtering variants with a minor allele frequency (MAF) above 5% (18 individuals) and with LD \(r^2<0.99\). We further removed admixed individuals (<90% gene pool membership, see Genetic species delimitation) and retained 70,737 biallelic SNPs over 372 individuals. We used the genome-transcriptome alignments built for the design of probes sets to annotate called SNPs (Method S1).

SNP calling and filtering in Eschweilera

We conducted quality control and trimming of sequence reads as described for Symphonia. We kept only trimmed reads that mapped on-targets using bwa (Li and Durbin 2009). Because of problems with paralogs in Eschweilera clade Parvifolia due to a strong signature of a past genome duplication (Heuertz et al. 2020), we decided to build a de novo reference of successfully mapping reads using ipyrad with a very strict sequence similarity threshold of 0.95 for clustering within and across individuals (Eaton and Overcast 2020). We automatically tested numerous values of missing data and outgroup filtering, and kept strict filters with population structure congruent with botanical identification but keeping a maximum of individuals (Fig. 38, Fig. 39). We filtered individuals and SNPs for missing data allowing a maximum of 99% of missing data per SNP and at least 5 individuals represented per SNP using plink2 (Chen et al. 2019). Outgroups were automatically filtered out by clustering individuals in two groups using K-means in a genomic principal component analysis (Fig. 40). We obtained 418,793 biallelic SNPs over 257 individuals.

Analyses

Genetic species delimitation

For genetic species delimitation in the Symphonia species complex, we investigated population genetic structure using admixture (Alexander and Lange 2011), using 10 repetitions of K genetic groups with K varying from 1 to 10 and assessed the number of gene pools with cross validation. We defined individuals with a membership to gene pools below 90% as admixed and the remaining individuals as genetically pure. We further investigated admixture with the introgress R package (Gompert and Alex Buerkle 2010), using genetically pure individuals as parental populations and all individuals as the hybrid population. Due to an excess of missing data, admixture failed to infer genetic custers in Eschweilera. We thus used principal component analysis and K-means clustering to delimit Eschweilera species. We assessed the number of gene pools using the within-group sum of squares. We validated gene pool delimitation in Symphonia and Eschweilera by comparison with botanical identifications using a confusion matrix, and we conducted a second blind-identification of every collected individual of Symphonia in November 2019.

We used the genomic scan approach implemented in bayescan (Foll and Gaggiotti 2008) to detect high-differentiation outlier SNPs between the different Symphonia gene pools using the dataset of 70,737 MAF and LD-filtered SNPs with an a priori odds ratio of ten. We considered as outliers those SNPs for which the p-value corrected for multiple testing by false discovery rate was below 5%. We tested genic outlier SNPs for enrichment in gene ontology with the R package clusterProfiler (De La Cruz and Raska 2014).

Topography as a proxy for the distribution of soil water and nutrients

We used two relatively-independent topographic variables, the topographic wetness index (TWI) and the relative elevation (RE, orthogonality shown in Fig. 10), as proxies of the distribution of soil water and nutrients availability in Paracou. Waterlogging and topography have been highlighted as crucial for forest dynamics (Ferry et al. 2010b), species-habitats relationships (Engelbrecht et al. 2007), and phenotypic variation of Symphonia and Eschweilera (Schmitt et al. 2020). TWI and RE were derived from a 1-m resolution digital elevation model using SAGA-GIS (Conrad et al. 2015) based on a LiDAR campaign of the whole Paracou field station done in 2015. Principal component analysis (PCA) was further used to characterize the distribution of variables summarizing soil water and nutrients along TWI and RE in Paracou interpolating previous studies: we specifically assessed the co-variation in organic matter, carbon, nitrogen and phosphorous content, exchangeable cation charge (Soucémarianadin 2004, Roelens in press), water table depth (Gourlet-Fleury et al. 2004), hydromorphy, and waterlogging (Cantet 2004).

Topography effect on neutral and adaptive genetic variation

We did environmental association analyses (Rellstab et al. 2015) in each complex using general linear mixed models developed for genome wide association studies (GWAS). We used topographic wetness index (TWI) or relative elevation (RE) as the response variable and genetic structure (i.e. gene pools representing species) and relatedness (i.e. kinship matrix) as explanatory variables, as it is common practice (Rellstab et al. 2015). This analysis assumed that the topographic conditions where individuals have grown above 10-cm DBH are strongly correlated to the individual heritable phenotypes (e.g. Eckert et al. 2010). We used genetic species and individual kinship in an animal model (Wilson et al. 2010) to estimate genetic variance associated with topography. We used a lognormal likelihood given that distributions of environmental variables were positive and skewed. We defined species based on the gene pools identified with admixture and inferred individual kinship using KING (Manichaikul et al. 2010), as the method is robust to population structure. We set negative kinship values to null as they were confounding with population structure, and we further ensured that the matrix was positive-definite using the nearPD function from the R package Matrix. The environment \(y_{s,i}\) where individual \(i\) in species \(s\) grows was inferred with a lognormal distribution with the following formula:

\[\begin{equation} \begin{aligned} y_{s,i} \sim logN(log(\mu_s.a_{i}),\sigma^2_1) \\ a_{i} \sim MVlogN_N(log(1),\sigma^2_2.K) \end{aligned} \tag{5} \end{equation}\]

where \(\mu_s\) is the mean environment of species \(s\), \(a_i\) is the breeding value of individual \(i\) and \(\sigma^2_1\) is the shape parameter of the lognormal distribution. Individual breeding values \(a_i\) are defined following a multivariate lognormal law \(\mathcal{MVlogN}\) of co-shape matrix defined as the product of the kinship matrix \(K\) with estimated individual genotypic variation \(\sigma^2_2\). To estimate variances on a normal scale, we log-transformed species fixed effects, genetic additive values, and we calculated conditional and marginal \(R^2\) (Nakagawa and Schielzeth 2013). A Bayesian method was used to infer parameters using stan language (Carpenter et al. 2017, code available in Model S1) and rstan package (Stan Development Team 2018) in the R environment (R Core Team 2020) using the No-U-Turn Sampler alogirthm (NUTS, Hoffman and Gelman 2014), which performs better for estimating genetic parameters and breeding values (Nishio and Arakawa 2019).

Results

Symphonia globulifera was previously recognized as composed of two morphotypes in French Guiana, but cross validation of admixture analyses for different levels of clustering advocated for K = 3 different gene pools (Fig. 8a, Fig. 41, Fig. 42). The three gene pools correspond to the two previously identified morphotypes (70-80% of match), S. globulifera and S. sp.1, where the S. globulifera morphotype itself comprises two gene pools. These two gene pools correspond to two locally defined sub-morphotypes herafter called S. globulifera type Paracou and S. globulifera type Régina (Fig. 8b). Symphonia sp.1 grows preferentially in drier plateaux and slopes (Schmitt et al., submitted) and has a light grey thin and smooth bark associated with smaller leaves, flowers, and fruits (Fig. 8b). Symphonia globulifera type Paracou grows in wet habitats of bottomlands, but drier than those of S. globulifera type Régina, and has a dark bark that is intermediate between that of the other morphotypes with respect to thickness and rugosity (Fig. 8b). Symphonia globulifera type Régina grows in the wettest bottomlands of Paracou and has a thick and lashed bark (Fig. 8b). The second blind-identification correctly assigned the difference between all S. globulifera type Paracou and S. globulifera type Régina and reclassified correctly 87% of previous mismatches between gene pools and previous botanical identification. Among the three populations, 25 individuals (7%) showed admixed genotypes, 15 between S. globulifera type Régina and S. sp.1 and 10 between S. globulifera type Paracou and S. sp.1, whereas no admixed individuals were found between S. globulifera type Régina and S. globulifera type Paracou (confirmed with introgress R package, Fig. 43). Admixed individuals were removed from subsequent analyses.

Figure 8: Admixture plot (a) of Symphonia individuals for the best K, K = 3 and corresponding trunk morphology (b). The three morphotypes are identified with their bark with S. sp.1 (left, red) having a light grey thin and smooth bark, the S. globulifera type Paracou (center, light blue) having a dark that is intermediate between that of the other morphotypes with respect to thickness and rugosity and the S. globulifera type Régina that has a thick and lashed bark (right, dark blue).

Genomic scans for Symphonia revealed high-differentiation outliers among gene pools (5.7% of SNPs, Fig. 44), and outliers were enriched in the gene ontology term “response to water deprivation” (\(p=0.01106\)). Based on our results, Symphonia globulifera is locally and genetically structured in three sympatric species, corresponding to three distinct morphotypes (Fig. 8). Eschweilera clade Parvifolia was locally and genetically structured in three genepools separated based on a principal component analysis. The gene pools corresponded to botanical species E. coriacea, E. sagotiana, and E. decolorans (65-80% of match, Fig. 9, Tab. 10). Symphonia and Eschweilera gene pools showed intermediate genetic differentiation, \(F_{ST}=0.15\) and \(F_{ST}=0.35\), respectively.

Figure 9: Principal component analysis (PCA) of single nucleotide polymorphisms (SNP) from Eschweilera clade Parvifolia. The colours represent the clusters detected with Kmeans for the best K, K = 3.

We confirmed topographic wetness index (TWI) and relative elevation (RE, Fig. 11) as descriptors of microenvironmental variation along the topographic gradient exploited by the two species complexes in Paracou. Topographic wetness index (TWI) was related to soil water distribution, opposing water tables deeper than 100 cm to those between the surface and 100 cm deep, which corresponded to hydromorphic soils with the most pronounced waterlogging (variable waterlog5, Fig. 10). Relative elevation (RE) was related to the distribution of nutrients, opposing hilltop soils with the lowest degree of waterlogging (waterlog1) and increased soil organic compounds (soil carbon, nitrogen, organic matter, exchangeable cation charge) to those further down the slopes with intermediate levels of waterlogging and higher soil fertility (higher availability of phosphorus) (Fig. 10).

Figure 10: Principal component analysis (PCA) of topographic wetness index (TWI), relative elevation (RE) with edaphic and hydrological properties. Edaphic and hydrological properties include organic matter (OM), carbon (C), nitrogen (N), exchangeable cation charge (C) (Soucémarianadin 2004, Roelens in press), water table depth (WTD in cm, Gourlet-Fleury et al. 2004), hydromorphy and waterlogging (waterlog in 5 classes, Cantet 2004)). Data were obtained through interpolation of soil sample data on TWI and RE from LiDAR data.

The distribution of Symphonia and Eschweilera species along the topographic gradient suggested that both topographic wetness index (TWI) and relative elevation (RE) drove individual survival across gene pools (\(0.06<R_m^2<0.39\), Fig. 11, Tab. 5). The spatial distribution of Symphonia species was driven by both the topographic wetness index (TWI, Fig. 11c, Tab. 5), especially between S. globulifera type Régina and S. globulifera type Paracou, and the relative elevation (RE, Fig. 11d, Tab. 5), especially between S. sp.1 and the two other species. The spatial distribution of Eschweilera species was more driven by relative elevation (RE, Fig. 11b, Tab. 5) than by topographic wetness index. Variance partitioning revealed that the total genetic variance for topographic position was higher in Symphonia (44 and 46% for TWI and RE, respectively) than in Eschweilera (27 and 25%, respectively) and that the variance component corresponding to species was much higher in Symphonia (35 and 39%) than in Eschweilera (6 and 13%, respectively). Within-species additive genetic variance explained between 7 and 21% of total variance in topographic position in all models of the two species complexes.

Figure 11: Species distribution and individual variance partitioning along topography. Each subfigure shows density of individuals per species (density plot) against the variable considered (topographic wetness index (TWI, left subfigures) or the logarithm of relative elevation (log(RE+1), right subfigures) and the partitioning of individual variation in the positioning along topography (barplot) partitioned with the Animal model into among-species (dark grey), among-genotypes (intermediate grey) and residual (light grey). Species complexes correspond to Eschweilera clade Parvifolia (top subfigures) and Symphonia (bottom subfigures). Colours represent species determined based on the genetic structure of capture datasets, which was not necessarily in agreement with original botanical identification.

Table 5: Parameters and associated standard deviation of the animal model for Symphonia and Eschweilera with topographic wetness index (TWI) and relative elevation (RE). Parameters include species mean (\(\mu_s\)), variance (\(V_{species}\)), genotype variance (\(V_{genotype}\)), and residuals variance (\(V_{residuals}\)). Parameters have been infered with bayesian modelling of the animal model (Model S1) independently for the two species complexes and the two topographic variables.

Species complex	Variable	Parameter	Species	Estimate	Standard deviation
Symphonia	TWI	\(\mu_s\)	S. globulifera Paracou	1.076	0.056
Symphonia	TWI	\(\mu_s\)	S. globulifera Régina	1.302	0.115
Symphonia	TWI	\(\mu_s\)	S. sp.1	0.570	0.021
Symphonia	TWI	\(V_{species}\)		0.111	0.013
Symphonia	TWI	\(V_{genotype}\)		0.038	0.031
Symphonia	TWI	\(V_{residuals}\)		0.171	0.029
Symphonia	RE	\(\mu_s\)	S. globulifera Paracou	1.153	0.038
Symphonia	RE	\(\mu_s\)	S. globulifera Régina	1.109	0.068
Symphonia	RE	\(\mu_s\)	S. sp.1	1.916	0.047
Symphonia	RE	\(V_{species}\)		0.063	0.006
Symphonia	RE	\(V_{genotype}\)		0.014	0.012
Symphonia	RE	\(V_{residuals}\)		0.086	0.012
Eschweilera	TWI	\(\mu_s\)	E. coriacea	0.998	0.059
Eschweilera	TWI	\(\mu_s\)	E. sagotiana	0.775	0.030
Eschweilera	TWI	\(\mu_s\)	E. decolorans	0.772	0.041
Eschweilera	TWI	\(V_{species}\)		0.012	0.006
Eschweilera	TWI	\(V_{genotype}\)		0.055	0.049
Eschweilera	TWI	\(V_{residuals}\)		0.133	0.050
Eschweilera	RE	\(\mu_s\)	E. coriacea	1.376	0.059
Eschweilera	RE	\(\mu_s\)	E. sagotiana	1.742	0.048
Eschweilera	RE	\(\mu_s\)	E. decolorans	1.928	0.072
Eschweilera	RE	\(V_{species}\)		0.015	0.005
Eschweilera	RE	\(V_{genotype}\)		0.021	0.023
Eschweilera	RE	\(V_{residuals}\)		0.077	0.023

Discussion

Despite the local abundance of closely-related tree species growing in sympatry in the Neotropics, little is known of the eco-evolutionary forces that shape their niches and allow their local coexistence. Here, we show that within tree species complexes, closely-related species have different realized optima for topographic niches and we provide evidence for genetic adaptations of species and of genotypes within species for the topographic position where they occur. Symphonia species are adapted to the distribution of water and nutrients captured through the topographic wetness index and the relative elevation, hence they coexist locally through exploiting a broad gradient of local habitats. Conversely, Eschweilera species are differentially adapted to soil chemistry, captured foremost through the relative elevation, and avoid the wettest, hydromorphic habitats. The greater adaptive signature to topographic position in Symphonia than in Eschweilera suggests that topographic niche differentiation may be the main factor explaining sympatric coexistence of Symphonia, while additional factors related to niche and/or fitness may be necessary to explain coexistence of Eschweilera species. Overall, our results suggest that genetic adaptations to different characteristics of topography stabilize local coexistence of species through reducing competition among closely-related species within the two tree species complexes.

Three species of Symphonia adapted to topography coexist in sympatry

We identified three distinct species within Symphonia globulifera sensu lato with different morphologies, distinct gene pools, and topographic niches. For the purpose of the study, we named the three species S. sp.1, S. globulifera type Paracou, and S. globulifera type Régina in reference to already existing morphotypes and local names. The three species showed a low but significant genetic differentiation (\(F_{ST} = 0.15\)), as shown previously at the local scale between the two previously recognized morphotypes, S. globulifera sensu stricto and S. sp.1 (\(F_{ST} = 0.114\) for genic SSRs; Olsson et al. 2017). Evidence for adaptive genetic differentiation among species was revealed by 5.7% of SNPs identified as high-differentiation outliers, with an enrichment in genes implicated in response to water deprivation.. Topography has already shown to be a driver of the distribution of the two previously recognized morphotypes (Allié et al. 2015) with S. globulifera growing in low-elevation and wet bottomlands and S. sp.1 growing in drier slopes and plateaux. In addition, our study revealed the existence of two Symphonia globulifera morphotypes segregated within wet areas, with Symphonia globulifera type Régina growing in the wettest bottomlands, such as swamp characterized by increased wetness and anoxia. Similar ecotypic differentiation with a designated swamp ecotype has been suggested in Africa (Budde et al. 2013). Individual topographic positions had a large plastic component (54-56% of variation) but genetics explained the remaining variation, with a ca. five-fold larger inter- than intra-specific component (Fig. 11). Thus, finescale topographic variation among individuals was driven by adaptive processes, revealing local adaptation to topography within and among the three species at the hectare-scale. Further studies at a broader regional scale are needed, but our results advocate for the definition of three distinct taxonomic species with differences in morphology (Fig. 8), habitat preference (Fig. 11), functional traits and growth trajectories (S. Schmitt unpublished results).

Three species of Eschweilera adapted to soil chemistry coexist in sympatry

Eschweilera species avoided hydromorphic soils and were differentially adapted to soil chemistry. At least three distinct species of Eschweilera clade Parvifolia with different morphologies, gene pools, and topographic niches coexist in sympatry. Eschweilera coriacea preferentially grows in fertile and phosphorus-rich valley bottoms but has a broad topographic niche, and E. sagotiana and E. decolorans preferentially exploit slopes and plateaus, characterized by low levels of waterlogging and soils rich in organic compounds (Fig. 11; Allié et al. 2015). Topography, through the dissolution of iron oxides, litter- and tree-fall transfers and waterlogging, shape soil nutrient distribution in tropical forests (John et al. 2007, Ferry et al. 2010a). The lack of power of Eschweilera data prevented us from detecting the less abundant Eschweilera species of the Parvifolia clade (e.g. E. grandiflora (Aubl.) Sandwith, E. wachenheimii (Benoist) Sandwith), but our results do not question their botanical definitions. Topographic position within Eschweilera was mainly plastic (73-75% of residual variation), but concerning genetic variation, relative elevation best captured differences among species (52%) whereas the topographic wetness index best captured genetic variation among genotypes (78%, Fig. 11).

Species niche differentiation and coexistence within species complexes

The interaction between topography and genes shapes the microgeographic genetic structure among and within species in the studied tropical tree species complexes. We found genomic evidence for a functional response to water deprivation between Symphonia species, confirming the role of water availability as a driver of the topographic niches in the complex. Similarly to European white oaks, Symphonia holds coexisting species with differential genetic adaptations that allow them to thrive in different ecological niches despite interspecific hybridization (Leroy et al. 2019). In Eschweilera, we evidenced a genomic signature of topographic niche differentiation but the low resolution of our data did not allow us to demonstrate hybridization. Although the regional-scale evolutionary history and the specific roles of interspecific hybridization and introgression remain unknown for both species complexes, our study shows that inter- and intra-specific adaptive processes to topographic niches stabilize local coexistence of species in these tropical tree species complexes, similarly to what is observed in the European white oaks syngameon. Symphonia and Eschweilera clade Parvifolia might thus represent syngamons (Cannon and Petit 2019) of tropical trees. It remains to be shown that interspecific gene flow benefits their constituent species, through adaptive introgression, or through preventing local extinction of congeners by enabling reproductive assurance and increasing their effective population size (Cannon and Lerdau 2015).

The distribution of Eschweilera species along the topographic gradient was more influenced by plasticity than in Symphonia and the role of intra-specific adaptive processes was stronger in Eschweilera, with intra-specific genotypic effects stronger or equal to inter-specific effects. As compared to Symphonia, Eschweilera exhibited globally a weaker role of topography on the structuring of its genetic variation, its species had a more marked topographic niche overlap, and the intra-specific adaptive signature to topography was larger. These patterns may be a consequence of the evolutionary history of the complex. Eschweilera species, despite similar morphology, have a higher genetic differentiation and their phylogenetic constraints may thus be relaxed as compared to Symphonia (Wiens et al. 2010), allowing greater topographic niche overlap possibly in conjunction with niche differentiation along other ecological gradients, e.g., related to root depth or light access (Rüger et al. 2009). But Eschweilera species may also coexist through a decrease in fitness differences among species allowing them to persist in a shared niche (Turcotte and Levine 2016). This explanation is plausible in a species complex like Eschweilera in which the risk of hybridization and consequent breakdown of inter-specific divergence may be more limited than in more closely-related, recently diverged species (Tobias et al. 2014). Furthermore, the signature of within-species genotypic adaptations to topographic niche may be associated with broad niches for TWI. Schmitt et al. (2020) showed that E. coriacea and the two major previously recognized Symphonia morphotypes had broader niches for TWI than the majority of tree species in Paracou. To further improve our understanding on the eco-evolutionary drivers of coexistence in these tropical tree species complexes, future work should improve our knowledge on the evolutionary history of both complexes, including aspects of interspecific gene flow and introgression, as well as consolidate knowledge on the factors that shape the niche space of both species complexes at broader geographic scale.

The literature includes numerous examples of niche differentiation among closely-related species growing in sympatry, along fine-scale topography (Gunatilleke et al. 2006, Engelbrecht et al. 2007, Kraft et al. 2008, Allié et al. 2015) or along other local ecological niche variables (Itoh et al. 2003, Yamasaki et al. 2013). On a larger scale, topography is also a driving force for tropical species diversity on the Andes (Mutke et al. 2014). In Paracou, several species complexes have been shown to present differentiated niches along topography [Schmitt et al. in prep]. Combined with our results, this abundant evidence of niche differentiation among closely related species calls into question the more general role of topography in the genetic adaptation and diversification of tropical tree species complexes, even for small topographic gradients such as those found in French Guiana. Our results suggest that ecological speciation as demonstrated along steep gradients, such as soil types or soil metal content (Fine et al. 2004, Paun et al. 2016), could possibly occur at a fine spatial scale with gradual habitat variation. Our work emphasizes the potential role of species complexes and syngameons in the origin and maintenance of Neotropical forest diversity.

Acknowledgements

We thank the University of Bordeaux for a PhD grant to Sylvain Schmitt. We are grateful to Pascal Petronelli and the CIRAD inventory team for their work on tree inventories and botanical identification. This study was funded through an Investissement d’Avenir grant of the ANR: CEBA (ANR-10-LABEX-0025), the “projet innovant” LOCOCAP granted through INRA and the LECYTOMICS project granted through the IRD-GUYAMAZON program.

Authors’ contributions

All authors conceived the ideas and designed methodology; SS, MH, and NT sampled individuals; SS and MH performed the gene capture experiments; SS, MH and BH analysed outputs; SS and MH led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Supporting information

Method S1. Design of the probes set for Symphonia.

Method S2. Design of the probes set for Eschweilera.

Model S1. Stan code for the animal model

Tab. 10. Eschweilera botanical species and genetic clusters

Fig. 36. Target selection for the capture experiment of Symphonia

Fig. 37. Target selection for the capture experiment of Eschweilera

Fig. 38. SNP abundance per library for Eschweilera

Fig. 39. Library abundance per SNP for Eschweilera

Fig. 40. Outgroup detection for Eschweilera

Fig. 41. Cross-validation for Symphonia population structure

Fig. 42. Symphonia population structure using admixture

Fig. 43. Symphonia population structure using hybrid index

Fig. 44. Outlier SNPs among Symphonia species

References

Abbott, R. et al. 2013. Hybridization and speciation. - Journal of Evolutionary Biology 26: 229–246.

Aguilos, M. et al. 2018. Interannual and Seasonal Variations in Ecosystem Transpiration and Water Use Efficiency in a Tropical Rainforest. - Forests 10: 14.

Alexander, D. H. and Lange, K. 2011. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. - BMC Bioinformatics in press.

Allié, E. et al. 2015. Pervasive local-scale tree-soil habitat association in a tropical forest community. - PLoS ONE 10: e0141488.

Auwera, G. A. et al. 2013. From FastQ Data to High-Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline. - Current Protocols in Bioinformatics 43: 483–492.

Baraloto, C. et al. 2007. Seasonal water stress tolerance and habitat associations within four Neotropical tree genera. - Ecology 88: 478–489.

Baraloto, C. et al. 2012a. Using functional traits and phylogenetic trees to examine the assembly of tropical tree communities. - Journal of Ecology 100: 690–701.

Bolger, A. M. et al. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. - Bioinformatics 30: 2114–2120.

Broad Institute 2018. Picard Tools.

Brousseau, L. et al. 2015. Neutral and Adaptive Drivers of Microgeographic Genetic Divergence within Continuous Populations: The Case of the Neotropical Tree Eperua falcata (Aubl.) (FA Aravanopoulos, Ed.). - PLOS ONE 10: e0121394.

Budde, K. B. et al. 2013. The ancient tropical rainforest tree Symphonia globulifera L. f. (Clusiaceae) was not restricted to postulated Pleistocene refugia in Atlantic Equatorial Africa. - Heredity 111: 66–76.

Cannon, C. H. and Lerdau, M. 2015. Variable mating behaviors and the maintenance of tropical biodiversity. - Frontiers in Genetics 6: 183.

Cannon, C. H. and Petit, R. J. 2019. The oak syngameon: more than the sum of its parts. - New Phytologist: nph.16091.

Cantet, L. 2004. Prédiction de l’engorgement hydrique de surface par les cortèges floristiques en forêt tropicale humide (Guyanefrançaise).

Caron, H. et al. 2019. Chloroplast DNA variation in a hyperdiverse tropical tree community. - Ecology and Evolution 9: ece3.5096.

Carpenter, B. et al. 2017. Stan : A Probabilistic Programming Language. - Journal of Statistical Software in press.

Chen, Z. L. et al. 2019. A high-speed search engine pLink 2 with systematic evaluation for proteome-scale identification of cross-linked peptides. - Nature Communications in press.

Chesson, P. 2000a. Mechanisms of maintenance of species diversity. - Annual review of Ecology and Systematics 31: 343–366.

Conrad, O. et al. 2015. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. - Geoscientific Model Development 8: 1991–2007.

De La Cruz, O. and Raska, P. 2014. Population structure at different minor allele frequency levels. - BMC proceedings 8: S55.

Doyle, J. and Doyle, J. 1987. Genomic plant DNA preparation from fresh tissue-CTAB method. - Phytochem Bull 19: 11–15.

Eaton, D. A. and Overcast, I. 2020. ipyrad: Interactive assembly and analysis of RADseq datasets. - Bioinformatics (Oxford, England) 36: 2592–2594.

Eckert, A. J. et al. 2010. Patterns of population structure and environmental associations to aridity across the range of loblolly pine (Pinus taeda L., Pinaceae). - Genetics 185: 969–982.

Engelbrecht, B. M. et al. 2007. Drought sensitivity shapes species distribution patterns in tropical forests. - Nature 447: 80–82.

Ewels, P. et al. 2016. MultiQC: Summarize analysis results for multiple tools and samples in a single report. - Bioinformatics 32: 3047–3048.

Ewédjè, E. E. B. K. et al. 2020. Species delimitation in the African tree genus Lophira (Ochnaceae) reveals cryptic genetic variation. - Conservation Genetics 21: 501–514.

Ferry, B. et al. 2010a. Higher treefall rates on slopes and waterlogged soils result in lower stand biomass and productivity in a tropical rain forest. - Journal of Ecology 98: 106–116.

Ferry, B. et al. 2010b. Higher treefall rates on slopes and waterlogged soils result in lower stand biomass and productivity in a tropical rain forest. - Journal of Ecology 98: 106–116.

Fine, P. V. A. et al. 2004. Herbivores promote habitat specialization by trees in Amazonian forests. - Science (New York, N.Y.) 305: 663–5.

Foll, M. and Gaggiotti, O. 2008. A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective. - Genetics 180: 977–993.

Gentry, A. H. 1988. Tree species richness of upper Amazonian forests. - Proceedings of the National Academy of Sciences 85: 156–159.

Gompert, Z. and Alex Buerkle, C. 2010. Introgress: A software package for mapping components of isolation in hybrids. - Molecular Ecology Resources 10: 378–384.

Gonzalez, M. A. et al. 2009. Identification of amazonian trees with DNA barcodes. - PLoS ONE 4: e7483.

Gourlet-Fleury, S. et al. 2004. Ecology and management of a neotropical rainforest : lessons drawn from Paracou, a long-term experimental research site in French Guiana Ecology and management of a neotropical rainforest : lessons drawn from Paracou, a long-term experimental research sit. - Paris: Elsevier.

Grant, P. R. and Grant, B. R. 2019. Hybridization increases population variation during adaptive radiation. - Proceedings of the National Academy of Sciences of the United States of America 116: 23216–23224.

Gunatilleke, C. V. S. et al. 2006. Species–habitat associations in a {Sri} {Lankan} dipterocarp forest. - Journal of Tropical Ecology 22: 371.

Heuertz, M. et al. 2020. The hyperdominant tropical tree Eschweilera coriacea (Parvifolia clade, Lecythidaceae) shows higher genetic heterogeneity than sympatric Eschweilera species in French Guiana. - Plant Ecology and Evolution. 153: 67–81.

Hérault, B. and Piponiot, C. 2018. Key drivers of ecosystem recovery after disturbance in a neotropical forest. - Forest Ecosystems 5: 2.

Hoffman, M. D. and Gelman, A. 2014. The no-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. - Journal of Machine Learning Research 15: 1593–1623.

Huang, Y. Y. et al. 2015. Toward a phylogenetic-based generic classification of neotropical lecythidaceae—I. Status of Bertholletia, Corythophora, Eschweilera and Lecythis. - Phytotaxa 203: 085–121.

Itoh, A. et al. 2003. Importance of topography and soil texture in the spatial distribution of two sympatric dipterocarp trees in a Bornean rainforest. - Ecological Research 18: 307–320.

John, R. et al. 2007. Soil nutrients influence spatial distributions of tropical tree species. - Proceedings of the National Academy of Sciences of the United States of America 104: 864–9.

Kraft, N. J. B. et al. 2008. Functional traits and niche-based tree community assembly in an Amazonian forest. - Science 322: 580–582.

Leroy, T. et al. 2019. Massive postglacial gene flow between European white oaks uncovered genes underlying species barriers. - New Phytologist in press.

Levin, D. A. et al. 1996. Hybridization and the extinction of rare plant species. - Conservation Biology 10: 10–16.

Li, H. and Durbin, R. 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. - Bioinformatics 25: 1754–1760.

Li, H. et al. 2009. The Sequence Alignment/Map format and SAMtools. - Bioinformatics 25: 2078–2079.

Lortie, C. J. et al. 2004b. Rethinking plant community theory. - Oikos 107: 433–438.

Manichaikul, A. et al. 2010. Robust relationship inference in genome-wide association studies. - Bioinformatics 26: 2867–2873.

Molino, J.-F. and Sabatier, D. 2001. Tree Diversity in Tropical Rain Forests: A Validation of the Intermediate Disturbance Hypothesis. - Science 294: 1702–1704.

Mori, S. A. et al. 1990. Lecythidaceae, Part 2. The Zygomorphic-Flowered New World Genera (Couroupita, Coroythophora, Bertholletia, Couratari, Eschweilera, & Lecythis), With a Study of Secondary Xylem of Neotropical Lecythidacea. - Flora Neotropica 21: 1–373.

Mori, S. A. et al. 2016. Observations on the phytogeography of the LECYTHIDACEAE clade (Brazil nut family). 30: 1–85.

Mutke, J. et al. 2014. Diversity patterns of selected Andean plant groups correspond to topography and habitat dynamics, not orogeny. - Frontiers in Genetics in press.

Nakagawa, S. and Schielzeth, H. 2013. A general and simple method for obtaining R2 from generalized linear mixed-effects models. - Methods in Ecology and Evolution 4: 133–142.

Nishio, M. and Arakawa, A. 2019. Performance of Hamiltonian Monte Carlo and No-U-Turn Sampler for estimating genetic parameters and breeding values. - Genetics Selection Evolution 51: 1–12.

Olsson, S. et al. 2017. Development of genomic tools in a widespread tropical tree, Symphonia globulifera L.f.: a new low-coverage draft genome, SNP and SSR markers. - Molecular Ecology Resources 17: 614–630.

Paun, O. et al. 2016. Processes Driving the Adaptive Radiation of a Tropical Tree (Diospyros, Ebenaceae) in New Caledonia, a Biodiversity Hotspot. - Systematic Biology 65: 212–227.

Pernès, J. and Lourd, M. 1984. Organisation des complexes d’espèces. - Gestion de ressources genetiques des plantes in press.

Petit, R. J. et al. 2002. Identification of refugia and post-glacial colonisation routes of European white oaks based on chloroplast DNA and fossil pollen evidence. - Forest Ecology and Management 156: 49–74.

Pillon, Y. et al. 2014. Cryptic adaptive radiation in tropical forest trees in New Caledonia. - New Phytologist 202: 521–530.

Pinheiro, F. et al. 2018. Plant Species Complexes as Models to Understand Speciation and Evolution: A Review of South American Studies. - Critical Reviews in Plant Sciences 37: 54–80.

Quinlan, A. R. and 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.

Rellstab, C. et al. 2015. A practical guide to environmental association analysis in landscape genomics. 24: 4348–4370.

Roelens, J. Cartographie p{édologique de six parcelles exploitées.: 25.

Rundell, R. J. and Price, T. D. 2009. Adaptive radiation, nonadaptive radiation, ecological speciation and nonecological speciation. - Trends in Ecology & Evolution 24: 394–399.

Runemark, A. et al. 2019. Eukaryote hybrid genomes. - PLOS Genetics 15: e1008404.

Rüger, N. et al. 2009. Response of recruitment to light availability across a tropical lowland rain forest community. - Journal of Ecology 97: 1360–1368.

Sabatier, D. et al. 1997. The influence of soil cover organization on the floristic and structural heterogeneity of a Guianan rain forest. - Plant Ecology 131: 81–108.

Scheffer, M. and Van Nes, E. H. 2006. Self-organized similarity, the evolutionary emergence of groups of similar species. - Proceedings of the National Academy of Sciences of the United States of America 103: 6230–6235.

Schmitt, S. et al. 2020. Topography consistently drives intra- and inter-specific leaf trait variation within tree species complexes in a Neotropical forest. - Oikos: oik.07488.

Seehausen, O. 2004. Hybridization and adaptive radiation. - Trends in Ecology & Evolution 19: 198–207.

Soucémarianadin, L. 2004. Recherche de critères du sol influançant la structure et la composition floristique d’une forêt tropicale humide (Paracou, Guyane française).

Stan Development Team 2018. Rstan: the R interface to Stan.

Steege, H. ter et al. 2013. Hyperdominance in the Amazonian Tree Flora. - Science 342: 1243092–1243092.

Suarez-Gonzalez, A. et al. 2018. Adaptive introgression: a plant perspective. - Biology Letters 14: 20170688.

Taylor, E. B. et al. 2006. Speciation in reverse: Morphological and genetic evidence of the collapse of a three-spined stickleback (Gasterosteus aculeatus) species pair. - Molecular Ecology 15: 343–355.

Tigano, A. and Friesen, V. L. 2016. Genomics of local adaptation with gene flow. - Molecular Ecology 25: 2144–2164.

Tobias, J. A. et al. 2014. Species coexistence and the dynamics of phenotypic evolution in adaptive radiation. - Nature 506: 359–363.

Torroba-Balmori, P. et al. 2017. Altitudinal gradients, biogeographic history and microhabitat adaptation affect fine-scale spatial genetic structure in African and Neotropical populations of an ancient tropical tree species. - PLOS ONE 12: e0182515.

Turcotte, M. M. and Levine, J. M. 2016. Phenotypic Plasticity and Species Coexistence. - Trends in Ecology & Evolution 31: 803–813.

Vargas, O. M. et al. 2019. Target sequence capture in the Brazil nut family (Lecythidaceae): Marker selection and in silico capture from genome skimming data. - Molecular Phylogenetics and Evolution 135: 98–104.

Weber, M. G. and Strauss, S. Y. 2016. Coexistence in Close Relatives: Beyond Competition and Reproductive Isolation in Sister Taxa. - Annual Review of Ecology, Evolution, and Systematics 47: 359–381.

Weiher, E. and Keddy, P. A. 1995. Assembly rules, null models, and trait dispersion: new questions from old patterns. - Oikos 74: 159–164.

Whitney, K. D. et al. 2010. Patterns of hybridization in plants. - Perspectives in Plant Ecology, Evolution and Systematics 12: 175–182.

Wiens, J. J. et al. 2010. Niche conservatism as an emerging principle in ecology and conservation biology. - Ecology Letters 13: 1310–1324.

Wilson, A. J. et al. 2010. An ecologist’s guide to the animal model. - Journal of Animal Ecology 79: 13–26.

Yamasaki, N. et al. 2013. Coexistence of two congeneric tree species of Lauraceae in a secondary warm-temperate forest on Miyajima Island, south-western Japan. - Plant Species Biology 28: 41–50.