Complete mitochondrial genome sequences of the northern spotted owl (Strix occidentalis caurina) and the barred owl (Strix varia; Aves: Strigiformes: Strigidae) confirm the presence of a duplicated control region

Introduction

The chicken (Gallus gallus) was the first avian species with a complete mitochondrial genome assembly (Desjardins & Morais, 1990). Subsequently, researchers assembled the mitochondrial genomes of members of the Paleognathae (e.g., ostriches, emus, kiwis) and other members of the Galloanserae (ducks, chicken-like birds) and recovered the same gene order found in the mitochondrial genome of the chicken, which led to the conclusion that the mitochondrial genome of the chicken is representative of the ancestral avian gene order (Desjardins & Morais, 1990; Mindell, Sorenson & Dimcheff, 1998a; Haddrath & Baker, 2001; Gibb et al., 2007). Almost a decade after publication of the chicken mitochondrial genome, Mindell, Sorenson & Dimcheff (1998a) described an alternative or, to use their terminology, “novel” avian gene order from that of the chicken, which included a different positioning of tRNA^Pro, ND6, and tRNA^Glu relative to the control region sequence as well as an additional noncoding segment that they hypothesized was a degraded copy of the control region. A few years later, researchers first described the presence of an intact, duplicated control region in the mitochondrial genomes of Amazona parrots (Eberhard, Wright & Bermingham, 2001) and the common buzzard Buteo buteo (Haring et al., 2001).

Mindell, Sorenson & Dimcheff (1998a) detected their novel avian gene order in the mitochondrial genomes of taxa in multiple avian orders that spanned a significant portion of Neoaves, but did not detect it in the single owl species that they studied, Otus asio (Mindell, Sorenson & Dimcheff, 1998a). However, further investigation of owl (Strigiformes) mitochondrial genomes has revealed several surprises.

First, at least three wood owl species (Strix aluco, S. uralensis (Brito, 2005), and S. varia (Barrowclough et al., 2011)) contain the novel mitochondrial gene order of Mindell, Sorenson & Dimcheff (1998a) as well as duplicate control regions. The researchers’ use of a primer in tRNA^Thr to amplify a fragment of the control region suggests that the novel gene order is present in two additional wood owl species, S. occidentalis (Barrowclough, Gutierrez & Groth, 1999) and S. nebulosa (Hull et al., 2010). However, the novel gene order was not reported as present in the mitochondrial genome of S. leptogrammica (Liu, Zhou & Gu, 2014).

Second, some species of eagle-owls (genus Bubo) have a large control region (up to ∼3,800 nucleotides) relative to Strix, their putative sister genus (Fuchs et al., 2008; Wink et al., 2009), largely due to a tandem repeat structure in the 3′ end of the control region (Omote et al., 2013). Such control region tandem repeat blocks appear to be widespread in Strigidae (Xiao et al., 2006; Omote et al., 2013). These results suggest that the structures of owl mitochondrial genomes are surprisingly dynamic and in need of further investigation, particularly for species of conservation concern for which portions of the control region are used in population genetic studies (Barrowclough, Gutierrez & Groth, 1999; Haig et al., 2004; Hull et al., 2010; Hull et al., 2014).

We here provide the complete mitochondrial genome sequence of both a northern spotted owl (Strix occidentalis caurina) and barred owl (S. varia). The spotted owl (S. occidentalis) is a large and charismatic denizen of dense forests whose range includes the Pacific coast of North America from southwestern British Columbia to southern California and extends eastward into the deserts of the Southwestern United States and southward to central Mexico. The range of the northern spotted owl (S. o. caurina) subspecies includes the Pacific Northwest portion of the S. occidentalis range from British Columbia south to the Golden Gate strait, California. The US Fish and Wildlife Service has listed S. o. caurina as “threatened” under the Endangered Species Act since 1990.

The barred owl (S. varia), formerly native east of the Rocky Mountains (Mazur & James, 2000), has extended its range into the western US in the last 50–100 years and, from British Columbia to southern California, has become broadly sympatric with the northern spotted owl in the last 50 years. Barred and spotted owls hybridize and successfully backcross (Haig et al., 2004; Kelly & Forsman, 2004; Funk et al., 2007). Mitochondrial DNA sequencing has served as a valuable tool in ascertaining the maternal lineage of western birds, especially in potential hybrids (Zink, 1994; Haig et al., 2004; Barrowclough et al., 2005; Ruegg, 2008; Krosby & Rohwer, 2009; Williford et al., 2014).

