Mitonuclear co-introgression opposes genetic differentiation between phenotypically divergent songbirds

TL;DR: In this article, the mtDNA introgression process can select for co-introgression at nuclear genes with mitochondrial functions (mitonuclear genes), which may contribute to continued hybridization between yellowhammers and pine buntings despite their clear morphological and genetic differences.

Abstract: Comparisons of genomic variation among closely related species often show more differentiation in mitochondrial DNA (mtDNA) and sex chromosomes than in autosomes, a pattern expected due to the relative effective population sizes of these genomic components. Differential introgression can cause some species pairs to deviate dramatically from this pattern. The yellowhammer (Emberiza citrinella) and the pine bunting (E. leucocephalos) are hybridizing avian sister species that differ greatly in appearance but show no mtDNA differentiation. This discordance might be explained by mtDNA introgression--a process that can select for co-introgression at nuclear genes with mitochondrial functions (mitonuclear genes). We investigated genome-wide nuclear differentiation between yellowhammers and pine buntings and compared it to what was seen previously in the mitochondrial genome. We found clear nuclear differentiation that was highly heterogeneous across the genome, with a particularly wide differentiation peak on the sex chromosome Z. We further tested for preferential introgression of mitonuclear genes and detected evidence for such biased introgression in yellowhammers. Mitonuclear co-introgression can remove post-zygotic incompatibilities between species and may contribute to the continued hybridization between yellowhammers and pine buntings despite their clear morphological and genetic differences. As such, our results highlight the potential ramifications of co-introgression in species evolution.

Summary

Introduction

• Evolution in eukaryotes is shaped by changes in multiple genomic components that differ in their modes of inheritance: mitochondrial DNA is usually inherited through the matrilineal line, autosomes are inherited through both parental lines and sex chromosomes are inherited differentially depending on the sex of both parent and offspring (Avise, 2000) .
• There is often much variation among these genomic components in the degree of genetic differentiation between related populations or species (reviewed in Coyne & Orr, 2004; reviewed in Price, 2008) , suggesting that their dynamics differ during the process of speciation of a single species into two or more.
• Differentiation across autosomes, which tends to be lower than on mtDNA and sex chromosomes, can be highly heterogeneous.
• Research investigating coevolution between mitochondrial and nuclear genomes is relatively novel as mtDNA was often treated as a neutral marker in past evolutionary research (Avise, 2000) .
• Genomic work has identified mitonuclear discordance between allopatric yellowhammers and pine buntings (Irwin et al. 2009) as they possess almost no mtDNA divergence but show moderate differentiation in nuclear AFLP (Amplified Fragment Length Polymorphism) markers.

DNA extraction and genotyping-by-sequencing

• DNA was extracted from samples using a standard phenol-chloroform method.
• The authors then divided the DNA samples into four genotyping-by-sequencing (GBS) libraries (Elshire et al. 2011) .
• The 109 samples included in this study were sequenced in libraries together with 226 yellowhammer, pine bunting and hybrid DNA samples collected near and within the sympatric zone as part of a larger project (Nikelski et al. in prep) .
• The libraries were prepared as per the protocol described by Alcaide et al. (2014) with the modifications specified by Geraldes et al. (2019) except that the authors maintained a 300-400 bp fragment size during size selection.
• Paired-end sequencing was completed by Genome Québec using the Illumina HiSeq 4000 system, producing more than 1.2 billion reads, 150 bp in length, across the four GBS libraries.

Genotyping-by-sequencing data filtering

• Processing of GBS reads for the samples analyzed in this study was done in conjunction with reads from the samples included in the larger project mentioned above.
• The resulting VCF files were filtered using VCFtools and GATK to remove indels, sites with more than two alleles, sites with more than 60% missing genomic data, sites with MQ values lower than 20 and sites with heterozygosities greater than 60% (to avoid potential paralogs).
• Use of these filters simplified calculations in downstream analyses and ensured that these analyses were restricted to sites with sufficient data.
• All following statistical analyses were completed using R version 3.6.2 (R Core Team, 2014).

Variant site analyses

• The genome-wide "variant site" VCF file was analyzed using modified versions of the R scripts described in Irwin et al. (2018) .
• A total of 374,780 SNPs were identified among allopatric yellowhammers and pine buntings.
• For each of these SNPs the authors calculated sample size, allele frequency, and Weir and Cockerham's FST (Weir & Cockerham, 1984) .
• The PC1 loadings were also graphed as a Manhattan plot using the package qqman (Turner, 2018) .
• The remaining SNPs did not possess known genomic locations and, therefore, could not be included in the plot.

