Stabilizing selection and mitochondrial heteroplasmy in the Canada lynx (Lynx canadensis)
Abstract
Mitochondrial DNA is commonly used in population genetic studies to investigate spatial structure, intraspecific variation, and phylogenetic relationships. The control region is the most rapidly evolving and largest non-coding region, but its analysis can be complicated by heteroplasmic signals of genome duplication in many mammals, including felids. Here, we describe the presence of heteroplasmy in the control region of Canada lynx (Lynx canadensis) through intra-individual sequence variation. Our results demonstrate multiple haplotypes of varying length in each lynx, resulting from different copy numbers of the repetitive sequence RS-2 and suggest possible heteroplasmic single nucleotide polymorphisms in both repetitive sequences RS-2 and RS-3. Intra-individual variation was only observed in the repetitive sequences while inter-individual variation was detected in the flanking regions outside of the repetitive sequences, indicating that heteroplasmic mutations are restricted to these repeat regions. Although each lynx displayed multiple haplotypes of varying length, we found the most common variant contained three complete copies of the RS-2 repeat unit, suggesting copy number is regulated by stabilizing selection. While genome duplication offers potential for increased diversity, heteroplasmy may lead to a selective advantage or detriment in the face of mitochondrial function and disease, which could have significant implications for wildlife populations experiencing decline (e.g., bottlenecks) as a result of habitat modification or climate change.
Introduction
Mitochondrial DNA is widely used in molecular studies, particularly in animals, due to its unique features including maternal inheritance without recombination, compact size, and high density in tissues (Avise 1994; Jae-Heup et al. 2001; Schwartz et al. 2007). One of the most useful and commonly studied regions of the mammalian mitochondrial genome in population genetics is the control region or D-loop, as it is the largest non-coding region with a high mutation rate (Fig. 1; Hoelzel et al. 1991; Brown et al. 1993; Jae-Heup et al. 2001). However, the benefits of using the control region when investigating population structure, intraspecific variation, and phylogenetic relationships (Avise 1994) can be complicated when two or more mitochondrial genomes exist in the same cell leading to variation within individuals, i.e., heteroplasmy (Wilkinson and Chapman 1991).
Fig. 1.
Mitochondrial heteroplasmy appears to be widespread across vertebrates and more common in mammals than previously thought (Hoelzel et al. 1994; Fumagalli et al. 1996; Wilkinson et al. 1997; Rensch et al. 2016; Päckert et al. 2019). Among the Felidae, mitochondrial heteroplasmy has been documented in cheetah (Acinonyx jubatus) (Freeman et al. 2001), five species of Panthera (Jae-Heup et al. 2001), Eurasian lynx (Lynx lynx) (Sindičić et al. 2012), and bobcat (Lynx rufus) (Croteau et al. 2012). Two of these studies reported heteroplasmic copy number variation within individuals (Jae-Heup et al. 2001; Sindičić et al. 2012), whereas others chose to exclude repetitive regions to avoid problems with heteroplasmy or sequence quality (e.g., Hellborg et al. 2002; Gugolz et al. 2008; Croteau et al. 2012; Anderson et al. 2015).
The control region in felids is generally split into three domains: the left domain, which includes the first hyper-variable sequence (HVS-1) and repetitive sequence 2 (RS-2); the right domain, which includes the second hyper-variable sequence (HVS-2) and repetitive sequence 3 (RS-3); and the central conserved region (CCR), situated between the two, which has relatively low variation compared to the others (Fig. 2A; Jae-Heup et al. 2001). RS-2 comprises an array of ∼80 bp tandem repetitive units that have been identified in several mammals (Fig. 2B; Wilkinson and Chapman 1991; Hoelzel et al. 1994), including domestic cat (Felis catus) (Lopez et al. 1996) and Eurasian lynx (Hellborg et al. 2002). The number of complete copies of these repetitive units can range from 1 to 4 across several felid species including cheetah (Freeman et al. 2001) and great cats in the genus Panthera (Jae-Heup et al. 2001), though strangely, the RS-2 repetitive sequence is either absent or highly reduced in the clouded leopard (Neofelis nebulosa) (Wu et al. 2007). Jae-Heup et al. (2001) also found single nucleotide polymorphisms (SNPs) in the RS-2 repeat unit within individual cats, suggesting that repeat copy number is not the only source of heteroplasmic variation.
Fig. 2.
RS-3 is a more complex repetitive sequence compared to RS-2 and is found in most carnivore lineages (Hoelzel et al. 1994; Ketmaier and Bernardini 2005). It varies in size and number of repeats, but among felids, it has a repeat unit of 6–10 bp with a core motif of TACAC or TACGC followed by various suffixes (Fig. 2C; Wu et al. 2007; Sindičić et al. 2012). Despite high levels of variation, this sequence has thus far proven to be uninformative for phylogenetic analyses (Jae-Heup et al. 2001; Ketmaier and Bernardini 2005).
These repetitions in RS-2 and RS-3 are thought to be the result of misalignment and slippage due to the D-loop’s involvement in replication and are expected to be under stabilizing selection (Buroker et al. 1990; Hoelzel et al. 1994). For example, multiple repeats may be selected for if they aid in binding proteins required for replication or compensate for deleterious mutations, but a larger genome will also result in slower replication and likely be selected against (Wilkinson et al. 1997; Jae-Heup et al. 2001; Clark et al. 2012). Additionally, the proportions of heteroplasmic mutations vary depending on the individual, age, and type of tissue, and may be influenced by tissue-specific selection for new variants or as a by-product of accumulated oxidative damage (Matthews et al. 1994; Michikawa et al. 1999; Jae-Heup et al. 2001).
Mutations in the mitochondrial genome can cause severe disease in humans (Kopinski et al. 2019; Jackson et al. 2020), inbred lab mice (Lechuga-Vieco et al. 2022), domestic dogs (Li et al. 2006; Baiker et al. 2009), and at least one domestic cat (Dell'Era et al. 2021). Two mutations in the tRNA-Leu gene were linked to selective symmetrical necrotizing encephalopathy in a house cat resulting in a progressive reduction in cognitive and motile function (Dell'Era et al. 2021). While this disease is not a direct result of heteroplasmy, the two causative mutations do present in a heteroplasmic fashion (Dell'Era et al. 2021), which in turn, can affect the clinical presentation of the disease and make diagnosis challenging (Jackson et al. 2020). Clinical disease in humans can either vary in severity according to the proportion of mutant DNA in the cell or present suddenly when a threshold of 60% is reached (Jackson et al. 2020). Furthermore, wild felid species with low genetic diversity can suffer inbreeding depression from abnormal sperm morphology and impaired sperm metabolic function (Pukazhenthi et al. 2006; Terrell et al. 2010). Since impairment of mitochondrial metabolism or defects in the mitochondrial genome reduce sperm motility in humans (Amaral et al. 2013; Irigoyen et al. 2022), there may also be a vital link between mitochondrial heteroplasmy and evolutionary fitness in wild felids.
Here, we aim to investigate intra-individual variation in the mitochondrial control region of Canada lynx (Lynx canadensis) which has not been described despite previous studies quantifying population-level variation (e.g., Rueness et al. 2003; Prentice et al. 2019). The Canada lynx is part of a keystone predator–prey cycle in the boreal forest and undergoes dramatic fluctuations in population size (Krebs et al. 2018). Thus, Canada lynx represents a valuable opportunity to study how mitochondrial heteroplasmy is maintained or influences fitness in small populations which will eventually assist in conservation and management decisions. With this in mind, we provide the first characterizations of heteroplasmic variation in copy number of the repetitive sequence RS-2 and high inter-individual variation in copy number of repetitive sequence RS-3 in this species and discuss the implications for stabilizing selection on genome evolution.
Materials and methods
Amplifying the left domain
We extracted DNA from 24 Canada lynx muscle samples from Yukon, Canada using the DNeasy Blood and Tissue kit (Qiagen, Valencia, USA) and amplified the left domain of the control region with primers mtu-F (5′-CTTTGGTCTTGTAAACCAAAAAA-3′) and R3-R (5′-TAAGAACCAGATGCCAGGTA-3′) previously used by Rueness et al. (2003) and Prentice et al. (2019) to target an ∼750 bp product containing HVS-1, RS-2, and some of the CCR (Fig. 2A). Reaction volumes were 50 µL and contained 2 µL of template DNA, 0.4 µmol/L of each primer, 25 µL of HotStart Taq 2X Mastermix (New England Biolabs [NEB], Ipswich, USA), and a final concentration of 2.5 mmol/L of MgCl2. PCR cycling conditions were 95 °C/5 min, (95 °C/30 s, 56 °C/30 s, and 72 °C/1 min) × 30 cycles, and a final extension of 72 °C/2 min. Products were visualized on a 2% agarose gel, resulting in multiple evenly-spaced bands of varying lengths within individuals, similar to Foran et al. (1997) who saw the same banding pattern in one cougar (Puma concolor) and one Canada lynx using conserved mammalian primers.
We then designed new primers Lc-mtd-F (5′-ACCGCCTCCTCAAATGAAGA-3′) and Lc-mtd-R (5′-TGTGTGATCATGGGCTGATT-3′) encompassing the previous pair using the Canada lynx mitochondrial reference genome (GenBank NC_028313; Li et al. 2016), which predicted a product size of 862 bp (Fig. 2A). The PCR reaction mixture and cycling conditions for the new primers were identical to the former and agarose gel electrophoresis revealed the same pattern with multiple evenly-spaced bands of varying lengths (Fig. 3A).
Fig. 3.
Isolating individual length variants
We isolated length variants amplified with the external primers (Lc-mtd-F and Lc-mtd-R) using a band-stab approach (Bjourson and Cooper 1992; Drew and Brindley 1997; Wilton et al. 1997). A 2% agarose gel pre-stained with SYBRSafe (Invitrogen, Waltham, USA) was loaded with 25 µL of each PCR product mixed with 5 µL of 6X purple loading dye (NEB). Gels were electrophoresed for 45 min at 120 V, then 135 min at 85 V, and finally 30 min at 65 V for optimal separation and resolution of bands. We stabbed several times across the width of each band using a 200 µL pipet tip to remove agarose plugs containing PCR products. Plugs were transferred to 50 µL reaction mixtures as described above but with primers mtu-F and R3-R and 3.5 mmol/L final concentration of MgCl2. Agarose plugs were allowed to rest for at least 30 s in the new PCR tube, before being swirled using the pipette tip and pipetted up and down three times (25 µL volume). Products were reamplified using the same cycling conditions as the initial PCR and the band-stab procedure was repeated 2–3 times as necessary to achieve pure isolates of each product (Fig. 3B).
Left domain fragments from one individual were also isolated by PCR cloning. Original products amplified from mtDNA using primers mtu-F and R3-R were purified using the Monarch PCR & DNA Cleanup Kit (NEB) to remove leftover PCR components. Cleaned PCR products were cloned into the pCR2.1-TOPO vector and transformed into OneShot TOP10 competent E. coli cells following the manufacturer’s protocol for chemical transformation with the TOPO TA Cloning Kit (Invitrogen). Transformants were screened with X-Gal on LB plates containing 100 µg/mL Ampicillin. White colonies were picked with a 200 µL pipet tip, which was then swirled into a 50 µL PCR reaction mixture. Reaction components followed initial PCR except using 0.4 µmol/L of the M13 universal primer pair. Cycling conditions for cell lysis and initial denaturation were 95 °C for 10 min, followed by 35 cycles of 95 °C for 30 s, 50 °C for 30 s, and 72 °C for 1 min, before a final extension at 72 °C for 10 min. Fragment isolation was confirmed by agarose gel electrophoresis as described above.
Amplifying the right domain
We amplified the right domain of the control region including the RS-3 repeat array using a new primer Lc-RS3-F (5′-GGAAGCTTAATCACCTGGTC-3′), which we designed from the same lynx reference genome as the left domain, paired with F20-R (5′-GACTCATCTAGGCATTTTCAG-3′; Wu et al. 2007) to target a 1078 bp product that partially overlapped with our sequence from the left domain (Fig. 2A). PCR reaction volumes were 25 µL consisting of 2 µL of template DNA, 0.8 µmol/L of each primer, and 12.5 µL of HotStart Taq 2X Mastermix (NEB) containing a final concentration of 1.5 mmol/L MgCl2, with the same cycling conditions as the left domain. Visualization of products on a 2% agarose gel showed a single band per individual so we proceeded directly with sequencing.
Sanger sequencing
All products from the left and right domains (including products separated by band-stab or cloning and original PCR products that were not isolated) were Sanger sequenced in both directions by the University of Alberta’s Molecular Biology Services Unit and a consensus sequence was assembled for each product in Geneious 10.2 (Kearse et al. 2012). All consensus sequences were aligned to compare length and nucleotide polymorphisms and uploaded to GenBank (left domain OR338890-OR338894, right domain OR338895-OR338915). To ensure our sequences represented the mitochondrial genome and not mitochondrial DNA transposed into the nuclear genome (NUMT), we aligned our 724 bp sequence from the left domain (including 3.5 copies of the RS-2 repeat; OR338893) to the Canada lynx nuclear (GCF_007474595) and mitochondrial (NC_028313) reference genomes using BLAST in GenBank and BWA-mem (Li and Durbin 2009) with default parameters.
Nanopore sequencing
We constructed six barcoded sequencing libraries with original PCR products (without isolation) from the left domain of six individuals by following the Native Barcoding Kit 24 V14 protocol (Oxford Nanopore Technologies, Oxford, UK) with 400 ng purified amplicons from each individual. Briefly, amplicons were A-tailed and ligated to unique dT-tailed native barcode adapters and bead cleaned with a 1:1 DNA to bead ratio using PCR DX beads (Aline Biosciences). Barcoded samples were pooled and ligated to the Oxford Nanopore sequencing adapter. A final size selection using 0.6:1 bead to library ratio selected against smaller molecules. In total, 100 ng of final library products were sequenced on a PromethION P48 instrument using the R10 M flow cell for 96 h by the BC Cancer Research Institute’s Genome Sciences Centre. DNase I (Invitrogen) nuclease flush was performed after 24 h by reloading the flow cell with the same library mix. Reads were base-called and demultiplexed using dorado 0.3.4 (Oxford Nanopore 2023), then aligned to the complete 80 bp repeat motif using minimap 2.2.26 (Li 2018) with the -x ont setting. We classified each read according to the number of times it aligned to the repeat motif to determine the number of repeat copies in each read, then counted the read types in each category within each individual. We used an Analysis of variance with Tukey’s honestly significant difference post hoc test to determine if the percentage of sequencing reads within each individual differed among categories according to the number of RS-2 repeat copies. Similarly, we aligned these sequences to the Canada lynx nuclear (GCF_007474595) and mitochondrial (NC_028313) reference genomes using BWA-mem (Li and Durbin 2009).
Results
Left domain
Through agarose gel electrophoresis of products from the left domain, we counted a mean of four bands per individual (maximum of seven, Fig. 3A) after multiple rounds of band-stab and reamplification. We also obtained four bands in the same size range after cloning fragments from one individual. We successfully sequenced 17 bands from five unique individuals representing replicates of the five most common band sizes across both isolation methods. All lynx showed evidence of heteroplasmy with more than one band size. Furthermore, Sanger sequencing showed that differences in sequence length were strictly the result of variable copy numbers of the 80 bp RS-2 repeat unit, ranging from 0–4 complete copies resulting in total lengths of 484, 564, 644, 724, and 804 bp after trimming primer sequences (Fig. 4). Each repeat loosely matched the 80 bp motif described by Lopez et al. (1996) and each sequence contained the 44 bp partial repeat at the 3′ end of the array. While we did not detect any SNPs among the RS-2 repeat units from isolated sequences within the same individual, we saw potential evidence of heteroplasmy in some individuals in the form of “heterozygous” chromatogram peaks within the RS-2 repeat units when sequencing the original PCR amplicons (i.e., without isolation). Furthermore, we observed inter-individual SNPs outside of the RS-2 repeat array in clones, band-stab isolates, and original products without isolation, though three SNPs found among two clones were not present in band-stab isolates or the original product from the same individual.
Fig. 4.
When BLASTing Sanger data, only the first 364 bp from the 5′ end (up to the start of RS-2) of our 724 bp sequence matched with chromosome D2 from the lynx nuclear reference genome (the remaining 360 bp did not align) and this match was only 80% identical as our sequence included 61 nucleotide variants and three insertions of 1, 2, and 7 bp each compared to D2. In contrast, the full 724 bp of our left domain sequence was 98% identical with the lynx mitochondrial reference genome. Alignments with BWA were similar; the 724 bp sequence fully aligned to the mitochondrial genome, but was only partially aligned to the nuclear genome.
The number of Nanopore sequencing reads ranged from 6.5 M to 10.6 M per individual and the percentage of reads within each individual varied among categories with different copy numbers of the RS-2 repeat unit. Reads with three complete copies were significantly more frequent 38.7 ± 6.9% (95% CI) than all other read categories (F6,35 = 60.2, P < 0.0001, Fig. 5). Reads with two and four copies were significantly more frequent than reads containing zero, one, five, and six copies (all P < 0.001), which did not differ from each other. Reads with zero complete repeats still had the 44 bp half repeat. Nanopore reads only partially aligned to the nuclear reference genome with 36% of reads appearing chimeric, while 100% of reads completely aligned with the mitochondrial genome. These findings imply that our sequences were from the mitogenome and not NUMTs.
Fig. 5.
Right domain
We sequenced 21 unique haplotypes in the right domain among 24 individual lynx. Length polymorphisms among individuals resulted in haplotypes varying from 908–1010 bp due to variable copy numbers of the RS-3 repetitive sequence, proving this region to be more variable among individuals than RS-2. Unlike RS-2, length polymorphism within individuals was not observed during gel electrophoresis, and there was no evidence of multiple fragments of different sizes during sequencing as forward and reverse sequences aligned perfectly. However, we found the same SNP occurring in one or two repeat units within two individuals that possessed a “heterozygous” peak in both their forward and reverse chromatograms (similar to RS-2), suggesting heteroplasmic mutations in this region.
The pattern of the RS-3 repeat array can be described in three sections that were composed of the repeat motifs TACAC and TACGC (Fig. 2C). The first section began with a single occurrence of the motif TACGC with a –G suffix followed by a variable number of copies of the repeat motif TACAC with a –G suffix until this changed to a –GTA suffix. This entire series then repeated itself (minus the TACGC start) multiple times over and continued with the same pattern of a variable number of TACAC(G) repeats followed by TACAC(GTA), until the second section which was a short string of three TACAC repeats always ending with –A, –G, and –ACG consecutively. Finally, the third section began with a variable number of TACAC(G) repeats followed by a variable-sized string of TACAC(GTA) repeats until the end of the array. Interestingly, TACAC was the dominant motif with a mean occurrence of 97.8 ± 0.5% (SD) among all repeats in each sequence while TACGC only occurred once in each of 23 individuals (though one lynx had two copies) and was consistently the first occurrence of a core motif from the 5′ end. The two most common suffixes for the repeat motifs were –G and –GTA with occurrences of 53.8 ± 5.7% and 41.0 ± 5.4%, respectively. The remaining suffixes occurred much less frequently, including –A (2.2 ± 0.4%), –ACG (2.5 ± 1.4%), –GTATACG (0.1 ± 0.5%), –RTA (0.3 ± 1.0%), and –GTG (0.1 ± 0.4%). Of these, –A would sometimes replace –G, while the rest would occur in lieu of –GTA, and every lynx had at least one occurrence of repeats ending with –A and –ACG.
Discussion
The control region is commonly sequenced for population and phylogenetic studies, but in many mammals, particularly felids, variable copy numbers of the RS-2 repeat unit can indicate heteroplasmy and distort sequence quality (Hoelzel et al. 1991; Brown et al. 1993; Jae-Heup et al. 2001). Although an electrophoretic signature of heteroplasmy was reported in Canada lynx (Foran et al. 1997), this phenomenon has never been confirmed by genetic sequence data, nor has it been reported in subsequent molecular studies (Rueness et al. 2003; Prentice et al. 2019). Hence, our study is the first to sequence isolates with different copy numbers and describe this phenomenon in this species. We were able to isolate and sequence five length polymorphisms from the left domain corresponding to the copy number of RS-2 repeat units that can exist simultaneously within individual lynx. Furthermore, we documented high haplotype diversity in the right domain through variability in the RS-3 repeat array among individuals.
We easily separated and counted four or five bands per individual through initial gel electrophoresis, but band-stab PCR allowed us to detect up to seven different band sizes from the left domain. The bands in the middle of the size range (corresponding to two and three complete repeats) were the brightest and most common while the longest band was only detected in one individual, and we were unable to isolate this product. Furthermore, we found that sequencing original PCR products (without band-stab or cloning) consistently resulted in sequences with three complete repeat units across all samples, suggesting that three copies is the most common number of repeats within and among individuals (Fig. 2B). Nanopore sequencing confirmed that amplicons with three complete repeat copies were the most frequent among individuals, and that frequency declined symmetrically as the number of repeat copies differed from three (Fig. 5). Our success at consistently sequencing three complete repeat copies from the original product also potentially explains why several previous studies never reported variation in repeat units or evidence of heteroplasmy (e.g., Rueness et al. 2003; Prentice et al. 2019). The pattern in our observed copy number frequencies is similar to Eurasian lynx which had two and three complete RS-2 repeats in 23% and 74% of individuals, respectively (Sindičić et al. 2012). Our pattern of declining repeat frequency from a mode of three supports the idea that the copy number of repeat units is under stabilizing selection (Buroker et al. 1990; Hoelzel et al. 1994; Lopez et al. 1996) and suggests that three complete repeats is the optimal copy number to balance replication efficacy and efficiency.
While we did not find SNPs within individuals when sequencing isolated fragments, we did find evidence of polymorphisms in the form of “heterozygous” chromatogram peaks when sequencing original PCR products without isolation, possibly as a result of SNPs among different mitochondrial genomes. Polymorphisms within the RS-2 repeat units are not uncommon and have been found in several species (Freeman et al. 2001; Jae-Heup et al. 2001; see raw data from Rueness et al. 2003), supporting the notion that heteroplasmic SNPs are possible in this region. Oddly, Sindičič et al. (2012) did not find any SNPs in the RS-2 repeat units of the closely related Eurasian lynx, though this may be a consequence of low genetic diversity across their sampling area. All inter-individual SNPs occurring outside of the RS-2 array were consistent across multiple band-stab isolates from the same individual, meaning heteroplasmic mutations seem to be minimal outside of the RS-2 region. However, two of the four clones introduced apparent SNPs outside of the RS-2 region, which were not present in band-stab or original PCR sequences for the same individual. While these SNPs are likely artifacts of the cloning procedure and not heteroplasmic mutations, we make note of any possibility of mitochondrial mutations as they can lead to serious mitochondrial disease in domestic cats (Dell'Era et al. 2021), domesticated dogs (Li et al. 2006; Baiker et al. 2009), inbred lab mice (Lechuga-Vieco et al. 2022), and humans (Kopinski et al. 2019; Jackson et al. 2020).
While mitogenome duplication has been described as mitochondrial heteroplasmy in Eurasian lynx and bobcat (Croteau et al. 2012; Sindičić et al. 2012), research on other species has suggested that intra-individual mitogenome sequence variants may also be attributed to NUMTs (Lopez et al. 1994; Antunes et al. 2007; Calabrese et al. 2017). Several NUMTs that contain the control region have been discovered in the domestic cat genome including on chromosome D2 (Lopez et al. 1994; Antunes et al. 2007), leading to the possibility that NUMTs may also exist in the Canada lynx genome. When comparing to the Canada lynx reference genomes, our sequences aligned better with the mitochondrial genome than the nuclear genome providing confidence to reject NUMTs in our data. However, we cautiously point out that Canada lynx likely has NUMTs on more than one chromosome, as we also found several partial alignments between the mitochondrial and nuclear reference genomes themselves without including our sequence data.
While beyond the scope of this paper, additional approaches for elucidating heteroplasmy versus NUMTs include isolating mitochondria for sequencing, long-read sequencing to capture large overlapping portions of the mitogenome to exclude the possibility of NUMTs (Patterson et al. 2023), or filtering tools that are designed to remove NUMTs from mitochondrial data, although, generally these programs require the usage of protein-coding genes, which is not the case for our data (Andújar et al. 2021; Graham et al. 2021; Porter and Hajibabaei 2021). Ultimately, we cannot entirely exclude the possibility of NUMTs in our data, but regardless of whether mitochondrial heteroplasmy or NUMTs are the cause of genome duplication in Canada lynx, the potential for increased genetic diversity and influence of stabilizing selection remains.
The RS-3 repeat array in Canada lynx generally matches the patterns seen in many other felid species. The most dominant core motif was TACAC, which appears to be unanimous among felines. Furthermore, like tiger (Panthera tigris), leopard (P. pardus), clouded leopard, and bobcat, Canada lynx also showed the rarer core motif TACGC (Jae-Heup et al. 2001; Sindičič et al. 2012). Interestingly, unlike the aforementioned species, the placement and occurrence rate of TACGC was remarkably consistent among individual lynx, appearing exclusively as the first motif from the 5′ end in all but one sample. Canada lynx also showed similar fundamental patterns in the suffixes at the end of the repeat motifs but still differed from even closely related congeners (bobcat and Eurasian lynx) upon closer inspection. Canada lynx and Eurasian lynx both showed –G as the most common suffix, but Canada lynx had far more –GTA repeats (41% occurrence), which usually appeared only at the 3′ end in Eurasian lynx (Sindičić et al. 2012). The most common suffix in bobcat was –ACG which was completely absent in Eurasian lynx and generally found only once in the second section for Canada lynx. Both bobcat and Eurasian lynx also had occurrences of –ACGTA which were not found at all in Canada lynx. Hence, the overall sequence repeat pattern of RS-3 appears to be unique in Canada lynx, as it is among all other felids investigated so far. While the full haplotype sequence of the RS-3 array was previously deemed uninformative for phylogenetic analysis (Jae-Heup et al. 2001; Ketmaier and Bernardini 2005), the general pattern of suffix repeats may be indicative of historical speciation among felids. One example could be the –ACG suffix commonly found in other felids, including bobcat, but is highly reduced or absent in other members of genus Lynx (Sindičić et al. 2012). Previous phylogenetic analyses have suggested that bobcat is the most basal member of the Lynx lineage (Johnson et al. 2004; Li et al. 2016) therefore, a reduction of the –ACG suffix could be an indicator of speciation within Lynx.
Here, we document heteroplasmy in the Canada lynx mitochondrial genome by isolating and sequencing products with heterogeneous repeat copy numbers within individuals, thereby adding Canada lynx to the list of felids exhibiting genome duplication. Notably, this research is the first description of intra-individual control region variation in this species, as this phenomenon has not been reported in previous population studies that have successfully sequenced across the left domain (Rueness et al. 2003; Prentice et al. 2019). Our results also demonstrate the efficacy of the band-stab method for isolating products from gels separated by as little as 80 bp. While the repeat sequence RS-3 has previously been deemed uninformative for phylogenetic analyses (Jae-Heup et al. 2001; Ketmaier and Bernardini 2005), sequencing different RS-2 length polymorphisms within individuals, particularly the first complete repeat from the 3′ end, may be of value to phylogenetic or population genetic analyses (Wilkinson and Chapman 1991; Jae-Heup et al. 2001). Ultimately, mitochondrial heteroplasmy provides evolutionary potential for increased genomic diversity which appears to be selected for among felids, possibly under stabilizing selection as a trade-off between increased transcription and slower replication (Wilkinson et al. 1997; Jae-Heup et al. 2001; Clark et al. 2012). Insight into the evolutionary advantage of this phenomenon will likely require a greater understanding of feline mitochondrial function and disease, though currently, there is little evidence for severe mitochondrial disease other than a singular case of selective symmetrical necrotizing encephalopathy in domestic cat (Dell'Era et al. 2021) and a potential link with teratospermia in wild felids (Pukazhenthi et al. 2006; Terrell et al. 2010). Perhaps this is because most heteroplasmic mutations in felids are restricted to non-coding regions, and thus have little detrimental effect. Further studies might investigate how individual factors known to influence heteroplasmy (i.e., age or tissue type) may affect the amplification and degree of heteroplasmy. Alongside this, a wider study documenting the extent of heteroplasmy across the species’ range may prove useful in further understanding this molecular phenomenon. Finally, we caution future studies to be aware of this genome duplication and to take care when inspecting sequence data.
Acknowledgements
The authors thank the First Nations, fur trappers, and government biologists, especially Caitlin Willier, Rudy Boonstra, Yasmine Majchrzak, and Tom Jung, who collected the samples used in this study, Joseph Nowlan, Scott Britney, and Ashton Sies for help and general commiseration in the lab, and three anonymous reviewers for helpful comments on earlier versions of the manuscript. Sanger sequencing was performed by the Molecular Biology Services Unit at the University of Alberta and Nanopore sequencing by Canada’s Michael Smith Genome Sciences Centre at the BC Cancer Research Institute.
References
Amaral A., Lourenço B., Marques M., Ramalho-Santos J. 2013. Mitochondria functionality and sperm quality. Reproduction, 146(5): R163–R174.
Anderson C.S., Prange S., Gibbs H.L. 2015. Origin and genetic structure of a recovering bobcat (Lynx rufus) population. Can. J. Zool. 93(11): 889–899.
Andújar C., Creedy T.J., Arribas P., López H., Salces-Castellano A., Pérez-Delgado A.J., et al. 2021. Validated removal of nuclear pseudogenes and sequencing artefacts from mitochondrial metabarcode data. Mol. Ecol. Resour. 21(1): 1772–1787.
Antunes A., Pontius J., Ramos M.J., O'Brien S.J., Johnson W.E. 2007. Mitochondrial introgressions into the nuclear genome of the domestic cat. J. Hered. 98(5): 414–420.
Avise J.C. 1994. Molecular markers, natural history, and evolution, Chapman and Hall, New York.
Baiker K., Hofmann S., Fischer A., Gödde T., Medl S., Schmahl W., et al. 2009. Leigh-like subacute necrotising encephalopathy in Yorkshire terriers: neuropathological characterisation, respiratory chain activities and mitochondrial DNA. Acta Neuropathol. 118(5): 697–709.
Bjourson A.J., Cooper J.E. 1992. Band-stab PCR: a simple technique for the purification of individual PCR products. Nucleic. Acids. Res. 20(17): 4675.
Brown J.R., Beckenbach A.T., Smith M.J. 1993. Intraspecific DNA sequence variation of the mitochondrial control region of white sturgeon (Acipenser transmontanus). Mol. Biol. Evol. 10(2): 326–341.
Buroker N.E., Brown J.R., Gilbert T.A., O'Hara P.J., Beckenbach A.T., Thomas W.K., Smith M.J. 1990. Length heteroplasmy of sturgeon mitochondrial DNA: an illegitimate elongation model. Genetics, 124(1): 157–163.
Calabrese F.M., Balacco D.L., Preste R., Diroma M.A., Forino R., Ventura M., Attimonelli M. 2017. NUMTs colonization in mammalian genomes. Sci. Rep. 7(1): 16357.
Clark K.A., Howe D.K., Gafner K., Kusuma D., Ping S., Estes S., et al. 2012. Selfish little circles: transmission bias and evolution of large deletion-bearing mitochondrial DNA in Caenorhabditis briggsae nematodes. PLoS One, 7(7): e41433.
Croteau E.K., Heist E.J., Nielsen C.K., Hutchinson J.R., Hellgren E.C. 2012. Microsatellites and mitochondrial DNA reveal regional population structure in bobcats (Lynx rufus) of North America. Conserv. Genet. 13(6): 1637–1651.
Dell'Era E., Polidori M., Bernardini M., Capomaccio S., Cappelli K., Balducci F., Mandara M.T. 2021. Selective symmetrical necrotizing encephalopathy secondary to primary mitochondrial disorder in a cat. J. Vet. Intern. Med. 35(5): 2401–2408.
Drew A.C., Brindley P.J. 1997. A retrotransposon of the non-long terminal repeat class from the human blood fluke Schistosoma mansoni. Similarities to the chicken-repeat-1-like elements of vertebrates. Mol. Biol. Evol. 14(6): 602–610.
Foran D.R., Crooks K.R., Minta S.C. 1997. Species identification from scat: an unambiguous genetic method. Wildl. Soc. Bull. 25(4): 835–839. Available from https://www.jstor.org/stable/3783732.
Freeman A.R., MacHugh D.E., McKeown S., Waltzer C., McConnell D.J., Bradley D.G. 2001. Sequence variation in the mitochondrial DNA control region of wild African cheetahs (Acinonyx jubatus). Heredity, 86(3): 355–362.
Fumagalli L., Taberlet P., Favre L., Hausser J. 1996. Origin and evolution of homologous repeated sequences in the mitochondrial DNA control region of shrews. Mol. Biol. Evol. 13(1): 31–46.
Graham N.R., Gillespie R.G., Krehenwinkel H. 2021. Towards eradicating the nuisance of numts and noise in molecular biodiversity assessment. Mol. Ecol. Resour. 21(6): 1755–1758.
Gugolz D., Bernasconi M.V., Breitenmoser-Würsten C., Wandeler P. 2008. Historical DNA reveals the phylogenetic position of the extinct alpine lynx. J. Zool. 275(2): 201–208.
Hellborg L., Walker C.W., Rueness E.K., Stacy J.E., Kojola I., Vladmann H., et al. 2002. Differentiation and levels of genetic variation in northern European lynx (Lynx lynx) populations revealed by microsatellites and mitochondrial DNA analysis. Conserv. Genet. 3(2): 97–111.
Hoelzel A.R., Hancock J.M., Dover G.A. 1991. Evolution of the cetacean mitochondrial D-loop region. Mol. Biol. Evol. 8(4): 475–493.
Hoelzel A.R., Lopez J.V., Dover G.A., O'Brien S.J. 1994. Rapid evolution of a heteroplasmic repetitive sequence in the mitochondrial DNA control region of carnivores. J. Mol. Evol. 39(2): 191–199.
Irigoyen P., Pintos-Polasky P., Rosa-Villagran L., Skowronek M.F., Cassina A., Sapiro R. 2022. Mitochondrial metabolismo deteremines the functional status of human sperm and correlates with semen parameters. Front. Cell. Dev. Biol. 10: 926684.
Jackson C.B., Turnbull D.M., Minczuk M., Gammage P.A. 2020. Therapeutic manipulation of mtDNA heteroplasmy: A shifting perspective. Trends Mol. Med. 26(7): 698–709.
Jae-Heup K., Eizirik E., O'Brien S.J., Johnson W.E. 2001. Structure and patterns of sequence variation in the mitochondrial DNA control region of the great cats. Mitochondrion, 1(3): 279–292.
Johnson W.E., Godoy J.A., Palomares F., Delibes M., Fernandes M., Revilla E., O'Brien S.J. 2004. Phylogenetic and phylogeographic analysis of Iberian lynx populations. J. Hered. 95(1): 19–28.
Kearse M., Moir R., Wilson A., Stones-Havas S., Cheung M., Sturrock S., et al. 2012. Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics, 28(12): 1647–1649.
Ketmaier V., Bernardini C. 2005. Structure of the mitochondrial control region of the Eurasian otter (Lutra lutra; Carnivora; Mustelidae): patterns of genetic heterogeneity and implications for conservation of the species in Italy. J. Hered. 96(4): 318–328.
Kopinski P.K., Janssen K.A., Schaefer P.M., Trefely S., Perry C.E., Potluri P., et al. 2019. Regulation of nuclear epigenome by mitochondrial DNA heteroplasmy. Proc.Natl. Acad. Sci. U.S.A. 116(32): 16028–16035.
Krebs C.J., Boonstra R., Boutin S. 2018. Using experimentation to understand the 10-year snowshoe hare cycle in the boreal forest of North America. J. Anim. Ecol. 87(1): 87–100.
Lechuga-Vieco A.V., Latorre-Pellicer A., Calvo E., Torroja C., Pellico J., Acín-Pérez R., et al. 2022. Heteroplasmy of wild-type mitochondrial DNA variants in mice causes metabolic heart disease with pulmonary hypertension and frailty. Circulation, 145(14): 1084–1101.
Li H., Durbin R. 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 25: 1754–1760.
Li F-Y., Cuddon P.A., Song J., Wood S.L., Patterson J.S., Shelton G.D., Duncan I.D. 2006. Canine spongiform leukoencephalomyelopathy is associated with a missense mutation in cytochrome b. Neurobiol. Dis. 21(1): 35–42.
Li G., Davis B.W., Eizirik E., Murphy W.J. 2016. Phylogenomic evidence for ancient hybridization in the genomes of living cats (Felidae). Genome Res. 26: 1–11.
Li H. 2018. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34(18): 3094–3100.
Lopez J.V., Yuhki N., Masuda R., Modi W., O'Brien S.J. 1994. Numt, a recent transfer and tandem amplification of mitochondrial DNA to the nuclear genome of the domestic cat. J. Mol. Evol. 39(1): 174–190.
Lopez J.V., Cevario S., O'Brien S.J. 1996. Complete nucleotide sequences of the domestic cat (Felis catus) mitochondrial genome and a transposed mtDNA tandem repeat (NUMT) in the nuclear genome. Genomics, 33(2): 229–246.
Matthews P.M., Hopkin J., Brown R.M., Stephenson J.B., Hilton-Jones D., Brown G.K. 1994. Comparison of the relative levels of the 3243 (A→G) mtDNA mutation in heteroplasmic adult and fetal tissues. J. Med. Genet. 31(1): 41–44.
Michikawa Y., Mazzucchelli F., Bresolin N., Scarlato G., Attardi G. 1999. Aging-dependent large accumulation of point mutations in the human mtDNA control region for replication. Science, 296(5440): 774–779.
Päckert M., Giacalone G., Lo Valvo M., Kehlmaier C. 2019. Mitochondrial heteroplasmy in an avian hybrid form (Passer italia: Aves, Passeriformes). Mitochondrial DNA B: Resour. 4(2): 3809–3812.
Patterson E.C., Lall G.M., Neumann R., Ottolini B., Batini C., Sacchini F., et al. 2023. Mitogenome sequences of domestic cats demonstrate lineage expansions and dynamic mutation processes in a mitochondrial minisatellite. BMC Genome, 24(1): 690.
Porter T.M., Hajibabaei M. 2021. Profile hidden Markov model sequence analysis can help remove putative pseudogenes from DNA barcoding and metabarcoding datasets. BMC Bioinform. 22(1): 256.
Prentice M.B., Bowman J., Murray D.L., Klütsch C.F.C., Khidas K., Wilson P.J. 2019. Evaluating evolutionary history and adaptive differentiation to identify conservation units of Canada lynx (Lynx canadensis). Glob. Ecol. Conserv. 20: e00708.
Pukazhenthi B.S., Neubauer K., Jewgenow K., Howard J., Wildt D.E. 2006. The impact and potential etiology of teratospermia in the domestic cat and its wild relatives. Theriogenology,66(1): 112–121.
Rensch T., Villar D., Horvath J., Odom D.T., Flicek P. 2016. Mitochondrial heteroplasmy in vertebrates using ChIP-sequencing data. Genome Biol. 17(1): 139.
Rueness E.K., Stenseth N.C., O'Donoghue M., Boutin S., Ellegren H., Jakobsen K.S. 2003. Ecological and genetic spatial structuring in the Canadian lynx. Nature, 425(6953): 69–72.
Schwartz M.K., Pilgrim K.L., McKelvey K.S., Rivera P.T., Ruggiero L.F. 2007. DNA markers for identifying individual snowshoe hares using field-collected pellets. Northwest Sci. 81(4): 316–322.
Sindičić M., Gomerčić T., Galov A., Polanc P., Huber Đ., Slavica A. 2012. Repetitive sequences in Eurasian lynx (Lynx lynx L.) mitochondrial DNA control region. Mitochondrial DNA, 23(3): 201–207.
Terrell K.A., Wildt D.E., Anthony N.M., Bavister B.D., Leibo S.P., Penfold L.M., et al. 2010. Evidence for compromised metabolic function and limited glucose uptake in spermatozoa from teratospermic domestic cat (Felis catus) and cheetah (Acinonyx jubatus). Biol. Reprod. 83: 833–841.
Wilkinson G.S., Chapman A.M. 1991. Length and sequence variation in evening bat D-loop mtDNA. Genetics, 128(3): 607–617.
Wilkinson G.S., Mayer F., Kerth G., Petri B. 1997. Evolution of repeated sequence arrays in the D-loop region of bat mitochondrial DNA. Genetics, 146(3): 1035–1048.
Wilton S.D., Lim L., Dye D., Laing N. 1997. Bandstab: a PCR-based alternative to cloning PCR products. BioTechniques, 22(4): 642–645.
Wu X., Zheng T., Jiang Z., Wei L. 2007. The mitochondrial genome structure of the clouded leopard (Neofelis nebulosa). Genome, 50(2): 252–257.
Information & Authors
Information
Published In
Genome
Volume 67 • Number 12 • December 2024
Pages: 493 - 502
History
Received: 11 September 2023
Accepted: 8 August 2024
Accepted manuscript online: 3 September 2024
Version of record online: 12 November 2024
Copyright
© 2024 The Author(s). This work is licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.
Data Availability Statement
Sanger sequences were uploaded to GenBank (left domain OR338890-OR338894, right domain OR338895-OR338915).
Key Words
Authors
Author Contributions
Conceptualization: KJF, MAB, KN, JCG
Data curation: KJF, MAB, KN, JKJ
Formal analysis: KJF, MAB, KN, EWH, JKJ, JCG
Funding acquisition: JKJ, WLH, JCG
Investigation: KJF, MAB, KN, EWH
Methodology: KJF, MAB, KN, EWH
Project administration: WLH, JCG
Resources: EWH, JKJ, WLH, JCG
Software: KJF, MAB, KN, JKJ, JCG
Supervision: JCG
Validation: KJF, MAB, KN, JKJ, JCG
Visualization: KJF, MAB, KN, JCG
Writing – original draft: KJF, MAB, JCG
Writing – review & editing: KJF, MAB, KN, EWH, JKJ, WLH, JCG
Competing Interests
The authors declare there are no competing interests.
Funding Information
Natural Sciences and Engineering Research Council (NSERC): RGPIN-2018-06764, RGPIN-2020-04475
Northwest Boreal Landscape Conservation Cooperative
This project was funded by Discovery Grants from the Natural Sciences and Engineering Research Council (NSERC) of Canada to JKJ (RGPIN-2020-04475) and JCG (RGPIN-2018-06764) and a grant from the Northwest Boreal Landscape Conservation Cooperative to WLH. KJF and MAB were supported by NSERC undergraduate student research awards, while EWH was supported by a Mitacs Accelerate Fellowship.
Metrics & Citations
Metrics
Other Metrics
Citations
Cite As
Krystyn J. Forbes, McIntyre A. Barrera, Karsten Nielsen-Roine, Evan W. Hersh, Jasmine K. Janes, William L. Harrower, and Jamieson C. Gorrell. 2024. Stabilizing selection and mitochondrial heteroplasmy in the Canada lynx (Lynx canadensis). Genome.
67(12): 493-502. https://doi.org/10.1139/gen-2023-0094
Export Citations
If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.
There are no citations for this item
View Options
View options
Login options
Check if you access through your login credentials or your institution to get full access on this article.