Population-level studies of the genetics of S. occidentalis and S. varia have mainly used two mitochondrial markers, a partial control region sequence (Barrowclough, Gutierrez & Groth, 1999; Haig et al., 2004; Barrowclough et al., 2005) and cytochrome b (cyt b) (Haig, Mullins & Forsman, 2004), although a phylogeographic study of S. varia also utilized portions of ND6 and COIII (Barrowclough et al., 2011). The sequences of the complete genomes of the mitochondria of these two species will aid researchers in utilizing additional mitochondrial markers in population genetic studies of these owls.

It is well known that mitochondrial genes can transfer to the nuclear genome; such regions of the nuclear genome are sometimes called Numts (Lopez et al., 1994; Sorenson & Quinn, 1998). As a high-quality nuclear genome of S. o. caurina is available (Hanna et al., 2017), we were able to explore the incidence of Numts within the nucleus and investigate which mitochondrial genes have most often transferred. Furthermore, by assessing divergence between mitochondrial genes and their descendent Numts, we ascertained the likelihood of them posing problems for phylogenetic and other types of studies.

Methods

Results

We deposited the complete mitochondrial genome sequences of Strix occidentalis caurina sample CAS:ORN:98821 and Strix varia sample CNHM<USA-OH>:ORNITH: B41533 as NCBI GenBank Accessions MF431746 and MF431745, respectively. The lengths of the final S. o. caurina and S. varia mitochondrial genome assemblies were 19,889 nt and 18,975 nt, respectively. As for all typical avian mitochondrial genomes, they are circular and code for two rRNAs, 22 tRNAs, and 13 polypeptides (Figs. S1 and S2). The annotations produced by MITOS identified a 1 nt gap that split ND3, which is consistent with the untranslated nucleotide and translational frameshift seen in ND3 in some other bird species (Mindell, Sorenson & Dimcheff, 1998b), including owls (Strigiformes) (Fuchs et al., 2008).

Both the S. o. caurina and S. varia mitochondrial genomes contain a duplicated control region (Fig. 1). In both genomes, CR1 and CR2 each include a C-rich sequence near the 5′ end, the goose hairpin (Quinn & Wilson, 1993), which is identical across the two species and across CR1 and CR2. The S. o. caurina CR1 contains a 70 nt motif repeated 6.8 times near the 3′ end while CR2 includes two sets of tandem repeats near the 3′ end of the region, a 70 nt motif repeated four times followed by 9.5 repetitions of a 78 nt motif (Table 3).

The S. o. caurina CR1 and CR2 share a conserved central block of 1,222 nt with only two mismatches between CR1 and CR2 (Fig. 2). This conserved block includes 202 nt of the 5′ portion of the repetitive regions. The S. varia CR1 and CR2 share a conserved 1,041 nt central sequence stretch containing five mismatches. In CR1, this conserved block begins in the 3′ 57 nt of the CR1 repetitive region, but in CR2 it does not extend into the repetitive region. The 5′and 3′ regions surrounding the conserved central blocks of the control regions in both S. o. caurina and S. varia are more divergent from each other.

We obtained an alignment (88.37% identity) of 1,429 nt from the 5′ ends of the S. o. caurina and S. varia CR1 sequences, but it included fifteen gaps (Fig. 3). In contrast, the more 3′ repetitive sections of the S. o. caurina and S. varia CR1 sequences yielded an uninformative alignment with numerous, long gap regions. Similarly to CR1, the 5′ ends of the S. o. caurina and S. varia CR2 sequences aligned well (90.62% identity), yielding a 1,300 nt alignment that included four gaps. However, the alignment of the 3′ region of the CR2 sequences was uninformative with numerous, long gaps due to conflicts between the 78 nt motif repetitive regions of the two CR2 sequences. We found no evidence of mitochondrial pseudogenes in the control region sequences of either S. o. caurina or S. varia.

