A target capture approach for phylogenomic analyses at multiple evolutionary timescales in rosewoods (Dalbergia spp.) and the legume family (Fabaceae)

Understanding the genetic changes associated with the evolution of biological diversity is of fundamental interest to molecular ecologists. The assessment of genetic variation at hundreds or thousands of unlinked genetic loci forms a sound basis to address questions ranging from micro-to macro-evolutionary timescales, and is now possible thanks to advances in sequencing technology. Major difficulties are associated with i) the lack of genomic resources for many taxa, especially from tropical biodiversity hotspots, ii) scaling the numbers of individuals analyzed and loci sequenced, and iii) building tools for reproducible bioinformatic analyses of such datasets. To address these challenges, we developed a set of target capture probes for phylogenomic studies of the highly diverse, pantropically distributed and economically significant rosewoods (Dalbergia spp.), explored the performance of an overlapping probe set for target capture across the legume family (Fabaceae), and built a general-purpose bioinformatics pipeline. Phylogenomic analyses of Dalbergia species from Madagascar yielded highly resolved and well supported hypotheses of evolutionary relationships. Population genomic analyses identified differences between closely related species and revealed the existence of a potentially new species, suggesting that the diversity of Malagasy Dalbergia species has been underestimated. Analyses at the family level corroborated previous findings by the recovery of monophyletic subfamilies and many well-known clades, as well as high levels of gene tree discordance, especially near the root of the family. The new genomic and bioinformatics resources will hopefully advance systematics and ecological genetics research in legumes, and promote conservation of the highly diverse and endangered Dalbergia rosewoods.

at both ends until an alignment site shows nucleotides across a specified minimum fraction of 292 aligned sequences, along with a specified maximum nucleotide diversity (i.e., the mean 293 number of base differences between all sequence pairs). In addition, internal trimming is 294 performed by only keeping sites with nucleotides in a specified minimum fraction of aligned 295 sequences. Potential mis-assemblies or mis-alignments at contig ends are further resolved 296 using a sliding window approach that identifies and masks sequences with large deviations 297 To mitigate potential biases arising from the reference sequences used, a new set of target 319 region reference sequences can be generated based on the target region alignments 320 generated in the two previous steps. For this purpose, a consensus sequence is generated for 321 each alignment as in step 6, but separate consensus sequences can be generated for different 322 specified taxon groups (see step 4). These sets of taxon group specific consensus sequences 323 are then aligned, and representative consensus sequences are generated as in step 6 (see 324 Supplementary Methods). These taxon-specific reference sequences are the main output of 325 step 7 and can be used to refine mapping, assembly and alignment by reiterating steps 1 to 5. 326 327 2.4.8 | Alignment assessment and filtering 328 We characterized all non-overlapping trimmed alignments for the number of gaps, gap ratio 329 (i.e, the fraction of non-nucleotides in the alignment), total nucleotide diversity, average 330 nucleotide diversity per site, and alignment length, as well as the number and proportion of 331 segregating and parsimony informative sites. We then visually filtered alignments using 332 filter.visual.alignments.R, which allows to apply filtering thresholds that are informed by 333 visualizations of alignment statistics (see Supplementary Methods). We used the filtered 334 alignments after the second iteration of step 6 for phylogenetic analyses. 335 336

| Phylogenetic analyses 337
We performed phylogenetic analyses with both the subfamily and species sets, using a 338 supermatrix (concatenation) approach and a gene tree summary approach. For the 339 supermatrix approach, we ran maximum likelihood searches on the concatenated alignments 340 using RAXML version 8.2.11 (Stamatakis, 2014) with rapid bootstrap analysis and search for 341 the best-scoring tree in the same run (-f a option), 100 bootstrap replicates, and the GTRCAT 342 Detarioideae) as ingroups; Figure S2), using the tree scoring option in ASTRAL-III and a file with 352 the assignment of taxa to subfamilies or the outgroup. All phylogenetic trees were displayed 353 using GGTREE version 2.0.2 (Yu et al., 2016). 354 355

