Whole-Genome Analysis of Introgression Between the Spotted Owl and Barred Owl (Strix occidentalis and Strix varia, Respectively; Aves: Strigidae) in Western North America

As the barred owl (Strix varia; Aves: Strigiformes: Strigidae) expands throughout western North America, hybridization between barred and spotted owls (Strix varia and S. occidentalis, respectively), if abundant, may lead to genetic swamping of the endangered spotted owl. We analyzed low-coverage, whole-genome sequence data from fifty-one barred and spotted owls to investigate recent introgression between these two species. Although we obtained genomic confirmation that these species can and do hybridize and backcross, we found no evidence of widespread introgression. Plumage characteristics of western S. varia that suggested admixture with S. occidentalis appear unrelated to S. occidentalis ancestry and may instead reflect local selection.

KEYWORDS bird avian hybridization raptor population genomics admixture Over the past century, humans have introduced several non-native vertebrate species in western North America into the native range of closely related species and generated moving hybrid swarms. For example, in California, genes of the non-native barred tiger salamander (Ambystoma tigrinum) are spreading into the range of the California tiger salamander (A. californiense) (Fitzpatrick et al. 2009(Fitzpatrick et al. , 2010. In the Flathead River system of Montana and British Columbia, the non-native rainbow trout (Oncorhynchus mykiss) is rapidly hybridizing with the native westslope cutthroat trout (O. clarkii lewisi) (Muhlfeld et al. 2014).
In addition to hybridization resulting from intentional introductions of non-native species, changing global climatic conditions and the documented movement of species ranges have led many species to invade novel geographic regions (Parmesan et al. 1999;Parmesan 2006) and establish broad contact with related taxa (Rieseberg et al. 2007).
The spotted owl (Strix occidentalis) is a large wood owl inhabitant of western North American forests. The U.S. Fish and Wildlife Service listed the northern spotted owl (S. o. caurina) as "threatened" under the Endangered Species Act (ESA) in 1990 (Thomas et al. 1990) and the species remains protected due to continuing population declines (Dugger et al. 2015;Davis et al. 2016). While researchers considered habitat loss the primary threat to the northern spotted owl in 1990 (Forsman et al. 1984;Anderson and Burnham 1992), recent research has confirmed a second major threat to its persistence: the invasion of the congeneric barred owl 50-100 years (Dark et al. 1998;Livezey 2009aLivezey , 2009b. At present, sympatric populations of spotted and barred owls exist from British Columbia to southern California (Taylor and Forsman 1976;Haig et al. 2004;Livezey 2009a).
Strix occidentalis and S. varia are approximately 13.9% divergent in the mitochondrial control region (Haig et al. 2004), 10.74% divergent in non-tRNA mitochondrial genes (Hanna et al. 2017a), and 0.7% divergent across the nuclear genome (Hanna et al. 2017c). Strix occidentalis and S. varia hybridize and backcross (Haig et al. 2004;Kelly and Forsman 2004;Funk et al. 2007), with heterospecific matings and F 1 hybrids commonly reported in areas where S. varia is rare and S. occidentalis is common (Kelly and Forsman 2004). Recently collected adult western S. varia specimens in museum collections display striking morphological variation. Strix varia from the eastern Klamath Mountains in Siskiyou County, California, have darker plumage overall, more spotting on the belly, and are smaller than barred owls from the Coast Range ( Figure 1 and Figure S1). These differences suggest either local selection for this phenotype or possible introgression of spotted owl genes. Hybridization of these species creates a potential for a loss of biodiversity in western North America due either to replacement of the spotted owl by the barred owl or to collapse of the boundaries of the two species (Huxel 1999).
For this study, we obtained fifty-one low-coverage whole-genome sequences (median 0.723X coverage) from Strix varia and S. occidentalis sampled outside and across their contact zone in western North America ( Figure 2). Existing methods available for detecting and quantifying admixture in genomic data require either specific genotype calls or accurate genotype likelihood estimates (e.g., Skotte et al. 2013;Wall et al. 2016;Corbett-Detig and Nielsen 2017;Meisner and Albrechtsen 2018). For most individuals, the vast majority of variants had a depth of coverage of 0 or 1. In this context, genotype likelihoods or called genotypes are not meaningful and we were unable to utilize existing methods that rely on them. This necessitated that we develop a new method for quantifying admixture that only requires single reads at various sites. We utilized available medium and high-coverage wholegenome sequences from an eastern S. varia (15.549X coverage) and a pre-contact S. occidentalis (60.815X coverage) to identify variant sites potentially fixed between S. varia and S. occidentalis. For each lowcoverage individual, we determined the genome-wide average ancestry and searched for windows of ancestry that were outliers from the average to detect rare, introgressed regions. We used these data to identify the extent of introgression between S. varia and S. occidentalis in western North America.

Samples
We obtained fifty-one samples from museum collections that included eleven Strix occidentalis samples (two samples predated contact with S. varia), thirty-eight S. varia samples (including five from eastern North America), and two samples identified by other researchers as probable hybrid S. varia x occidentalis individuals (Tables S1-S2). We mapped the samples using QGIS version 2.18.2 (Quantum GIS Development Team 2017) with raster and vector files from Natural Earth (http://www.naturalearthdata.com) (File S1 section 1.1).

Sequence data
Hybridization of Strix varia and S. occidentalis has previously been investigated using a set of four microsatellite (Funk et al. 2007) and fourteen amplified fragment-length polymorphism (Haig et al. 2004) markers, which the authors found useful for diagnosing F 1 and F 2 hybrids (Haig et al. 2004;Funk et al. 2007). We utilized a whole-genome sequencing approach to enable us to detect any introgression that has taken place over the last 50-70 years that S. varia and S. occidentalis have been in contact in western North America (Taylor and Forsman 1976;Livezey 2009a).

Alignment and filtering
For the sequence data of the reference samples Strix occidentalis CAS: ORN:98821 and S. varia CNHMB41533, which Hanna et al. (2017c) generated for their study, we followed the sequence processing methods described in Hanna et al. (2017c) to trim and filter the genomic sequence data. For all other samples we used Trimmomatic version 0.32 (Bolger et al. 2014) to remove adapter sequences and perform quality trimming of all of the low-coverage, short-read data (File S1 section 1.3). We used BWA-MEM version 0.7.12-r1044 (Li 2013) to align the processed reference and low-coverage sequences to the repeat-masked S. o. caurina genome "StrOccCau_1.0_nuc_masked" (Hanna et al. 2017d(Hanna et al. , 2017c. We merged the alignments, sorted the alignments, and marked duplicate sequences using Picard version 1.104 (http:// broadinstitute.github.io/picard; accessed 2017 Oct 1) (File S1 section 1.4.1-1.4.2). We filtered the alignment files to only retain alignments of high quality using the Genome Analysis Toolkit ( with the alignment files for all samples included as simultaneous inputs (File S1 section 1.5.1). We used the vcf_qual_filter.sh script from SPOW-BDOW-introgression-scripts version 1.1.1 (Hanna et al. 2017b) to exclude indels and low genotyping quality sites while retaining only biallelic sites where CAS:ORN:98821 (the source of the StrOcc-Cau_1.0_nuc_masked reference genome) was homozygous for the reference allele and CNHMB41533, the Strix varia reference sample, was homozygous for the alternative allele (File S1 section 1.6.1). Of the remaining variable sites, we excluded those with excessively high coverage [greater than the mean plus five times the standard deviation (s), as recommended by the GATK documentation (https://software. broadinstitute.org/gatk/guide/article?id=3225; accessed 2017 Oct 1)] (File S1 section 1.6.2). We used DP_means_std_dev.sh from SPOW-BDOW-introgression-scripts version 1.1.1 to calculate the mean and standard deviation (s) of the depth of coverage for each sample across the final set of variant sites.

Ancestry and diversity analyses
For each sample at each of the final variant sites, we calculated a percentage spotted owl ancestry, which was the percentage of the coverage that supported the CAS:ORN:98821 (the Strix occidentalis reference sequence) allele. We calculated the mean and standard deviation of the spotted owl ancestry of each sample across all variant sites (File S1 section 1.6.3). We tested for significant differences between the mean spotted owl ancestries in populations using Welch's t-test (Welch 1947) as the populations had unequal numbers of samples and then applied a Bonferroni adjustment (Dunn 1961) when we evaluated significance (File S1 section 1.6.4).
We estimated the probabilities of observing an introgressed region greater than 50,000 nt, 100,000 nt, or 150,000 nt in length if Strix varia and S. occidentalis hybridized in 1945, approximately the earliest date of their potential contact (Livezey 2009a), using the formula from Racimo et al. (2015). For the recombination rate, we used 1.5 centimorgans/ million nucleotides (cM/Mnt), which Backström et al. (2010) estimated for the zebra finch (Taeniopygia guttata). For the number of generations since the earliest potential date of hybridization, we assumed a generation time of two years (Gutiérrez et al. 1995;Mazur and James 2000) even though S. o. caurina is able to breed in its first year and others have used ten years as the generation time for S. o. caurina (Noon and Biles 1990;USDA Forest Service 1992). With that generation time, approximately thirty-five generations have potentially elapsed since the two species first contacted in 1945 and 2014, the date of our most recent sample.
In order to probe further for evidence of introgression in the samples that did not appear as hybrids from their genome-wide average spotted owl ancestry, we attempted to identify regions that were outliers from the genome-wide ancestry average by conducting a sliding window analysis. We examined adjacent windows of 50,000 nucleotides (nt) where a sample had data for at least ten variant sites within that window and calculated the average spotted owl ancestry for the window. We assumed that, if a region was introgressed from the other species, the average should be close to 0.5. Thus, in samples with an average genome-wide ancestry close to 0, we called a window an outlier if the average spotted owl ancestry was .= 0.4. Inversely, in samples with an average genomewide ancestry close to 1, we called a window an outlier if the average spotted owl ancestry was ,= 0.6 (File S1 sections 1.6.5-1.6.6).
In order to estimate the genome-wide diversity harbored by Strix varia and S. occidentalis populations, we considered all biallelic variant sites (not just those fixed between our S. varia and S. occidentalis references) and calculated p Within , the number of nucleotide differences within populations, and p Between , the number of nucleotide differences between populations using the countFstPi script from SPOW-BADOintrogression-scripts (Hanna et al. 2017b). We also used countFstPi to calculate the fixation index (F ST ) (Hudson et al. 1992) in order to estimate the differentiation of S. varia and S. occidentalis populations (File S1 section 1.6.7).

RESULTS
After filtering, the final set of variable sites fixed between the Strix varia and S. occidentalis reference individuals included 5,816,692 sites. The median genome coverage per individual was 0.723X (Table S3). Except for the two putative hybrid samples that we included as a test of our methodology, the genome-wide average spotted owl ancestry for all samples was close to either 0 or 1, indicating that they were either pure S. varia or S. occidentalis, respectively ( Figure 3 and Table S3). A genome-wide average spotted owl ancestry of 0.538 confirmed the F 1 hybrid (S. varia x occidentalis) identity of a sample from Humboldt County, California. We calculated a spotted owl ancestry of 0.359 for the second hybrid sample from Benton County, Oregon, which suggested that this individual was likely a F 2 hybrid (F 1 x S. varia backcross). The "Specimen Identifier" column provides the voucher specimen codes. The "Other Sample Identifier" column provides an abbreviated sample code. Column "Sample Set" refers to the round of sequencing that produced the sequence data for a given sample. The main and supplemental methodology sections provide details of the production of these two sets of sequence data. Column "SRA ACCN" provides NCBI Sequence Read Archive (SRA) run accessions in which the raw sequences for each sample are archived.
The mean genome-wide spotted owl ancestry of the Siskiyou County Strix varia population was 0.0696 whereas the mean was 0.0699 for the rest of the western S. varia (Table S4). There was no significant difference in spotted owl ancestry between these two populations (Table S5). When we combined all S. varia from western North America together (0.0698 mean spotted owl ancestry) and compared their spotted owl ancestry with that of the eastern S. varia (0.0676 mean spotted owl ancestry), we found no significant difference in ancestry between the western and eastern S. varia after applying a Bonferroni adjustment (Tables S4 and S5). There was also no significant difference in spotted owl ancestry between S. occidentalis individuals sampled from populations not in contact with S. varia and those from populations already in contact with S. varia (mean ancestries of 0.9930 and 0.9952, respectively) (Tables S4 and S5).
The average spotted owl ancestry in the Strix varia samples ranged from approximately 6.55-7.28% greater than the 0% value at which our methodology set the reference S. varia (Table S3). The S. occidentalis samples ranged from approximately 0.43-0.94% less than the 100% value for the reference S. occidentalis. The standard deviation in the S. varia samples was consistently more than two times greater than the standard deviation in the S. occidentalis samples. The genome-wide average spotted owl ancestry values for the Strix varia individuals deviated more from those of the reference S. varia than did the S. occidentalis individuals from the S. occidentalis reference due to the greater amount of genetic variation within S. varia (Hanna et al. 2017c). The sites fixed between our reference S. varia and S. occidentalis samples were not fixed across all of the low-coverage S. varia and S. occidentalis individuals. Further high-coverage sequencing of whole-genomes for both species will help to more clearly identify the fixed genetic differences between the two species.
Based upon an estimate of thirty-five generations as the maximum number of generations since contact of Strix varia and S. occidentalis (Gutiérrez et al. 1995;Mazur and James 2000;Livezey 2009a) and the recombination rate of Taeniopygia guttata (Backström et al. 2010), we estimated that the probability of observing a tract . 50,000 nt resulting from hybridization during the initial contact of S. varia and S. occidentalis was 97.41%, the probability of observing a track . 100,000 nt was 94.89%, and the probability of observing a track . 150,000 nt was 92.43%.
Of the forty-nine samples for which we conducted an outlier window analysis, we detected outlier windows in thirty-nine samples (79.6%). Across all samples, we detected 316 outlier windows of length 50,000 nt, forty-one of length 100,000 nt, and only three of length 150,000 nt and none exceeded this length ( Figure S2). In all samples the outlier windows represented , 1.01% of the analyzed windows. For thirty-six of the thirty-nine samples with outliers, the number of outlier windows was , 0.08% of the analyzed windows. There were three samples for which the outlier windows represented between 0.1% and 1.01% of the analyzed windows. However, the increased proportion of outlier windows in these samples appeared to be related to exceptionally low sequence coverage as these three Strix varia samples had the lowest coverage (0.036-0.118X) and, consequently, the fewest number of analyzed windows of any of the samples in which we detected outlier windows ( Figure S3). A S. occidentalis sample with 0.017X coverage was the only sample with lower coverage than those three, but our analyses did not recover any outlier windows for it.
We found little evidence of differentiation between the Siskiyou Strix varia and the other western S. varia, recovering a low F ST (0.008) and very similar levels of nucleotide diversity in the two populations (Table 2). Similar levels of nucleotide diversity also exist in the S. varia populations from western and eastern North America. We additionally estimated a low F ST value (0.051) between western and eastern S. varia, which suggests a low level of differentiation between these populations. Strix occidentalis populations pre and post-contact with S. varia exhibited similar levels of nucleotide diversity and appeared weakly differentiated (F ST = 0.022). We estimated approximately 14X greater nucleotide diversity in S. varia than S. occidentalis and a high level of divergence (F ST = 0.833) between the species.

DISCUSSION
Our genome-wide average spotted owl ancestry analysis confirmed that our two positive control hybrids from Humboldt County, California, and Benton County, Oregon, were an F 1 and F 2 (F 1 x Strix varia) backcross, respectively. Apart from those hybrids, our genome-wide average spotted owl ancestry analysis indicated that all individuals were either pure S. occidentalis or pure S. varia (Figure 3 and Table S3). Our global analysis found no evidence for admixture, but we cannot rule out low (e.g., , 2%) levels of admixture. We also implemented a sliding window approach to determine whether any such regions existed in our data. Scanning for ancestry windows that were outliers from a given individual's genome-wide average ancestry using a sliding window approach corroborated the genome-wide average results and provided no evidence of introgression between Strix varia and S. occidentalis within the past 50-70 years of their contact in western North America (Taylor and Forsman 1976;Livezey 2009a). Although our test found short windows of outlier ancestry, these represented a small proportion of the total windows analyzed for each individual. Thus, we can confidently exclude the possibility of introgression within the past ten generations. Hybridization that has occurred in the last thirty-five generations (assuming a generation time of two years for both species and erring conservatively on the side of overestimating the maximum number of generations of contact) should have yielded much longer outlier blocks than we found. Even with this conservative estimate, there is a . 97% probability of introgressed regions being larger than the 50,000 nt windows that we used to check for potential introgression and a . 92% probability of the introgressed regions being larger than the 150,000 nt length of the longest outlier window that we detected with our sliding window analysis. Figure 3 Plot of coverage vs. genome-wide average spotted owl ancestry. The average spotted owl (Strix occidentalis) ancestry of all of the samples for which we collected low-coverage, whole-genome sequence data. We plotted DNA sequence coverage on the y-axis to display that the average percentage of spotted owl ancestry was independent of the amount of coverage for a given sample.
Since Strix varia's zone of contact with S. occidentalis in western North America began in British Columbia and expanded southward to the southern Sierra Nevada, California (Taylor and Forsman 1976;Haig et al. 2004;Livezey 2009a), we expected S. varia individuals in the southern portion of the zone of sympatry to have the highest chance of being admixed. With this prediction in mind, we focused our sampling on S. varia populations in California ( Figure 1) and targeted our sampling to include the morphologically anomalous western S. varia population in Siskiyou County, California. It is notable that we found no evidence of admixture even though these populations visually appeared intermediate in plumage between S. varia and S. occidentalis. Range expansion simulations suggest that we should predict asymmetric introgression into S. varia even when the hybridization rate is less than 2% (Currat and Excoffier 2011). Coupled with these predictions, our findings suggest that, although hybridization between S. varia and S. occidentalis occurs, it has either been vanishingly rare on the edge of the S. varia expansion wave or other processes, such as selection or migration, are effectively removing introgressed genetic material from S. varia and S. occidentalis populations.
We estimated that Strix varia has more than ten times greater nucleotide diversity than S. occidentalis and we calculated a high F ST between the species (Table 2), closely matching results from highcoverage genomes of the two species (Hanna et al. 2017c). We estimated similar levels of nucleotide diversity in the Siskiyou S. varia population and the population comprised of other western S. varia, which was consistent with our having found no difference in spotted owl ancestry between these populations (Tables S4 and S5). Similarly, S. occidentalis populations pre and post-contact with S. varia exhibited similar levels of nucleotide diversity, appeared weakly differentiated, and did not differ in spotted owl ancestry.
We were surprised to find similar levels of nucleotide diversity in western and eastern North American Strix varia populations. We expected western S. varia populations to harbor lower genetic diversity than the eastern S. varia after having been subjected to successive founder effects and corresponding reductions in nucleotide diversity (Austerlitz et al. 1997). Simulations have suggested that long-distance dispersal by individuals of a species undergoing a range expansion can inhibit the loss of genetic diversity in the newly formed populations on the edge of the range (Ray and Excoffier 2010). Engler et al. (2015) suggested that this explains why some populations retained genetic diversity in an Old World warbler, Hippolais polyglotta, experiencing a range expansion. Recent simulations have also suggested that longdistance dispersal in an invading taxon can counteract introgression of local genetic material into the invader by inhibiting the "surfing" of introgressed genetic regions (Amorim et al. 2017). Livezey (2009b) reported the mean natal dispersal distance of Strix varia as 41.3 km, but mentioned that some individuals have dispersed as far as 488.1 km. Even if long-distance dispersal has only been occurring at low levels during the S. varia range expansion, this could account both the lack of reduction in genetic diversity in western S. varia and for the lack of large-scale introgression of S. occidentalis genetic material into western S. varia populations (Ray and Excoffier 2010;Amorim et al. 2017). Long-distance dispersal would have been especially capable of countering introgression of S. occidentalis material if non-introgressed S. varia were dispersing to the front of the expansion wave (Amorim et al. 2017). Long-distance dispersal may also lead to high rates of intraspecific gene flow in western S. varia, which could both maintain S. varia genetic diversity and counter introgression of S. occidentalis genetic material (Ray et al. 2003;Currat et al. 2008;Petit and Excoffier 2009).
Although our results provide genomic confirmation that hybridization and backcrossing does occur, we found no evidence of widespread admixture between Strix varia and S. occidentalis in western North America. The distinctive plumage of the S. varia individuals collected in Siskiyou County, California, (Figure 1 and Figure S1) does not appear to be a result of hybridization with S. occidentalis. We conclude that some plumage characteristics that appear intermediate between S. varia and S. occidentalis do not in fact indicate hybridization. Previous investigators have issued similar cautionary statements after their genetic studies of hybridization in these taxa (Haig et al. 2004;Funk et al. 2007). The lack of spotted owl ancestry in these oddly plumaged western S. varia suggests that some western S. varia may be undergoing drift or local selection, which has affected plumage and size. Coupled with demographic studies (Kelly et al. 2003;Dugger et al. 2015;Diller et al. 2016), our results indicate that the expansion of S. varia into the range of S. occidentalis in western North America is following a pattern of pure replacement, rather than inducing extinction through hybridization and introgression (Rhymer and Simberloff 1996). It seems unlikely that even introgressed remnants of the S. occidentalis genome will remain in areas in contact with S. varia if S. occidentalis is not able to persist.  The p Within statistic signifies the average number of pairwise differences between two individuals sampled from the same population. The p Between statistic denotes the average number of pairwise differences between two individuals sampled from different populations (Populations 1 and 2). "Pop 1" and "Pop 2" refer to Population 1 or 2 from columns 1 and 2, respectively. The "All Western Barred Owls" population is a combination of the "Western Barred Owls" and "Siskiyou Barred Owls" populations. The "Spotted Owls (pre-contact)" and "Spotted Owls (post)" populations indicate Strix occidentalis from populations not in contact or in contact with S. varia, respectively.

ACKNOWLEDGMENTS
Museum Center; Irby Lovette and the Cornell University Museum of Vertebrates; Maureen Flannery, Laura Wilkinson, and the California Academy of Sciences; and Carla Cicero, Theresa Barclay, Shelby Medina, Elizabeth Wommack, and the Museum of Vertebrate Zoology. We thank Anna Sellas for assistance with laboratory work. We generated genomic libraries at the Center for Comparative Genomics, California Academy of Sciences.