Across all of the 35 genes that were present in the previous, incomplete S. o. caurina assembly that was produced as a byproduct of the assembly of the S. o. caurina whole nuclear genome (Hanna et al., 2017), we only found one mismatch with our complete assembly, which occurred between the two ND1 sequences. This assembly improves upon the previous version by providing the complete sequences of ND6, tRNA^Pro, and the two control regions.

The S. o. caurina CR1 is 2,021 nt in length and the S. varia CR1 is 1,686 nt long. In both species, the 5′ end of CR1 borders tRNA^Thr and the 3′ end is adjacent to tRNA^Pro, then ND6, and then tRNA^Glu (Fig. 1). The initial 1,104 nt of the S. o. caurina CR1 are identical to a S. o. caurina partial control region sequence (GenBank Accession AY833630.1) (Barrowclough et al., 2005). All of the top 100 matches of the BLASTN searches of the S. o. caurina CR1 to NCBI-nt were to either S. occidentalis or S. varia control region sequences deposited by other researchers, as we expected.

CR2 follows tRNA^Glu and is 2,319 nt in length in S. o. caurina and 1,719 nt long in S.  varia. The initial 549 nt of the S. o. caurina CR2 matches the beginning of the D-loop sequence of an annotated complete genome of a Bubo bubo mitochondrion (GenBank Accession AB918148.1) (Hengjiu et al., 2016). One of the top 100 matches of the BLASTN searches of the S. o. caurina CR2 to NCBI-nt, which had the highest total score (2,177) and query coverage (96% versus 36–41% for the other matches) of the top 100 matches, was to a S. uralensis control region sequence (GenBank Accession AB743794.1) (Omote et al., 2013). The majority of the primers used by Omote et al. (2013) to PCR-amplify the control region in S. uralensis align within and around the S. o. caurina CR2. Four of the control-region-specific primers align to the middle of CR2 in our S. o. caurina sequence, which is identical to the middle of the S. o. caurina CR1 sequence. Perhaps most crucially, the primer L16728 aligns in the forward direction in tRNA^Glu such that it would amplify CR2, if present in the species.

As we mentioned in the methodology, our PCR-amplification of the S. o. caurina CR1 using primers that spanned from cyt b to ND6 yielded two products approximately 2,250 and 3,500 nt in length (Fig. S3). The sequences of these fragments were identical in the cyt b and ND6 portions as well as in the adjacent CR1 sections except when they entered the repetitive region at the 3′ end of CR1. We were only able to obtain sequence spanning the entirety of this repetitive region in the 2,250 nt fragment. This was largely due to the fact that the 3,500 nt fragment, in addition to the 70 nt motif repetitive section observed in the sequence of the 2,250 nt fragment, contained another repetitive region on the tRNA^Pro side of the 70 nt motif region with at least 13.1 copies of a 67 nt motif. We did not find any copies of tRNA^Pro or ND6 in the S. o. caurina nuclear genome, but we did find nuclear copies of cyt b and tRNA^Thr. With ND6 absent from the nuclear genome, PCR-amplification using primers in cyt b and ND6 should only generate mitochondrial genome fragments. Additionally, the cyt b and tRNA^Thr sequence in both the 2,250 nt and 3,500 nt fragments did not match the nuclear genome copies of these genes. In summary, we believe that both of these fragments were mitochondrial in origin and this evidence suggests that at least two different versions of the mitochondrial genome were present in this S. o. caurina individual.

The annotations of the mitochondrial genome sequences of the owls (Aves: Strigiformes) Tyto alba, Ninox novaeseelandiae, Strix leptogrammica, Glaucidium brodiei, and Asio flammeus indicate that those owls all share the same mitochondrial gene order as Gallus gallus, the ancestral avian mitochondrial gene order (Desjardins & Morais, 1990; Mindell, Sorenson & Dimcheff, 1998a; Haddrath & Baker, 2001) (Fig. 1). Our alignment of the S. leptogrammica mitochondrial genome to the mitochondrial genomes of other owls, including our S. o. caurina and S. varia assemblies, resulted in a poor, gap-filled alignment of the genes from the second half of the S. leptogrammica cyt b sequence through ND6 to tRNA^Phe. We could not obtain a reasonable alignment of the last 210 nt of the S. leptogrammica D-loop adjacent to the tRNA^Phe sequence to our S. o. caurina and S. varia assemblies or to the mitochondrial genomes of any of the other Strigiformes whose sequences we examined. Additionally, alignment of the S. leptogrammica mitochondrial genome with our S. o. caurina and S. varia assemblies yielded an ND5 alignment with seven gaps and numerous mismatches (85.60% and 84.82% identity to S. o. caurina and S. varia, respectively). Together, these results suggest that the S. leptogrammica sequence potentially contains significant errors in the sequences of the genes from ND5 through ND6 to tRNA^Phe.