Differentiation across the genome

• To thoroughly investigate genomic differentiation between allopatric yellowhammers and pine buntings, the authors performed further analysis on both variant and invariant loci within "info site" VCF files using R scripts described in Irwin et al. (2018) .
• The authors calculated Weir and Cockerham's FST, between-group nucleotide distance (! ! ) and within-group nucleotide diversity (! " ) for nonoverlapping windows of available sequence data across each chromosome.
• Significantly positive Tajima's D suggests that there are fewer rare alleles in a population than expected under neutrality, potentially stemming from balancing selection or rapid population contraction.
• Chromosome Z in particular showed a large peak in FST with several SNPs possessing values close to one.
• Additionally, the authors detected a weak negative correlation between the windowed averages of FST and !.

Phylogenetic comparison with other Emberizidae species

• Between allopatric yellowhammers and allopatric pine buntings as well as among these focal species and six other Emberizidae species (Emberiza aureola, Emberiza calandra, Emberiza cioides, Emberiza cirlus, Emberiza hortulana and Emberiza stewarti) to estimate a phylogeny.
• Values for each species pair was converted into a distance matrix and used to create an unrooted neighbour-joining tree.
• This tree was constructed using the ape package (Paradis & Schliep, 2019) and the BioNJ algorithm (Gascuel, 1997) with Emberiza aureola set as the outgroup (Alström et al. 2008) .
• Values between yellowhammers, pine buntings and six other Emberizidae species depicted similar species relationships as were identified previously using mitochondrial markers (Alström et al.
• In terms of the relative genetic distance (! ! ) between yellowhammers and pine buntings compared to the distance between each of those and E. stewarti, nuclear genetic distance was 11.4 times greater than mitochondrial genetic distance.

Signals of mitonuclear co-introgression

• To test for signals of mitonuclear gene introgression between allopatric yellowhammers and pine buntings, the authors compiled a list of mitonuclear genes and a list of 2000 bp putative introgression windows (hereafter referred to as "introgression windows") and then tested for an association between them.
• It should be noted that sharing of some introgression windows is expected given that the contribution of ! ! to window selection was identical for both taxa (in contrast, Tajima's D was calculated separately for yellowhammers and pine buntings).
• This finding was significant in a two-tailed binomial test (p = 0.04952) indicating that mitonuclear genes appeared in yellowhammer introgression windows more often than would be expected if they were assigned to windows randomly.
• Four mitonuclear genes appeared within pine bunting introgression windows-3.0% of the genes considered.
• All four genes appeared on separate autosomes with three of these genes encoding structural subunits of the ETC and one encoding a protein subunit of the mitochondrial ribosome.

Results

• When comparing allopatric yellowhammers and pine buntings, following filtering the authors identified 374,780 variable SNPs within their "variant site" VCF file and 13,703,455 invariant and 699,122 variant sites across thirty autosomes and the Z chromosome within their "info site" VCF files.
• In the latter "info site" files, the authors designated a total of 7187 genomic windows (of 2000 sequenced bp each) across the genome, with each window covering an average distance of about 139 kilobases.

Overall genetic differentiation

• Based on 374,780 SNPs considered all together, their genome-wide FST estimate was 0.0232 between allopatric yellowhammers and pine buntings.
• PC1 explained 3.6% of the variation among individuals while PC2 explained 2.9% of the variation.
• Further investigation into these outliers revealed that they were males from the same location.
• To explore the causes of these outliers, the authors temporarily removed one of them and re-ran the PCA.
• It is unclear what is responsible for these outliers, but the distinct yellowhammer and pine bunting genetic clusters remained intact in all the PCAs considered.

Discussion

• Yellowhammers and pine buntings show negligible mtDNA differentiation (Irwin et al. 2009 ) but are well differentiated phenotypically (Panov et al.
• This scenario is consistent with the observed extensive hybridization between these taxa (Panov et al.
• The large regions of the Z chromosome that have FST values near zero suggests that additional factors are involved in producing the large and wide island of differentiation on the Z.

in prep).

• While numerous "islands of differentiation" were observed between yellowhammers and pine buntings implying moderate genetic divergence between them, mtDNA introgression has the potential to homogenize their nuclear genomes at mitonuclear genes by selecting for cointrogression of compatible alleles (Beck et al.
• Because a comparable signal of introgression was not found in allopatric pine buntings, the authors suggest that mitonuclear co-introgression could have occurred in the direction of pine buntings into yellowhammers.
• Unlike the protein-protein interactions occurring within ETC complexes, mitonuclear interactions in the mitoribosome are between nuclear-encoded proteins and mitochondrialencoded RNA (Hill,. 2019) .
• As well, the inclusion of analyses that compare mtDNA and mitonuclear gene differentiation in a wider range of systems would help to clarify the potentially important role that mitonuclear interactions play in the merging or diverging of species.
• Each photo represents one of four phenotypic classes: PC, SC, PL and SL.

Figures

Figure 2. Unrooted neighbour-joining tree of Emberizidae species constructed based on average absolute between-population nucleotide diversity (!!). Sample sizes for each species are as follows: E. aureola = 1, E. calandra = 1, E. cioides = 1, E. hortulana = 1, E. cirlus = 6, E. stewarti = 4, E. citrinella = 53 and E. leucocephalos = 42.

Table 2. Summary statistics calculated while conducting mitonuclear co-introgression analysis. A total of 7187 windows, each of 2000 bp of obtained sequence, were considered when determining introgression windows. A total of 134 mitonuclear genes were investigated for signals of co-introgression. “*” indicates a significant p-value.

Figure 6. Mean absolute within-group nucleotide diversity (!") of allopatric yellowhammers (n = 53) and allopatric pine buntings (n =42) plotted against absolute between-group nucleotide diversity (!!). Each dot represents the average value taken from a 2000 bp window of sequenced data across the nuclear genome. The black line indicates where mean within-group nucleotide diversity equals between-group nucleotide diversity. Increasing values of relative differentiation (FST) calculated for each window are shown in darker shades of blue.

Table 3. Identities, chromosomal locations, windowed Tajima’s D values and functions of mitonuclear genes that appeared within 244 yellowhammer introgression windows. In the “Mitonuclear Gene Function” column, ETC stands for “Electron Transport Chain”. Mitonuclear gene names are written as they appear in Hill (2019).

Figure 5. Patterns of genetic variation comparing allopatric yellowhammers (n = 53) and allopatric pine buntings (n = 42) across three chromosomes (2, 5 and Z). Relative nucleotide differentiation (FST), absolute between-population nucleotide diversity (!!), absolute withinpopulation nucleotide diversity (!") and Tajima’s D (TajD) are shown as 2000 bp windowed averages across each chromosome. FST and !! are shown as purple lines to indicate that values were calculated as a comparison between allopatric yellowhammers and pine buntings. !" and TajD are shown as two separate lines (yellow = yellowhammers, brown = pine buntings) to indicate that values were calculated separately for each population.

Figure 4. Relative differentiation (FST) of 349,807 genome-wide SNPs identified between allopatric yellowhammers (n = 53) and allopatric pine buntings (n = 42), with chromosomes represented with alternating black and grey. Narrow regions of elevated differentiation can be seen on many autosomes, and there are broad regions of high differentiation on the Z chromosome.

Table 4. Identities, chromosomal locations, windowed Tajima’s D values and functions of mitonuclear genes that appeared within 222 pine bunting introgression windows. In the “Mitonuclear Gene Function” column, ETC stands for “Electron Transport Chain”. Mitonuclear gene names are written as they appear in Hill (2019).

Figure 3. PCA of genetic variation between allopatric yellowhammers (yellow; n = 53) and allopatric pine buntings (brown; n = 42), based on 374,780 genome-wide SNPs. PC1 and PC2 explain 3.6% and 2.9%, respectively, of the variation among individuals.

Table 1. Geographical locations and sample sizes of the sites included in this study. Sampling locations may include multiple sites that appeared too close together to be shown in detail in Figure 1A. Full details for the sites included in each sampling locations can be found in Supplementary Table 1. The sampling location numbers that appear in the “Sampling Location” column correspond to those that appear in red in Figure 1A. The “Sample Size” columns describes the total number of samples collected from a particular location.

Figure 7. Association between relative differentiation (FST) and absolute between-group nucleotide diversity (!!) of allopatric yellowhammers (n=53) and allopatric pine buntings (n = 42). Each black dot represents average values calculated from a 2000 bp window of sequenced data. A cubic spline fit between the variables is shown as a purple line.

Full text

1
Mitonuclear co-introgression opposes genetic differentiation between
1
phenotypically divergent songbirds
2
3
Ellen Nikelski
1, *
, Alexander S. Rubtsov
2
, and Darren Irwin
1
4
5
1
Department of Zoology, and Biodiversity Research Centre, 6270 University Blvd., University of British Columbia,
6
Vancouver, BC, Canada
7
2
State Darwin Museum, Moscow, Russia
8
*
Present address: Department of Ecology and Evolutionary Biology, University of Toronto, Toronto, ON, Canada
9
10
Running Title:
11
Co-introgression opposes differentiation
12
13
Corresponding Author:
14
Ellen Nikelski, Department of Ecology and Evolutionary Biology, University of Toronto,
15
Toronto, ON, Canada. Email: ellen.nikelski@mail.utoronto.ca
16
17
Keywords:
18
mitonuclear co-introgression, mtDNA, mitonuclear gene, genetic differentiation, chromosome Z,
19
Aves
20
21
22
23
24
25
26
27
28
29
30
31
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 8, 2021. ; https://doi.org/10.1101/2021.08.08.455564doi: bioRxiv preprint

2
Abstract
32
Comparisons of genomic variation among closely related species often show more differentiation
33
in mitochondrial DNA (mtDNA) and sex chromosomes than in autosomes, a pattern expected
34
due to the relative effective population sizes of these genomic components. Differential
35
introgression can cause some species pairs to deviate dramatically from this pattern. The
36
yellowhammer (Emberiza citrinella) and the pine bunting (E. leucocephalos) are hybridizing
37
avian sister species that differ greatly in appearance but show no mtDNA differentiation. This
38
discordance might be explained by mtDNA introgression—a process that can select for co-
39
introgression at nuclear genes with mitochondrial functions (mitonuclear genes). We investigated
40
genome-wide nuclear differentiation between yellowhammers and pine buntings and compared it
41
to what was seen previously in the mitochondrial genome. We found clear nuclear differentiation
42
that was highly heterogeneous across the genome, with a particularly wide differentiation peak
43
on the sex chromosome Z. We further tested for preferential introgression of mitonuclear genes
44
and detected evidence for such biased introgression in yellowhammers. Mitonuclear co-
45
introgression can remove post-zygotic incompatibilities between species and may contribute to
46
the continued hybridization between yellowhammers and pine buntings despite their clear
47
morphological and genetic differences. As such, our results highlight the potential ramifications
48
of co-introgression in species evolution.
49
50
Introduction
51
Evolution in eukaryotes is shaped by changes in multiple genomic components that differ
52
in their modes of inheritance: mitochondrial DNA (mtDNA) is usually inherited through the
53
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 8, 2021. ; https://doi.org/10.1101/2021.08.08.455564doi: bioRxiv preprint

3
matrilineal line, autosomes are inherited through both parental lines and sex chromosomes are
54
inherited differentially depending on the sex of both parent and offspring (Avise, 2000). There is
55
often much variation among these genomic components in the degree of genetic differentiation
56
between related populations or species (reviewed in Coyne & Orr, 2004; reviewed in Price,
57
2008), suggesting that their dynamics differ during the process of speciation of a single species
58
into two or more. This variation can arise through differences in both the rate at which specific
59
DNA sequences evolve and the degree to which different components contribute towards genetic
60
incompatibilities that reduce gene flow between populations. A common pattern observed
61
between speciating taxa is clear differentiation in mtDNA (eg. Hebert et al. 2004; Kerr et al.
62
2007), moderate differentiation in sex chromosomes (eg. Thornton & Long, 2002; Borge et al.
63
2005; Lu & Wu, 2005; Harr, 2006; Ruegg et al. 2014; Sackton et al. 2014), and comparatively
64
modest differentiation across autosomes (Harr, 2006; Nadeau et al. 2012; Irwin et al. 2018).
65
Measures of mtDNA differentiation are often used to identify and classify genetically
66
distinct populations (eg. Hebert et al. 2004; Kerr et al. 2007) and to infer their histories (Moore,
67
1995; Zink & Barrowclough, 2008). Due to its uniparental inheritance, mtDNA has one quarter
68
the effective population size and coalescence time of autosomal nuclear DNA (Moore, 1995).
69
This characteristic combined with mtDNA’s relatively high mutation rate (Lynch et al. 2006)
70
mean that genetic differences arise and fix relatively quickly, creating patterns of clear mtDNA
71
differentiation between recently diverged populations.
72
Sex chromosomes are another genomic region that often shows higher between-
73
population genetic differentiation compared to autosomes between speciating taxa, in both Z/W
74
(Borge et al. 2005; Ruegg et al. 2014; Sackton et al. 2014) and X/Y systems (Thorton & Long,
75
2002; Lu & Wu, 2005; Harr, 2006). To explain this “faster Z/X effect,” researchers have noted
76
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 8, 2021. ; https://doi.org/10.1101/2021.08.08.455564doi: bioRxiv preprint

4
that, because beneficial recessive mutations on the Z or X chromosome are immediately exposed
77
to selective forces in the heterogametic sex, fixation of these mutations should proceed faster
78
than if the mutations appeared on autosomes (reviewed in Meisel & Connallon, 2013; Irwin,
79
2018). Also contributing to genetic differentiation on the Z and X chromosomes are the lower
80
effective population sizes of these chromosomes compared to autosomes (Mank et al. 2010;
81
reviewed in Irwin, 2018). A lower effective population size allows for the fixation of a greater
82
number of slightly deleterious mutations due to less effective purifying selection and a larger role
83
of genetic drift. It is likely that both forces—the faster Z/X effect and less effective purifying
84
selection—contribute to the moderate amount of genetic differentiation seen between the sex
85
chromosomes of diverging taxa (Thorton & Long, 2002; Borge et al. 2005; Lu & Wu, 2005;
86
Harr, 2006; Ruegg et al. 2014; Sackton et al. 2014).
87
Differentiation across autosomes, which tends to be lower than on mtDNA and sex
88
chromosomes, can be highly heterogeneous. In fact, many researchers report “islands of
89
differentiation” on autosomes where peaks of high relative differentiation are found against a
90
background of low relative differentiation (e.g., Harr, 2006; Nadeau et al. 2012; Hejase et al.
91
2020). Explanations for these “islands” usually invoke reduced gene flow (reviewed in Wu,
92
2001) and/or repeated bouts of selection (Cruickshank and Hahn, 2014; Irwin et al. 2018). In the
93
former scenario, differentiation peaks are hypothesized to house the loci responsible for
94
reproductive barriers between interacting taxa and, as a result, they are resistant to the gene flow
95
that homogenizes the rest of the nuclear genome. In contrast, explanations invoking repeated
96
selection hypothesize that differentiation islands are areas of the genome that experienced
97
repeated reductions in genetic diversity as a result of selection or selective sweeps in both
98
ancestral and daughter populations.
99
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 8, 2021. ; https://doi.org/10.1101/2021.08.08.455564doi: bioRxiv preprint

5
Despite the general patterns of differentiation discussed above, an increasing number of
100
studies report remarkably low differentiation between populations at what are normally highly
101
divergent genetic components when compared to other genetic regions or observable phenotypes
102
(e.g., Irwin et al. 2009; Yannic et al. 2010; Bryson et al. 2012). In a number of cases, mtDNA
103
shows dramatically low differentiation when compared to differentiation of the nuclear genome,
104
a pattern referred to as “mitonuclear discordance” (reviewed in Toews & Brelsford, 2012).
105
Discordance between marker types may be explained by hybridization and introgression between
106
populations, perhaps due to a selective advantage of the introgressing genetic region. For
107
example, Hulsey et al. (2016) documented low mtDNA differentiation—likely due to
108
introgression—and clear differentiation in nuclear DNA (nucDNA) between two hybridizing
109
cichlid species (Hulsey & García de León, 2013). The researchers further reported high mtDNA
110
differentiation between isolated populations of cichlids at genetic sites associated with thermal
111
tolerance and a significant correlation between mtDNA divergence and water temperature
112
(Hulsey et al. 2016). Altogether, these results suggest that mtDNA introgression produced the
113
discordance seen between marker types and that this outcome was potentially driven by adaptive
114
selection for tolerance of extreme water temperatures.
115
The hypothesis of adaptive introgression increases in complexity if we consider the
116
potential for coevolution between genomic components. Research investigating coevolution
117
between mitochondrial and nuclear genomes is relatively novel as mtDNA was often treated as a
118
neutral marker in past evolutionary research (Avise, 2000). Nevertheless, recent empirical and
119
theoretical work has provided greater context regarding how mitonuclear coevolution may
120
influence the progression of differentiation and speciation between taxa (Hill, 2019).
121
.CC-BY-NC-ND 4.0 International licenseavailable under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprintthis version posted August 8, 2021. ; https://doi.org/10.1101/2021.08.08.455564doi: bioRxiv preprint