| Population genetic analyses 356
We carried out population genetic analyses for the population dataset. We mapped quality-357 filtered reads against the target region reference sequences that were representative of the 358 species set after the second iteration using BWA-MEM. We verified efficient recovery of target 359 regions by plotting heatmaps of coverage statistics, removed PCR duplicates using PICARD TOOLS 360 version 2.21.3 (Broad Institute, 2019), and capped excessive coverage to 500 using 361 biostar154220.jar (Lindenbaum, 2015). We then called variants using FREEBAYES version 1.1.0- covariance matrix, as well as to calculate a neighbor-joining (NJ) tree (Saitou & Nei, 1987) on 372 Nei's genetic distances, as implemented in POPPR version 2.8.1 (Kamvar et al., 2014). We also 373 used the allele table to create a SNP subset for population clustering analysis using STRUCTURE 374 version 2.3.4 (Pritchard et al., 2000). Specifically, we kept SNPs with genotype data in at least 375 95% of individuals, and we randomly sampled up to three SNPs per target region for linkage 376 disequilibrium pruning and computational ease. STRUCTURE analyses were performed for one 377 to ten demes (K), using 110,000 iterations, including a burn-in period of 10,000 iterations, with 378

| Two probe sets for target capture across legumes and Dalbergia 383
We obtained 0.13 to 13.76 (median: 1.56) million raw read pairs per sample, of which we 384 retained 86.55% to 99.34% (median: 93.82%) after quality trimming (Tables S1 -S3). In the 385 first iteration applied to twelve representative samples, reads mapped to 6,519 or 6,287 of 386 the 6,555 initial target regions in the subfamily or species set, respectively (step 1). Of these 387 we retained 3,436 or 4,908 target regions, which showed adequate coverage across taxon 388 groups. After assembly (step 2) and orthology assessment (step 3), 2,710 or 4,181 target 389 regions passed the region specificity and sensitivity filters of lower stringency (step 4). 390 Following alignment and trimming (step 5), overlapping portions in 207 or 377 regions were 391 successfully merged, resulting in 2,468 or 3,736 non-overlapping trimmed alignments (step 6). 392 20 Longer and more representative consensus sequences were generated from these target 393 regions (step 7) and used as references for mapping quality-trimmed reads of the complete 394 taxon sets (step 1, see Tables S1 and S2). We retained 1,917 or 3,418 target regions with 395 adequate coverage ( Figures S5 and S6), of which 1,020 or 2,407 passed the specificity and 396 Other relationships among subfamilies remained unresolved using both phylogenetic 434 methods ( Figure 1, Figure S1). In particular, a clade comprising Caesalpinioideae, 435 Cercidoideae, Detarioideae and Dialioideae as sister group to Papilionoideae was not 436 supported in the supermatrix tree, and was recovered in only 47% of quartet trees. We Madagascar (location 5, see Figure S4 and Table S3)     In this study, we enriched DNA libraries from three taxon sets spanning micro-529 evolutionary (populations) to macro-evolutionary (family) timescales, using a single set of 530 26 12,049 RNA probes targeting 6,555 genomic regions conserved across five Meso-531 Papilionoideae genomes and a Dalbergia transcriptome. We then identified 2,396 and 1,005 532 target regions with high capture specificity and sensitivity within the species-rich genus 533 Dalbergia (Dalbergia2396 probe set) and more broadly across legumes (Fabaceae1005 probe 534 set). We used our CAPTUREAL pipeline to refine phylogenomic and population genomic analyses 535 using taxon-specific and longer reference sequences. This procedure has both benefits and 536 drawbacks. An advantage is that different but overlapping probe sets amenable for efficient 537 target capture in different focal groups can be identified, and that a single enriched DNA 538 library can be included in multiple data sets spanning different evolutionary timescales. On 539 the other hand, bioinformatic analyses took longer due to the iterative refinement, and only 540 a portion of captured sequence data was ultimately used for phylogenomic or population 541 genomic analyses in each focal group (see Tables S1 -S3). Higher costs per used sequence