We found 29,520 nt of Numt sequences in the draft S. o. caurina nuclear genome assembly spanning nine Numts (Table 4). The Numts ranged in length from 226–19,522 nt and had an average length of 3,280 nt. The Numts provided evidence of nuclear copies of all mitochondrial genes, except tRNA^Pro, ND6, and tRNA^Glu, the three genes between CR1 and CR2. Numt #9 (Table 4) aligns to both CR1 and CR2 with the alignments extending into the conserved block shared by the control regions. The portion of genome scaffold-294 aligned to CR2 for this Numt is 519 nt, whereas the length aligned to CR1 is 592 nt. As we could not be sure of which control region was incorporated into the nuclear genome, we have provided information for both alignments and derived the length of the Numt from the alignment to CR1 (Table 4).

Strix occidentalis caurina and S. varia display an average of 10.74% (8.68% uncorrected p-distance) base substitutions per site across the 2 rRNA genes and 13 polypeptide genes (the non-tRNA mitochondrial genes) (Table 5). The lowest number of base substitutions per site occurs within 16S and the highest within ATP8 (Table 5).

Discussion

Sequences of most mitochondrial genes can often be recovered from high-throughput short-read sequencing data if genome complexity is not too great. Algorithms using short-read data have more difficulty assembling low-complexity or repetitive regions due to an inability to span these regions. Thus, assembly of complete mitochondrial genome sequences can be more difficult when such genomes include regions of low-complexity. The sequence of the avian control region can both contain blocks of tandem repeats (Omote et al., 2013) and be duplicated (Eberhard, Wright & Bermingham, 2001; Haring et al., 2001). Moreover, the presence of multiple controls regions that are similar or identical, which has been observed in snakes (Kumazawa et al., 1996), can cause problems with assembly. In such situations, additional types of sequencing data that complement short-read data may be necessary in order to obtain an accurate and complete assembly of the mitochondrial genome. This proved to be the case in our study where the longer Sanger-derived sequence data were crucial in obtaining the complete sequence of the lengthy, repeat-rich control regions in S. o. caurina and S. varia. Although Brito (2005) and Barrowclough et al. (2011) inferred the presence of a duplicated control region structure in the mitochondrial genomes of at least three wood owl species, Strix aluco, S. uralensis, and S. varia, they did not sequence complete mitochondrial genomes. They likely deduced that a duplication was present from the appearance of multiple bands on agarose gels resulting from PCR-amplification of portions of the mitochondrial control region. Here we describe the first complete genome sequences of the mitochondrion of an owl (Aves: Strigiformes) with a duplicate control region.

The mitochondrial genomes of S. o. caurina and S. varia exhibit the novel avian gene order first described by Mindell, Sorenson & Dimcheff (1998a) for several bird orders, but not reported by them as present in the owl Otus asio. As mentioned above, this duplicated control region structure and novel gene order has previously been reported in the mitochondrial genome of S. varia (Barrowclough et al., 2011) and the congeners S. aluco and S. uralensis (Brito, 2005). The novel gene order was previously implied for S. occidentalis by the placement of primer N1 in tRNA^Thr by Barrowclough, Gutierrez & Groth (1999) to PCR-amplify the control region (CR1) fragment used in their study. Hull et al. (2010) also used the Barrowclough, Gutierrez & Groth (1999) N1 primer to PCR-amplify the control region in their study of S. nebulosa, so we can infer that the S. nebulosa mitochondrion also possesses the Mindell, Sorenson & Dimcheff (1998a) novel gene order. Notably, this mitochondrial gene order was not reported as present in S. leptogrammica (Liu, Zhou & Gu, 2014). However, our alignments of this mitochondrial genome to our S. o. caurina and S. varia sequences as well as the sequences of other owl mitochondrial genomes indicated problems with the S. leptogrammica sequence from cyt b through ND6 to tRNA^Phe. If we then leave aside the S. leptogrammica sequence, available evidence suggests that the novel gene order and duplicate control region structure is present across the genus Strix.

