scispace - formally typeset

Posted ContentDOI

Mitonuclear co-introgression opposes genetic differentiation between phenotypically divergent songbirds

08 Aug 2021-bioRxiv (Cold Spring Harbor Laboratory)-

AbstractComparisons 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.

Topics: Introgression (55%), Mitochondrial DNA (51%), Nuclear gene (50%)

Summary (4 min read)

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.

Did you find this useful? Give us your feedback

...read more

Content maybe subject to copyright    Report

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

Figures (10)
References
More filters

Journal Article
TL;DR: Copyright (©) 1999–2012 R Foundation for Statistical Computing; permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and permission notice are preserved on all copies.
Abstract: Copyright (©) 1999–2012 R Foundation for Statistical Computing. Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies. Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying, provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one. Permission is granted to copy and distribute translations of this manual into another language, under the above conditions for modified versions, except that this permission notice may be stated in a translation approved by the R Core Team.

229,202 citations


Journal ArticleDOI
Abstract: Summary: The Sequence Alignment/Map (SAM) format is a generic alignment format for storing read alignments against reference sequences, supporting short and long reads (up to 128 Mbp) produced by different sequencing platforms. It is flexible in style, compact in size, efficient in random access and is the format in which alignments from the 1000 Genomes Project are released. SAMtools implements various utilities for post-processing alignments in the SAM format, such as indexing, variant caller and alignment viewer, and thus provides universal tools for processing read alignments. Availability: http://samtools.sourceforge.net Contact: [email protected]

35,747 citations


Journal ArticleDOI
TL;DR: Burrows-Wheeler Alignment tool (BWA) is implemented, a new read alignment package that is based on backward search with Burrows–Wheeler Transform (BWT), to efficiently align short sequencing reads against a large reference sequence such as the human genome, allowing mismatches and gaps.
Abstract: Motivation: The enormous amount of short reads generated by the new DNA sequencing technologies call for the development of fast and accurate read alignment programs. A first generation of hash table-based methods has been developed, including MAQ, which is accurate, feature rich and fast enough to align short reads from a single individual. However, MAQ does not support gapped alignment for single-end reads, which makes it unsuitable for alignment of longer reads where indels may occur frequently. The speed of MAQ is also a concern when the alignment is scaled up to the resequencing of hundreds of individuals. Results: We implemented Burrows-Wheeler Alignment tool (BWA), a new read alignment package that is based on backward search with Burrows–Wheeler Transform (BWT), to efficiently align short sequencing reads against a large reference sequence such as the human genome, allowing mismatches and gaps. BWA supports both base space reads, e.g. from Illumina sequencing machines, and color space reads from AB SOLiD machines. Evaluations on both simulated and real data suggest that BWA is ~10–20× faster than MAQ, while achieving similar accuracy. In addition, BWA outputs alignment in the new standard SAM (Sequence Alignment/Map) format. Variant calling and other downstream analyses after the alignment can be achieved with the open source SAMtools software package. Availability: http://maq.sourceforge.net Contact: [email protected]

35,234 citations


Journal ArticleDOI
TL;DR: Timmomatic is developed as a more flexible and efficient preprocessing tool, which could correctly handle paired-end data and is shown to produce output that is at least competitive with, and in many cases superior to, that produced by other tools, in all scenarios tested.
Abstract: Motivation: Although many next-generation sequencing (NGS) read preprocessing tools already existed, we could not find any tool or combination of tools that met our requirements in terms of flexibility, correct handling of paired-end data and high performance. We have developed Trimmomatic as a more flexible and efficient preprocessing tool, which could correctly handle paired-end data. Results: The value of NGS read preprocessing is demonstrated for both reference-based and reference-free tasks. Trimmomatic is shown to produce output that is at least competitive with, and in many cases superior to, that produced by other tools, in all scenarios tested. Availability and implementation: Trimmomatic is licensed under GPL V3. It is cross-platform (Java 1.5+ required) and available at http://www.usadellab.org/cms/index.php?page=trimmomatic Contact: ed.nehcaa-htwr.1oib@ledasu Supplementary information: Supplementary data are available at Bioinformatics online.

26,464 citations


Journal ArticleDOI
TL;DR: The purpose of this discussion is to offer some unity to various estimation formulae and to point out that correlations of genes in structured populations, with which F-statistics are concerned, are expressed very conveniently with a set of parameters treated by Cockerham (1 969, 1973).
Abstract: This journal frequently contains papers that report values of F-statistics estimated from genetic data collected from several populations. These parameters, FST, FIT, and FIS, were introduced by Wright (1951), and offer a convenient means of summarizing population structure. While there is some disagreement about the interpretation of the quantities, there is considerably more disagreement on the method of evaluating them. Different authors make different assumptions about sample sizes or numbers of populations and handle the difficulties of multiple alleles and unequal sample sizes in different ways. Wright himself, for example, did not consider the effects of finite sample size. The purpose of this discussion is to offer some unity to various estimation formulae and to point out that correlations of genes in structured populations, with which F-statistics are concerned, are expressed very conveniently with a set of parameters treated by Cockerham (1 969, 1973). We start with the parameters and construct appropriate estimators for them, rather than beginning the discussion with various data functions. The extension of Cockerham's work to multiple alleles and loci will be made explicit, and the use of jackknife procedures for estimating variances will be advocated. All of this may be regarded as an extension of a recent treatment of estimating the coancestry coefficient to serve as a mea-

16,821 citations