| CAPTUREAL bioinformatics pipeline 547
The CAPTUREAL pipeline starts with the mapping of quality-trimmed reads to conserved target 548 regions identified during probe design, followed by assembly on a per-region basis, orthology 549 assessment, and filtering for target regions with high capture sensitivity and specificity for 550 downstream analyses. This approach differs from the PHYLUCE pipeline (Faircloth, 2016), 551 where quality-trimmed reads are first assembled, and then matched to target regions. 552 CAPTUREAL simplifies the assembly of reads specific to each locus, circumventing the 553 27 challenging task of de novo assembly of contigs from the large pool of sequencing reads 554 representative of thousands of loci (reviewed by Chaisson et al., 2015). Likewise, alignments 555 are conducted in clearly defined target regions in which overlap among individual contigs is 556 higher. However, assembly per region is more time-consuming and requires reference 557 sequences for the initial mapping step. This might introduce a reference bias when divergent 558 sequences are not mapped (Lunter & Goodson, 2011). We addressed this problem by 559 generating consensus sequences that are representative of a given taxon set and by limiting 560 analyses to target regions that can be efficiently recovered in all groups of that taxon set. 561 These set-specific reference sequences can then be used to iteratively refine mapping, 562 assembly, and target region filtering for any set of taxa. Our approach is conceptually similar 563 to the HYBPIPER pipeline (Johnson et al., 2016), which also employs a mapping-assembly 564 strategy, and uses depth of coverage and percent identity to the target region to choose 565 between multiple contigs, before it identifies intron/exon boundaries using target peptide 566 regions. It is therefore suitable for application in systems lacking high-quality annotated 572 reference genomes or transcriptomes. The main strengths of this pipeline are its modularity, 573 which allows for an iterative refinement of read mapping, assembly and alignment, its 574 flexibility given by user-defined parameters, the merging of alignments representing 575 28 physically overlapping target regions, and the visualization of key summary statistics and 576 alignments along the workflow to inform the user on critical analysis parameters. 577 578

| Macro-and micro-evolutionary patters in Dalbergia 579
Dalbergia species endemic to Madagascar were recovered as two large, well-supported and 580 fully resolved clades, each exclusively comprising Malagasy species. These two clades were 581 previously identified on the basis of three chloroplast markers, but phylogenetic relationships 582 within clades were not resolved, which exposed traditional DNA barcoding as insufficient for while the multispecies approach applied in this study is expected to return a phylogeny that 624 reflects nuclear evolution accounting for incomplete lineage sorting, conflicts in gene tree 625 topologies due to hybridization and chloroplast capture can further underlie the observed 626

differences. 627
Our target capture approach demonstrated great potential to facilitate the resolution 628 of several taxonomic conundrums within the genus, which likely resulted from few observable 629 and diagnostic morphological characters, insufficient collection effort, and the difficulty of 630 distinguishing between heritable and plastic trait variation within and among Dalbergia 631 species (Lachenaud, 2016). The integration of highly informative museum specimens,

| Phylogenetic analyses across legumes 654
At the family level, 1,005 merged regions of the 6,555 targeted regions passed our stringent 655 sensitivity and specificity filters, suggesting that many target regions were not efficiently 656 captured across taxa. However, phylogenetic analysis of 986 nuclear target regions recovered 657 multiple known clades within monophyletic subfamilies with strong bootstrap and quartet 658 support, providing excellent resolution comparable to that obtained in the recent nuclear 659 phylogenomic analysis of transcriptome and genome-wide data across legumes (Koenen et 660 al., 2020b). As in that study, we found high support for Cercidoideae and Detarioideae as sister 661 taxa, a relationship that was never inferred in analyses based on chloroplast genes (LPWG, Substantial gene tree incongruence was also found with respect to the relationships 682 among the three large clades within Meso-Papilionoidae. The sister relationship between 683 Dalbergioids s.l. and Genistoids s.l. received only slightly higher quartet support than the two 684 alternative hypotheses, which is consistent with previous results (Koenen et al., 2020b).

| Conclusions and perspectives 692
The resources developed here for Fabaceae and in particular the genus Dalbergia bridge 693 micro-and macro-evolutionary timescales and will hopefully facilitate community-driven 694 efforts to advance legume genomics. Comprehensive sampling and sequencing by target 695 capture of Dalbergia across its distribution range, and in particular from the hotspot of 696 diversity in Madagascar, can yield valuable insights into the origin and diversification of the 697 genus, thereby informing conservation policies and the taxonomic revision of Malagasy 698 Dalbergia. The obtained sequence data will further serve to build a reference library for 699 molecular identification of CITES-listed Dalbergia species, which would make a significant 700 contribution toward the conservation of the valuable and endangered rosewoods.