The primers developed by Barrowclough, Gutierrez & Groth (1999) to PCR-amplify a fragment of the control region (CR1) in S. occidentalis have been used extensively in additional genetic studies of owl species (Haig et al., 2004; Brito, 2005; Marthinsen et al., 2009; Hull et al., 2010; Barrowclough et al., 2011; Hausknecht et al., 2014). The Barrowclough, Gutierrez & Groth (1999) control region primers D16 (the most 3′ of their primers) and D20 (more 5′relative to primer D16) align to a region conserved between CR1 and CR2, although the length of the distance from the binding site of primer N1 in tRNA^Thr to the CR2 sites of primers D16 and D20 (3,742 nt and 3,392 nt, respectively, in our S. o. caurina assembly) likely reduces the probability of this second primer binding site causing problems in the PCR-amplification of the CR1 fragment.

The second control region appears intact, not degraded as found in some other avian taxa (Mindell, Sorenson & Dimcheff, 1998a). This gene order corresponds to the “Type D Duplicate CR gene order” of Gibb et al. (2007) and the “Duplicate CR gene order I” of Eberhard & Wright (2016). The goose-hairpin structure is typically found near the beginning of the control region in avian mitochondria (Marshall & Baker, 1997; Randi & Lucchini, 1998; Bensch & Härlid, 2000) and, in agreement with what we found, it appears in the beginning of the intact, duplicated control region sequences in the genomes of Amazona (Eberhard, Wright & Bermingham, 2001) and additional parrot mitochondria (Eberhard & Wright, 2016).

The lengths of the S. o. caurina CR1 and CR2 (2,021 nt and 2,319 nt, respectively) and of the S. varia CR1 and CR2 (1,686 nt and 1,719 nt , respectively) are all shorter than the length reported for the control regions of some species in the sister genus of owls, Bubo (Wink et al., 2009), which have lengths up to approximately 3,800 nt due to tandem repeats in the 3′ end of the control region (Omote et al., 2013). Similar tandem repeat blocks occur in the control regions of several other owl species in the family Strigidae (Xiao et al., 2006; Omote et al., 2013). The length of the tandem repeat motif unit is 78 nt in the 3′end of the control region sequences of Bubo blakistoni, Bubo virginianus, Strix uralensis (Omote et al., 2013), and Strix aluco (Xiao et al., 2006); 78 nt is also the length of the motif in the longest tandem repeat block in both the S. o. caurina and S. varia CR2 (Table 3).

As we previously mentioned, both S. uralensis and S. aluco exhibit a duplicated control region structure in their mitochondrial genomes (Brito, 2005). Neither Omote et al. (2013) nor Xiao et al. (2006) report the presence of a duplicated control region structure in either S. uralensis or S. aluco, respectively, in their discussions of the repetitive content of the control region sequences of these two species. It is not overtly clear from their methodologies which control region they sequenced. The precise primer combinations used for the PCR-amplification and sequencing of the control region of the Bubo species and S. uralensis are not provided by Omote et al. (2013), but mapping the primer sequences used by the researchers to our S. o. caurina genome suggests that, if the structure of the S. uralensis mitochondrial genome shares that of S. o. caurina, they likely sequenced CR2 in at least S. uralensis and in the Bubo species if a CR2 was present. We are unsure how placement of primers in cyt b and 12S, as reported in the methodology of Xiao et al. (2006) could PCR-amplify a single control region sequence for S. aluco, given the duplicated control region structure (Brito, 2005).

The duplicated control region structure is unreported in Strigiformes outside of Strix, but we infer that it is also likely present in Bubo due to the positioning of primers used for PCR-amplification of the control region in previous studies (Marthinsen et al., 2009; Omote et al., 2013). If also present in Bubo, then the duplicate control region structure appears to have arisen in the common ancestor of Strix and Bubo, but a proper phylogenetic test of this hypothesis with increased taxon sampling is warranted. Further work on the structure of control region sequences in Strix and related taxa is needed to elucidate the pattern of evolution of this region across the Strigidae phylogeny.

