Genome alignment-based probe design improves phylogenomic resolution in target enrichment datasets of ultra-diverse taxa

Abstract Target enrichment methods have provided unprecedented advances in phylogenomics. Targeting hundreds of conserved regions has proven to be a good tradeoff between cost and efficiency, while being useful for museomics and diversified non-model clades. Unfortunately, current methods used for identifying such regions involve high degrees of conservation within targeted elements, usually pushing researchers to rely on flanking data with little guarantee for homology. With a growing number of high quality genomes available throughout the Tree of Life emerges new opportunities to improve marker selection. In this study, we introduce GABBI, a new method for designing target capture probes by taking advantage of genome alignments, avoiding the selection of a single reference genome that can cause notable biases. We compare GABBI-derived markers to the most commonly used probe design method, PHYLUCE, at two taxonomic scales, the weevil superfamily Curculionoidea and the tribe Pachyrhynchini. At both taxonomic scales, results show that our new method allows identifying more variable loci that prove to be more phylogenetically resolutive than the PHYLUCE-derived ones. Doing so, we provide the first probe set specifically designed for weevils, targeting a wide set of 4,255 shared homologous regions, encouraging future research on systematics and macroevolution of one of the most diverse and economically important groups of insects. By providing GABBI as an automated and open-access pipeline, we hope to open new probe design opportunities to other taxonomic groups that face similar phylogenetic obstacles. keywords: target capture, genome alignment, weevil, UCE, ultraconserved elements, probe design, phylogenomics, museomics Figures & Tables Figure 1: Step by step comparative flowchart of PHYLUCE and GABBI pipelines. PHYLUCE = Name of Faircloth’s pipeline (2016), GABBI = Genome Alignment Based Bait Inference, developed in this study. a) The first step, where both pipelines differ the most, consists in the identification of targeted loci from a set of chromosome-level genomes to make temporary probes. b) The second step, identical in both pipelines, consists in the validation of the temporary probe set on a set of additional genomes to select the final set of targeted loci. Resulting loci were handed over to myBaits® for actual probe design and filtration. Numerical results and method names written in grey are specific to this study. Figure 2: Proportions of parsimony informative sites between PHYLUCE and GABBI datasets. Markers were splitted into core and flanking regions and analysed accordingly. Core regions annotated as coding were splitted into codon positions 1+2 (pos 1+2) and 3 (pos 3) separately, while non-coding core regions and unannotated flanking regions were left untouched. For each alignment, a proportion of parsimony informative sites (PPIS) is calculated, their distribution in each subcategory is represented with a boxplot and a background violin plot. Distributions of PPIS values of GABBI-derived alignments are shown on the left and PHYLUCE-derived alignments on the right of each pair of boxplots. Asterisks indicate significance levels of pairwise mean comparisons using Wilcoxon-Mann-Whitney tests corrected for multiple comparisons using Benjamini-Hochberg false-rate discovery: *** p < 0.001, **** p < 0.0001. Figure 3: Scaled comparison between Curculionoidea phylogenomic trees obtained with PHYLUCE and GABBI probe sets. Both tree inferences were generated in a ML framework in IQTREE3 using the -m MFP+MERGE option and restricted to markers represented by at least 70% of total taxa. The scale is exactly the same between both trees to highlight differences of tree size. Unsupported nodes are marked with (*) if only UFBS fail or (**) if SH-aLRT and UFBS fall below 80 and 95%, respectively. Main weevil groups are shaded with different colors to better compare topologies and tree lengths. a) Species tree of Curculionoidea obtained following in silico target capture of PHYLUCE-derived markers. UCE = ultra-conserved elements, nt = nucleotides in the supermatrix, MMS = mean marker size (combining core and flanking regions), PVS = proportion of variable sites, PPIS = proportion of parsimony informative sites, TTL = total tree length, IBL = internal branch length and proportion to TTL. b) Species tree of Curculionoidea obtained following in silico target capture of GABBI-derived markers. SHR = shared homologous regions. c) Comparative statistics between PHYLUCE and GABBI-derived trees depending on supermatrix completeness (i.e. the minimum percentage of taxa required to keep a locus). The red rectangle indicates the supermatrix completeness represented. Figure 4: Phylogenomic inference using GABBI weevil probe set combining the Curculionoidea dataset with in vitro target capture specimens. The inferred tree results from a ML framework using IQTREE3 on markers represented by at least 70% of total taxa. Evolution models were computed with the -m MFP+MERGE option. Unsupported nodes are marked with (*) if only UFBS fail or (**) if SH-aLRT and UFBS fall below 80 and 95%, respectively. Main weevil groups are shaded as in figure 3. Specimens obtained from in silico target capture of the Curculionoidea dataset are shown on dark branches whereas specimens obtained from in vitro target capture are shown on red branches. SHR = shared homologous regions, nt = nucleotides in the supermatrix, MMS = mean marker size (combining core and flanking regions), PVS = proportion of variable sites, PPIS = proportion of parsimony informative sites, TTL = total tree length, IBL = internal branch length and proportion to TTL. Credit for pictures: top: Nanodiscus transversus, Martin Galli; bottom: Graptus nictitans, Benjamin Zelvelder. Table 1: Phylogenomic statistics following in silico target capture of currently available Coleoptera probe sets and the new GABBI weevil probe set. The highest values of a given taxonomic dataset are represented in bold. Table 2: DNA extraction and phylogenomic statistics of samples processed with the GABBI-derived weevil probe set. Supplementary Material Figure S1: Scaled comparison between Curculionoidea phylogenomic trees obtained with PHYLUCE and GABBI probe sets. (50% Supermatrix completeness) Figure S2: Scaled comparison between Curculionoidea phylogenomic trees obtained with PHYLUCE and GABBI probe sets. (90% Supermatrix completeness) Figure S3: Scaled comparison between Pachyrhynchini phylogenomic trees obtained with PHYLUCE and GABBI probe sets. (50% Supermatrix completeness) Figure S4: Scaled comparison between Pachyrhynchini phylogenomic trees obtained with PHYLUCE and GABBI probe sets. (70% Supermatrix completeness) Figure S5: Phylogenomic tree obtained with Van Dam et al. ’s (2023) probe set on the Pachyrhynchini dataset. Figure S6: Phylogenomic tree obtained with the GABBI probe set on the Curculionoidea dataset for comparison with existing Coleoptera target capture sets. Figure S7: Phylogenomic tree obtained with the Coleoptera UCE probe set on the Curculionoidea dataset. Figure S8: Phylogenomic tree obtained with the Coleoptera AHE probe set on the Curculionoidea dataset. Figure S9: Phylogenomic tree obtained with the GABBI probe set on the Curculionoidea dataset and Haran et al. ‘s (2023) dataset obtained from Coleoptera AHE probe set. Zenodo supplementary files SVG Figures: - All main text figures are provided as SVG formats - GABBI_pipeline.svg = Figure representing the GABBI pipeline, as uploaded on github Curculionoidea_Dataset_Results: Final results obtained with the Curculionoidea dataset. - prefix Curculionoidea. = 60 genomes from the Curculionoidea dataset - prefix Curculionoidea+12 = 60 genomes from the Curculionoidea dataset + 12 target capture data on in vitro samples - prefix Curculionoidea+AHE = 60 genomes from the Curculionoidea dataset + 30 target capture data conducted with the Coleoptera AHE probe set in Haran et al. (2023b) - GABBI = new GABBI weevil probe set - PHYLUCE = PHYLUCE probe set from this study - UCE = Coleoptera UCE probe set from Faircloth 2017 - AHE = Coleoptera AHE probe set from Haddad et al. 2018 - cds_probes = targeted loci are annotated as coding or non-coding sequences and taken into account for alignment cleaning and supermatrix partitioning using the OMM MACSE pipeline - no_cds = targeted loci are not annotated, no sequences are cleaned with the OMM MACSE pipeline - 50/70/90 = supermatrix completeness, i.e. minimum percentage of taxa required to keep a locus in the supermatrix - suffix summary.txt = AMAS summary at each supermatrix completeness value - suffix MFPMERGE.iqtree = tree statistics given by IQTREE3 - suffix MFPMERGE.treefile = tree file resulting from IQTREE analysis - suffix MFPMERGE.treefile.pdf = graphical representation of the tree file using FigTree, rooting the tree on Anthribidae clade and showing SH-aLRT and UFBS bootstrap values at each node (respectively) - suffix .nex = NEXUS partition file - suffix .phylip = PHYLIP formatted supermatrix file - suffix genetrees_alpha.txt = Alpha parameter value computed for each gene tree with GTR+G8 model - suffix genetrees_summary.txt = AMAS summary of all gene trees - suffix raw_summary.txt = AMAS summary obtained on all alignments before computing gene trees (some alignments present in this file might be missing from gene trees if they have less than 4 taxa or are shorter than 20 bp) Final_Stats_and_Tables: - bbmap_reads.GABBI.tsv = results from "seqkit stats -a" command on bbmap-filtered GABBI reads (named cactus) - bbmap_reads.PHYLUCE.tsv = results from "seqkit stats -a" command on bbmap-filtered PHYLUCE reads (named cactus) - bbmap_reads.tsv = raw results from "seqkit stats -a" command on bbmap-filtered data - perc_stats.Curculionoidea.cds_probe.table = subtable obtained from Table S3 for figure 3 and S1-2 small plots - perc_stats.Pachyrhynchini.cds_probe.table = subtable obtained from Table S3 for figure S3-4 small plots - probe12_multiqc_general_stats.txt = results from multiqc command on in vitro target capture trimmed reads - scaled_trees.R = R script used to generate scaled trees for figures 3 and S1-4 - stats_codon.R = R script used to generate figure 2 and associated statistics - stats_final.sh = Bash script used to generate Table S3 from Curculionoidea_Dataset_Results and Pachyrhynchini_Dataset_Results - stats_final.txt = Unformatted Table S3 obtained from stats_final.sh - summaries.c123.txt = AMAS summaries of codon delimited alignments for figure 2 plot GABBI_Weevil_Probe_Set: - 41g_cactus_v1b.90.anc.loci.final.sorted.fasta = Sequences of all targeted loci and ancestral sequences as given to myBaits - prefix baits-Moderate-RM25pc-0MT-720061count.fas.clust-75-80 = name of the GABBI weevil final probe file generated by myBaits - suffix .loci.fasta = Final set of loci targeted by the GABBI weevil probe set - suffix .soft_cds.ann = Annotation file of the GABBI weevil probe set - GABBI_clust-70-80_mybaits_probes.loci.cons.fasta = Consensus sequences of each targeted locus used to delimit core and flanking regions Genome_Data: Genome assemblies and raw data newly generated in this study and used in the Curculionoidea dataset (see Table S1). Pachyrhynchini_Dataset_Results: Final results obtained with the Pachyrhynchini dataset following the same pattern as the Curculionoidea_Dataset_Results. PHYLUCE_Probe_Set: - PHYLUCE_clust-70-80_mybaits_probes = name of the PHYLUCE final probe file generated by myBaits - PHYLUCE_clust-70-80_mybaits_probes.loci.cons.fasta = Consensus sequences of each targeted locus used to delimit core and flanking regions - PHYLUCE_clust-70-80_mybaits_probes.loci.fasta = Final set of loci targeted by the PHYLUCE weevil probe set - PHYLUCE_clust-70-80_mybaits_probes.soft_cds.ann = Annotation file of the PHYLUCE weevil probe set - PHYLUCE.final.32.anc.loci.sorted.fasta = Sequences of all targeted loci and ancestral sequences as given to myBaits Pipeline_codes: Pipeline of commands used to conduct the analyses - pipeline_annotation.sh = Annotation of targeted loci as CDS or unknown states - pipeline_cactus.sh = GABBI weevil probe set design from genome alignments to in silico target capture of temporary probes. An updated, automated and improved version of this pipeline is available on github: https://github.com/bjzelvelder/GABBI - pipeline_in_silico_test.sh = In silico target capture and phylogenomic analyses of the Curculionoidea datasets - pipeline_Pachyrhynchini.sh = In silico target capture and phylogenomic analyses of the Pachyrhynchini datasets - pipeline_phyluce.sh = PHYLUCE weevil probe set design following the PHYLUCE pipeline Raw_Capture_Data: Raw reads from the targeted enrichment of 12 in vitro samples of weevils using the synthesized GABBI probe set Scripts: Custom scripts used throughout the analyses (called by pipelines in Pipeline_codes). - BBMap_1run.sh = Early version of run_BBMap.sh - clean_genome_headers.py = Remove extra stuff in genome headers - download_genomes_from_ncbi.sh = Download genomes listed in a NCBI table. A new version of this script is available on GABBI's github - get_tiled_probes.py = Manually generate probes using a sliding window - get_tiled_probes.slice.py = Same as get_tiled_probes.py, adapted to different sequence header - mafft_parallel_linsi.sh = Launch multiple MAFFT alignments in parallel on a Slurm server - make_consensus_from_mafft.sh = Launch make_consensus_from_mafft_v2.R on a Slurm server - make_consensus_from_mafft_v2.R = Make different consensus methods based on pairwise distances of a multiple alignment file. A new version of this script is available on GABBI's github - make_consensus.sh = Same as make_consensus_from_mafft.sh, but works in parallel with fixed options - make_partition_from_probe2.py = Generate a partition file based on the coordinates of a reference sequence - phylomera_v0.2.0.sh = Early version of Phylomera (see main text for a description) - phylomera_v0.6.1.sh = Early version of Phylomera (see main text for a description) - phylomera_v0.8.3.sh = Last version of Phylomera (see main text for a description) - postprocess_v2.sh = Script that centralizes a few functions to postprocess IBA results - regroup_matches_from_blastn.py = Regroup multiple pairwise cross-blastn outputs into a single table - regroup_uce.sh = Regroup all assemblies obtained with aTRAM from a given marker into a single fasta file - remove_flank.py = Remove flanking data from a multiple sequence alignment based on a given reference within the MSA - run_aTRAM_v1.sh = Early version of run_aTRAM_v2.sh - run_aTRAM_v2.sh = Run multiple instances of aTRAM in parallel on a Slurm server - run_BBMap.sh = Script to filter raw reads with BBMap - run_IBA_v2.sh = Early version of run_IBA_v3.sh - run_IBA_v3.sh = Run a custom version of IBA on a SLurm server - run_maffilter_phastcons.sh = Small pipeline to run maffilter and phastcons - split_fasta.py = Split a fasta file into X fasta files based on the number of sequences - split_fastq.sh = Split an interleaved paired-end fastq file into two seperate files - transfer_cons_files.sh = Early version of transfer_files.sh (in a loop) - transfer_files.sh = Transfer consensus files obtained after aTRAM into one file per marker, merging all individuals (in parallel)