Range-wide phylogenomics of the Great Horned Owl (Bubo virginianus) reveals deep north-south divergence in northern Peru

Taxon sampling

Twenty-seven samples of Bubo virginianus were assembled from across the range of the species. Specimen-vouchered tissue and toepad samples were loaned from sixteen museum collections (Table 1). This sampling includes thirteen of the sixteen recognized subspecies (Clements et al., 2022); subspecies assignment was based on collection locality and described distributions of subspecies in the literature (Fig. 1). When possible, sampling prioritized specimens taken during the regional nesting season (Table 1). Sequence data for three outgroup species (B. cinerascens, B. nipalensis, and B. scandiacus) were downloaded from the Sequence Read Archive (Salter et al., 2020; Table 1).

DNA extraction and sequencing

Two DNA extraction methods were used, depending on tissue type. DNA was extracted from toepad samples using a Promega Maxwell RSC automated instrument and the Maxwell RSC Tissue DNA Kit. This instrument minimizes potential for contamination, which is important for degraded DNA samples. Fresh tissue samples were extracted using a manual bead-extraction protocol (https://github.com/phyletica/lab-protocols/blob/master/extraction-spri.md), based on Rohland & Reich (2012), with elution of DNA in 1 X TE buffer.

After extraction, we sonicated all fresh tissue samples using a Covaris M220 sonicator (50 W peak incident power, 20% duty factor, and 200 cycles per burst for 65 s). DNA from toepad samples was already relatively degraded, and therefore did not need to be sonicated (McCormack, Tsai & Faircloth, 2016). We then enriched for ultraconserved elements (UCEs) following a standard protocol using the Mycroarray MYbaits kit for Tetrapods UCE 5K v1 (Faircloth et al., 2012). During pooling, on average, each fresh tissue sample received 2.2 million reads per sample and each toepad received 5 million reads per sample. All UCE libraries were sequenced using Illumina paired end 150 bp sequencing on a high output run of a NextSeq550 machine at the KU Genome Sequencing Core. Raw sequence files were deposited in the Sequence Read Archive (Table 1).

UCE sequence data assembly

All UCE data were processed using the standard Phyluce pipeline with the mapping and correction workflows (Faircloth et al., 2012). Briefly, we cleaned raw reads using Illumiprocessor (version 2.0.9), a wrapper for Trimmomatic (version 0.39, Bolger, Lohse & Usadel, 2014) and assembled clean reads using the trinity (version 2.8.5, Grabherr et al., 2011) assembler in Phyluce (version 1.6.8, Faircloth et al., 2012). We then mapped reads onto the resulting contigs using the Phyluce mapping workflow. Bases were then removed if they had a Phred score below 20, a depth below five, or greater than two alleles in called genotypes. These new corrected contigs were then matched with the tetrapods 5K probeset available at ultraconserved.org. These contigs were aligned using MAFFT (version 7.455, Katoh & Standley, 2013), and internally trimmed using gblocks (version 0.91B, Castresana, 2000). We then removed loci that were present in <75% of individuals. These data were split into two datasets, one from fresh tissues only and one with all samples. For the data that included all samples, we wrote a custom script that replaced phylogenetically uninformative sites in the alignment (https://github.com/emilyostrow/ReplaceUninformativeSites) due to potential data degradation in toepad samples. In addition to missing data, toepad samples tended to have erroneous bases sequenced near cut sites. This script replaced any base that is represented only once at a particular site with an N to address branch length concerns caused by degraded DNA.

After generating UCEs, we called single nucleotide polymorphisms (SNPs) using the consensus UCE sequences as a reference with the GATK pipeline. We used Geneious Prime (version 2023.0.4, https://www.geneious.com) to generate consensus sequences on all UCEs found in at least 75% of the individuals using a 50% strict cutoff for each site. To call SNPs, we first trimmed raw read files using AdapterRemoval (version 2.3.2, Schubert, Lindgreen & Orlando, 2016). We then built an index of the UCE consensus sequences in Bowtie2 (version 2.3.5.1, Langmead & Salzberg, 2012) and aligned reads to the reference using the very-sensitive-local algorithm. After aligning the reads, added read group information and marked duplicates in picard (version 2.20.3, http://broadinstitute.github.io/picard) and indexed the files using SAMtools (version 1.9 Li et al., 2009) We then used GATK (version 4.2.6.1, McKenna et al., 2010) to call haplotypes, genotype sequences, select SNPs, and filter SNPs. We hard filtered SNPs for quality by depth (QD) < 2.0, strand odds ratio (SOR) > 4.0, Fisher strand (FS) > 60.0, and RMS mapping quality (MQ) < 40.0. We then filtered the SNPs again using SNPfiltR (version 1.0.0, DeRaad, 2022) in R (version 4.1.1, R Core Team, 2021) for 90% or greater data completeness by SNP, phred score of 30 or greater and only included biallelic SNPs.

Mitochondrial genome assembly

As a byproduct of UCE sequencing, the mitochondrial genome is often also sequenced (Do Amaral et al., 2015). To obtain mitochondrial genomes, we used the program MITObim following the two-step procedure outlined in the manual (version 1.9.1, Hahn, Bachmann & Chevreux, 2013). We used cleaned reads from the Illumiprocessor step in Phyluce as input data and a B. scandiacus sample (NC_038220.1) from NCBI as a reference mitochondrial genome for alignment. Briefly, we interleaved reads one and two from the cleaned Illumiprocessor data. We created an initial read map using MIRA (version 4.0.2, Chevreux, Wetter & Suhai, 1999), then completed iterative mapping using MITObim. MITObim ran up to ten iterations or until it reached a stationary read number.

All individual mitochondrial genome alignments were then aligned to the reference B. scandiacus sample using MAFFT in Geneious (version 8.1.9, Kearse et al., 2012) to create multi-sample alignments. Gaps in alignment created by a single individual were removed by hand. Non-genic regions were removed from the multi-sample alignment because more variable areas such as the control region did not align well and may have led to inaccurate results.

Phylogenetics

Phylogenetic analyses were completed with the UCE data using maximum-likelihood, quartet-based, and neighbor-joining approaches. For the maximum likelihood method, we used the program SWSC-EN (Tagliacollo & Lanfear, 2018) to split UCE data into conserved core and more variable flanking regions for each UCE. These regions were then tested as potential partitions in PartitionFinder2 using a rclusterf search scheme (version 2.1.1, Lanfear et al., 2017). We used the best partitioning scheme according to the AICc score with a GTR+G model in RAxML-NG (version 1.1.0, Kozlov et al., 2019) using both the only-tissues dataset and the full sampling with no uninformative sites. We completed only 200 bootstrap replicates of the tissues-only matrix due to bootstrap convergence and 1,000 bootstrap replicates of the full sampling matrix. Our quartet-based approach used SVDquartets in PAUP* (version 4.0a, Chifman & Kubatko, 2015) on default settings with 1,000 bootstrap replicates. We generated a neighbor-joining tree using the dataset with all individuals and no uninformative sites in Geneious Prime. We completed 1,000 bootstraps and collapsed nodes with a support threshold below 50%.

Our mitochondrial phylogenetic analyses included both maximum likelihood and Bayesian approaches. We used PartitionFinder2 to test models for both tree-building approaches. We split all genes by codon and tested all models of evolution using a greedy search scheme. Models and partitions were evaluated using AICc scores. Our maximum likelihood analysis was completed with RAxML-NG, using the partitions from PartitionFinder2 and a GTR+G model with 200 bootstrap replicates, after which the tree converged. Our Bayesian analysis was completed using the program MrBayes (version 3.2.7a, Ronquist & Huelsenbeck, 2003). We used the optimal partitions and models from PartitionFinder2. Our analysis used four chains with 10 million generations, using a burn-in fraction of 25%, sampling every 1,000 generations. We checked to ensure convergence for the Bayesian tree using potential scale reduction factor and effective sample size estimates.

Gene flow

We tested for gene flow between populations using Dsuite (version 0.4 r42, Malinsky, Matschiner & Svardal, 2021). Specifically, we calculated the D-statistics and f4-ratios using the Dtrios algorithm using the vcf file and full sampling UCE tree. We then used the Fbranch algorithm to calculate the F-branch statistic for each individual and ancestor comparison. We also calculated p-values associated with the F-branch statistics and corrected for multiple testing using the False Discovery Rate Benjamini-Hochberg procedure (Benjamini & Hochberg, 1995). We visualized these results and support values using dtools.

Audio recording sonograms

We downloaded six audio recordings of B. virginianus songs with minimal background noise from Xeno-Canto (https://xeno-canto.org): XC428421 from Minnesota, USA, XC511167 from New York, USA, XC548133 from California, USA, XC76397 from Napo, Ecuador, XC212734 from Salta, Argentina, and XC494437 from Santiago, Chile. Audio recordings were chosen based on audio quality and geographic breadth. We imported these songs into Adobe Audition (Build 14.1.0.43) and iteratively used the denoise tool to reduce background noise in the sonograms until the song was clearly differentiated from the background noise. We then exported audio files from Adobe Audition and imported them into Raven Lite (version 2.0.1, K. Lisa Yang Center, 2023) where we visualized each song as a sonogram, adjusted contrast as needed for clarity, and saved each song as an image file.

UCE and mitochondrial data

We generated a mean of 4248 UCE loci per individual. Fresh tissue samples yielded a mean of 4318 UCE loci per individual, whereas toepad samples yielded a mean of 4127 loci per individual. The 75% matrix (i.e., loci present in ≥22 of the 30 individuals) included 4299 UCE loci, with a 1.26 Mbp overall alignment length and an average UCE length of 796 bp for fresh tissues and 299 bp for the toepad samples. Our unfiltered SNP dataset included 7584 SNPs. After filtering for SNP completeness, phred score, and biallelic sites, our final dataset included 6556 SNPs.

The mitochondrial baiting method was successful for all individuals in this study. MITObim contigs had a mean read depth of 101X and a mean contig length of 18,967 bp. Fresh tissue samples had a mean read depth of 95X, whereas toepad samples had a mean read depth of 111X, presumably due to increased sequencing effort dedicated to toepad samples. After trimming to genic regions only, we generated an 11,347 bp alignment, with 1.09% missing data. One individual in particular, B. v. virginianus (FMNH 396867), had 23.05% missing data and disproportionately contributed to the overall missing data in the mitogenome matrix. The effects of these missing data and degradation, contamination, or sequencing error at the edges of gaps are apparent in the longer branch in the mitochondrial tree.

Phylogenetics

Maximum-likelihood, quartet-based, and neighbor-joining methods using the nuclear data produced trees that were congruent at most moderate to highly supported nodes (bootstrap support above 90%); the trees mainly differed in topology at nodes with low support. The one topological difference with bootstrap support in the quartet-based tree (Fig. S1) was the placement of B. v. nigrescens (FMNH 102981), which was sister to B. v. nacurutu (KU 90879) + B. v. nigrescens (USNM 368745) in the SVDquartets tree whereas in all other full sampling trees, B. v. nigrescens (FMNH 102981) was sister to B. v. nigrescens (LSUMZ 30041). The one topological change between the maximum-likelihood tree and neighbor-joining tree was the placement of B. v. nacurutu (KU 90879) + B. v. nigrescens (USNM 368745), which was sister to B. v. nigrescens (FMNH 102981) + B. v. nigrescens (LSUMZ 30041) in the neighbor-joining tree (Fig. S2) but is within the mostly Central American clade in the maximum likelihood tree. Here, we present the maximum-likelihood tree (Fig. 2). The quartet-based SVDquartets tree had few nodes with high support (Fig. S1). However, both RAxML and SVDquartets trees showed strong support for a B. v. magellanicus clade sister to all other samples (hereafter referred to as the northern clade). Within the northern clade, lineages branching from earlier nodes belonged to more southern taxa (B. v. nacurutu and B. v. nigrescens) with less geographically concordant branching order in other northern subspecies.

Maximum-likelihood and Bayesian analyses using the mitochondrial dataset generally supported the same topology (Fig. 3, Fig. S3). Two discrepancies separated the trees, both at extremely shallow nodes: B. v. mayensis (FMNH 187122) and B. v. pinorum (UWBM 74120) switched positions, and B. v. elachistus (DMNS 21432) and B. v. pacificus (LACM 115875) switched positions. Support values from both analyses are mapped onto the RAxML tree at concordant nodes (Fig. 3).

The trees built using UCE and mitochondrial data were not entirely concordant. These inconsistencies might be due to the different inheritance of nuclear and mitochondrial DNA, gene flow, taxon sampling differences in the UCE trees, or potentially greater errors in toepad sequencing due to their degraded DNA. We attempted to address potential data errors by using the Phyluce correction pipeline and only using Phylogenetically informative sites to reduce spurious branch lengths. All trees shared a general pattern of southern lineages branching from earlier nodes, but the topologies contained strongly supported differences. The relationships of many of the taxa in southern Mexico, Central America, and northern South America are different in all reported trees. One major difference in the topologies is the relative placement of B. v. nacurutu (KU 90879) and B. v. nigrescens (LSUMZ 30041). In the fresh tissues only UCE tree, they are sequential sister taxa to the remainder of the northern clade, whereas in the all samples UCE tree, B. v. nacurutu (KU 90879) is found within the larger North American clade (Fig. 2). In the mitochondrial tree, B. v. nigrescens (LSUMZ 30041) and B. v. nacurutu (KU 90879) switch branching order as compared to the UCE trees (Fig. 3). Similar to the UCE data, most B. v. magellanicus in the mitochondrial trees form their own clade divergent from all other subspecies. This clade is 5.95% divergent from the northern clade using the mitochondrial alignment that included all genes. However, B. v. magellanicus MSB 35949 falls within the northern clade and is one of three northernmost B. v. magellanicus samples.

Gene flow

Significant F-branch statistics between B. v. nigrescens (LSUMZ 30041) from Ecuador and Peruvian B. v. magellanicus support Nuclear gene flow between B. v. magellanicus and non-magellanicus B. virginianus populations in the Andes (Fig. S4). We do not have the sampling to assess gene flow between B. v. magellanicus and other B. virginianus subspecies at the eastern limit of the B. v. magellanicus range. We also identified consistent gene flow between B. v. nigrescens (FMNH 102981) and much of the northern clade (Fig. S4). Other areas of the tree also showed evidence of gene flow, but gene flow is expected at the subspecific level (Price, 2008).

Audio recording sonograms

Sonograms of the six audio recordings from Xeno-Canto illustrate the vocal differences between populations of B. v. magellanicus and northern B. virginianus described by López-Lanús (2015) (Fig. 4). The two recordings from B. v. magellanicus had two notes followed by approximately a two second trill. The song from the four northern B. virginianus recordings had three to five notes. These songs are clearly distinguishable without formal vocal analyses due to the strong differentiation in amount and length of notes.