Although inconclusive and warranting further investigation, our evidence for two versions of the 3′ repetitive region of CR1 suggests that mitochondrial heteroplasmy is present in this S. o. caurina individual. Mitochondrial heteroplasmy due to tandem repeat variability in the control region has been shown to occur in other bird species (Berg, Moum & Johansen, 1995; Mundy, Winchell & Woodruff, 1996). Previous work has suggested that the most likely mechanism by which the gain and loss of such tandem repeat elements occurs in the mitochondrial control region is that the repetitive region forms a stable, single-stranded secondary structure and there is slippage during replication (Levinson & Gutman, 1987; Wilkinson & Chapman, 1991; Fumagalli et al., 1996; Faber & Stepien, 1998). Greater numbers of repeats may improve the stability of the secondary structure (Faber & Stepien, 1998). Utilizing sequence from the 3′ region of CR1 for population genetic study of S. o. caurina is not likely to be useful due to the variability (in terms of the number of copies of the tandem repeat motifs in this region) that is potentially present within a single individual.

The 29,520 nt of Numt sequence in the draft S. o. caurina nuclear genome assembly is more than triple the 8,869 nt of Numt sequence found in a Gallus gallus draft nuclear genome assembly (Pereira & Baker, 2004). The 3,280 nt average Numt size exceeds the average size in all of the eukaryotic genomes examined by Richly & Leister (2004). There are markedly fewer control region Numts in the S. o. caurina draft genome assembly than found in a Gallus gallus draft genome assembly (Pereira & Baker, 2004). We only found one control region Numt (Table 4). Indeed the longest Numt, Numt #1, extends through almost the entire mitochondrial genome sequence including from tRNA^Phe through tRNA^Thr, immediately adjacent to, but ending at CR1. The percentage identity of the nuclear pseudogenes with the true mitochondrial genes ranges from 77.5–87.81%, so care must be taken to insure that Numts are not PCR-amplified in place of mitochondrial gene sequences. As the control region is the mitochondrial segment that has been used most often in studies of the population genetics of Strix species (Barrowclough, Gutierrez & Groth, 1999; Haig et al., 2004; Barrowclough et al., 2005; Brito, 2005; Hull et al., 2010; Barrowclough et al., 2011), it is encouraging that only one, short Numt included CR1 or CR2 (Table 4). As long as researchers PCR-amplify sequences that span beyond the 592 nt Numt #9, they should have confidence in amplifying the true mitochondrial control regions.

The average pairwise sequence divergence between S. occidentalis and S. varia has been previously reported as 13.9% for a 524 nt section of CR1 (Haig et al., 2004). This exceeds the weighted average of 10.74% (8.68% uncorrected p-distance) that we calculated across the non-tRNA mitochondrial genes (Table 5), which is unsurprising as the control region is known to be rapidly evolving in birds (Quinn & Wilson, 1993). However, the pairwise sequence divergence between S. occidentalis and S. varia appears higher in ND3, ND4L, ND6, and ATP8 (Table 5) than in the 524 nt CR1 portion. Hull et al. (2010) found uncorrected p-distances of 13.73–13.93% for ND2 and 14.58–14.81% for the control region between S. nebulosa and S. occidentalis. Wink & Heidrich (1999) calculated uncorrected p-distances of 8.15–11.72% for cyt b sequence pairwise comparisons of six Strix species (S. aluco, S. butleri, S. nebulosa, S. rufipes, S. uralensis, and S. woodfordii). Of these six Strix species, only S. aluco x S. uralensis are known to hybridize (McCarthy, 2006), albeit in captivity (Scherzinger, 1982), and their cyt b pairwise divergence was 8.15%, the lowest of those calculated by Wink & Heidrich (1999) between Strix species. Our cyt b uncorrected p-distance value between S. o. caurin a and S. varia was 9.21%, which is also on the lower end of the range of the Wink & Heidrich (1999) Strix interspecific divergences. Overall, however, high levels of interspecific pairwise divergence of mitochondrial DNA seem to be typical for the genus Strix, even for species able to hybridize. We anticipate that these whole mitochondrial genome resources will be useful to those with an interest in developing new mitochondrial markers to study the genetics of S. o. caurina, S. varia, and related taxa.

Supplemental Information