Environmental pollution is a widespread problem that living organisms have to contend with on a global scale. In contaminated sites especially, wild populations undergo intense selective pressure that may result in phenotypic adaptations to pollutants (Hendry et al., 2008). The scientific article (Reid et al., 2016) discussed in this blogpost explores the genetic mechanisms that have allowed the rapid adaptation to industrial pollutants in wild Atlantic killifish populations.
Results
The genomic landscape of the killifish populations
Atlantic killifish (Fundulus heteroclitus) are non-migratory fish that are abundant along the US east coastline (Fig. 1A). Some killifish populations show inherited resistance to lethal levels of industrial pollutants in sites that have been contaminated for decades. For instance, the authors show that the percentage of larva that survive in increasing concentrations of a highly toxic pollutant called PCB 126, is higher in tolerant populations compared to the sensitive populations (Fig. 1B). To understand the genetic adaptations underlying the rapid adaptation to polluted sites in killifish populations, the authors sequenced the complete genomes from eight populations. Four tolerant populations that reside in highly polluted sites were sampled. Each one was paired with a sensitive population from a nearby site (Fig. 1A). The authors combined these genomic data with corresponding RNA sequencing (RNA-seq) to identify unique and shared pathways among tolerant populations as well as to uncover adaptive evidence in the populations.
The genomes from 43-50 individuals from each population were sequenced. One pair of tolerant and sensitive populations (T1 and S1) were sequenced to 7-fold coverage, while the remaining populations, to 0.6-fold coverage. These data indicate that the populations’ genetic variation is strongly by their geographical locations. Meanwhile, all tolerant and sensitive pairs of populations share the most similar genomic backgrounds and have low Fst values between them (0.01 – 0.08). Additionally, tolerant populations show a lower genome-wide nucleotide diversity (?) along with a positive-shifted Tajima’s D. Thus, the authors conclude that tolerant populations have recently and independently diverged from local ancestral populations.
Signatures of convergent evolution in tolerant killifish populations
To identify genomic regions responsible for conferring pollution tolerance in killifish, the authors scanned the populations’ genomes looking for signals of selective sweeps in 5 kb sliding windows. Candidate regions were defined as those showing low values of genetic diversity and a skewed allele frequency spectrum (0.1% tails of ? and Tajima’s D, respectively) as well as high allele frequency differentiation (99.9% tails for Fst). Each tolerant population showed prevalent selection signatures compared to their sensitive counterparts (as seen by ? and Fst). Most of these outlier regions are small (52 – 69 kb, up to ~1.8 Mb) and specific to each tolerant population. Nevertheless, the highest ranked outlier regions are shared between tolerant populations (Fig. 2A). The shared outlier regions harbour genes involved in the aryl hydrocarbon receptor (AHR) signaling pathway (AHR2a, AHR1a, AIP, and CYP1A) (Fig. 2B). These results suggest repeated convergent evolution of pollutant tolerance in the sampled killifish populations.
The authors then tested whether the genes located in outlier regions showed distinct expression profiles in tolerant killifish. Individuals from sensitive and tolerant populations were reared in a common, clean environment for two generations. Following this, embryos were challenged with the toxic pollutant PCB 126 and RNA was collected ~10 days post fertilization. Indeed, AHR-regulated genes were less induced in individuals from tolerant populations (Fig. 2C). Concomitantly, AHR-regulated genes were enriched (P < 0.0001) in the set of genes that were up-regulated in response to PCB 126 treatment in sensitive populations exclusively. Notably, some of the dominant pollutants at the sampled “T” sites bind AHR. Also, aberrant AHR signalling leads to embryo and larval lethality (Pohjanvirta, 2011). The authors thus conclude that the AHR signalling pathway is a key and repeated target of natural selection in polluted sites given the multiple, independent “desensitizing” events in tolerant killifish populations.
It is important to note that genome sequencing coverage seems to have an effect on the ranking of outlier regions. For instance, the regions that contain key AHR-signalling genes (AIP, CYP1A, AHR1a/2a, and ARNT) are very highly ranked in low-coverage populations whereas they are lowly ranked in the high-coverage population pair (T1-S1). Given that outlier regions are ranked based on Fst and nucleotide diversity, these measures must be impacted by low genome sequencing coverage. It would be interesting to determine the ranking of the outlier regions if the other populations were sequenced to higher coverage. However, despite being lowly ranked, these regions are classified as outliers in all four population pairs, giving strength to the argument that impaired AHR-signalling is key to pollution tolerance.
In-depth analysis of genetic variants in tolerant populations
There is evidence for selection of AHR pathway genes in tolerant killifish populations. Tolerant populations harbour distinct deletions spanning AHR2a and AHR1a. On the contrary, individuals from their sensitive counterparts are almost completely devoid of such mutations. Furthermore, RNA-seq data revealed the expression of a chimeric transcript, part AHR2a, part AHR1a in T4 individuals. Meanwhile, AIP (a regulator of AHR stability and cellular localization) is found within a region showing the strongest signals of selection that is shared between all tolerant populations. CYP1A (a transcriptional target of AHR) is also in located in top-ranking outlier regions in all tolerant populations (except for T1 where the region is ranked #401). Interestingly, CYP1A duplications are found in high frequencies in tolerant populations, without a concomitant increase in expression. The authors hypothesize that CYP1A duplications may function as a dosage-compensation mechanism in tolerant populations with impaired AHR signalling because it has been reported that AHR knockout decreases CYP1A expression in rodents (Schmidt et al., 1996). Finally, other AHR-related genes lie within population-specific outlier regions such as the tandem paralogs AHR1b and AHR2b in T3 and T4 and five other AHR pathway genes in T4. Together these observations led the authors to conclude that AHR pathway genes are indeed common and repeated targets of selection, a clear example of convergent evolution.
Genes outside of the AHR signalling pathway are also targets of selection. For example, two genes that are implicated in AHR-independent cardiotoxicity (KCNB2 and KCNC3) are within outlier regions in T4, where such cardiotoxic pollutants are abundant. Additionally, the authors found adaptations that may compensate for the potential costs of pollutant tolerance. AHR signalling is interconnected with multiple other pathways, such as estrogen and hypoxia signalling, as well as cell cycle and immune system regulation (Beischlag et al., 2008). Consequently, estrogen receptor 2b lies within an outlier region in T2, while estrogen receptor-regulated genes are enriched in the gene set of the outlier regions in all tolerant populations (P < 0.001). Furthermore, the estrogen receptor is inferred as an upstream regulator of differentially expressed genes between tolerant and sensitive killifish (Fig. 2C). Alternatively, the hypoxia-inducible factor 2? is in an outlier window in T3, and interleukin and cytokine receptors are in outlier regions in T4. Thus, the authors highlighted the possibility that compensatory adaptation selection may be common following rapid adaptive evolution.
Conclusions
In this article, we can appreciate that genetic adaptations to pollution in wild killifish populations are complex. The authors attribute this to two main factors. Firstly, sites are contaminated with a complex mix of pollutants. This may affect how the AHR- and other pathways are impacted. Therefore, adaptations in multiple pathways, at different genetic levels, may be necessary for tolerance to diverse pollutant mixes to arise. Second, AHR pathway genes are interconnected with other gene-regulatory pathways, thus these genes’ functions may be impaired upon aberrant AHR signalling, It follows that adaptations that compensate these genes’ functions may also be selected for in tolerant populations.
The authors argue that their data clearly reveal signals of convergent evolution. The AHR pathway genes are shown to be repeated targets of selection in distinct pollutant-tolerant killifish populations. This also suggests molecular constraints in the adaptation to pollution. However, in spite of this, multiple variants were favoured in different tolerant populations. The authors say that their data show evidence of selection of preexisting common variants in multiple tolerant populations. In other words, it seems that soft sweeps have been important for the emergence of pollutant tolerance in killifish. This conclusion is supported by several lines of evidence: 1) the sensitive-tolerant populations are genetically close, which suggests that the selected variants were part of the standing variation, 2) sensitive populations have some of the variants that tolerant populations have, and finally 3) these fish are low dispersal. Interestingly, the authors point out that Atlantic killifish have large population size and a wide range of standing genetic variation. These little critters were not only the first space-going-fish, they are one of the most genetically diverse vertebrates, which positioned them well to evolve pollutant-tolerance.
However, it is important to realize that not all species are as well poised to adapt to ever-changing, increasing pollution in their habitats. Research like this gives us key information on how natural populations are dealing with our pollution. The best chance we can give all life forms on earth is to curb or rates of worldwide pollution. Luckily, it is on our power to do so!
References
In this article the authors describe an evolutionary convergence in mammals, birds, and reptiles, based on genomic data from NCBI. The evolution of different species and lineages is due to mutations that can appear and accumulate in organisms over time. Those mutations need a high functional potential and have to be conserved in time in order to form new species. The conservation of mutations can occur via selection pressure, mutational compensation, and/or by the separation of members from the same species by geological and environmental events.
In this comprehensive study, the authors describe, a genomic landscape of the parallel evolution by analysing functional nodal mutations (fNMs) by using different types of DNA (mitochondrial and nucleic), the thermostability of mtDNA encoding RNA genes, and the structural proximity of proteins, using the available 3D structures from PDB database. Functional nodal mutations (fNMs) can be separated in single nodal (fSNMs), recurrent nodal mutations (fRNMs), occured independently in unrelated lineages and recurrent combinations of nodal mutations (fRCNMs) recurred independently along with other nodal mutations in combinations in more than a single lineage. The recurrent ones can be taken in consideration the most when we are talking about the convergent adaptive responses, that means the parallel evolution of different species. In this study, one of the aim is to find the best candidate for this adaptive mutations that was present in the evolution of the amniotes. The compensated ones are used to identify the adaptive mutations. The main explanation for the convergent evolution is the presence of the recurrent nodal mutations. Many fNMs are in combination with potential compensatory mutations in RNA and protein-coding genes. The compensation of a functional mutation is the co-occurrence with additional mutations that are “affecting” the original function.
Results
In the article it is claimed that the evidence for parallel evolution is mainly due to the presence of a high number of uncompensated reccurent fNMs. The best candidate to show the parallel evolution is the emergence of body thermoregulation in mammals and birds, that seems to be independent.
The mtDNA, the maternal genetic information was used to identify the fNMs in the amniotes. The study is based on mtDNA from 1003 species and nDNA from 91 species. The mtDNA was used for the structure-base alignment for 24 mtDNA-encoded RNA genes (tRNA and rRNa) and 13 protein-coding gene. To this, they added 4 more mtDNA proteins with the 3D structure: CO1-3 and Cytb, as the cytochromes are highly conserved proteins across various species. The mtDNA genes are usually the same, but what seems to be different it is the order of the genes, that are changed by evolutionary rearrangements. Because of this, they first aligned the genes individually and after this, they concatenated the 37 proteins to the human mtDNA gene order.
The sequence alignment revealed a number of 25234 nodal non-synonimous and RNA gene mutations. To see the potential of this mutations, there were calculating a score that include: evolutionary conservation, physical-properties (of non-synonymous changes) and the molecular thermostability (the free estimated energy (?G) for the two RNA sequences was calculated before and after the mutational event). The score, from 1 to 9 is depending to the level of conservation and physico-chemical properties of the tested amino acid.After calculating the potential function score of all the nodal mutations, there were 3262 non-synonimous fNMs, mainly in RNA genes with mutations related to disease-causing.
The next step was to identify the best candidate for adaptive fNMs by studying the compensated and non-compensated mutations, but the approach chosen by the authors cannot reveal the exact order of compensation process. Meanwhile, there are some compensatory mutations that could gain lower functionality scores than the co-occurring fNMs. In the Figure 1, we can see a demonstration of the potential compensation and a possible adaptation in a protein-coding gene (COX2) through different species. The panel b shows the locations of the fNMs (S155T) and different other co-occurring compensatory mutations. The S155T mutation appears as independently re-occurrent as well as compensatory co-occurring mutations. As we can see, this approach is pure theoretical, because cannot show all the compensations, only the best ones, that got fixed in evolution. The Figure 2 shows the prevalence of different types of mutations that could be compensated or not. The predictive results reveal a high probability of fRCNMs to be compensated for RNA and protein-coding genes. Here are introduced also the information from the nDNA, that is compared with mtDNA in term of prevalence of the compensatory and non-compensatory mutations. Because there was a big difference of the number of species involved in this approach, the evolutionary resolution was reduced. So, the authors decided to analyze the same 91 species for mtDNA and nDNA and reducing the bias. Because of the reduction in the resolution, they redid the analysis by using the most ancient mutations, that occurs in deeper nodes in the case of mtDNA, but this revealed almost the same proccent as they were working with the 91 species (37% for the ancient mutations and 34% by including the younger ones) (Figure 2e & Supplementary 5b,c). So, the older mutations appear to be less compensated and this give more uncompensated mutations that are best candidates in the ancient adaptative mutations. In the supplementary Figures, the authors are using the OXPHOS complexes to compare the fNMs in mtDNA and nDNA by using 91 species. For the intra-mtDNA the albeit is less prominent (31%).
For the nDNA data is used the whole genome of the species. So, the information is much more comprehensive by the presence of a higher number of genes. In comparison with the mtDNA, the compensation prevalence is lower, having a difference of 10%, but in both case the proccent of possible compensation is higher than can be explained by the mutation rate or the chance.
In the end, to determine the best adaptive mutations over the evolution, they used the fRNMs from mtDNA, but maybe because of the low number of the samples, the result did not show any proof of the impact of non-compensated fRNMs in being the main reason for the convergent evolution. Instead, the nDNA revealed a significant pattern with highest number of potential non-compensated fRNMs shared between birds and mammals (N=51). The best candidates resulted by being the mutations in the genes related to the thermoregulation in the birds and mammals.
Conclusion
In this comprehensive study, the authors merged several information, including different types of DNA, from many species, with various physico-chemical parameters. The results of this work reveal, that the ancient functional mutation are the best for being studied, because of their possibility to overcome negative selective. The best candidates for the adaptive nodal mutations are in the end the non-compensated fNMs, that are in a higher presence in the case of old fNM. This seems to be the main helper for the evolution of the thermoregulation in birds and mammals. The protein analysis reinforces the main conclusion: for enriching the adaptative mutations, the non-compensated mutations are the best candidates.
Taken together this study provides new insights into how different lineages and species might have developed over time. It also shows a new way how to combine data from different sources. However, the authors fail in giving an adequate explanation for the fNMs, together with the fact that they lack references that describe this term makes the article difficult to understand, especially for people that are not from the field and this is in fact the contrary of how scientific writing should be done.
Levin & Mishmar, 2017, The genomic landscape of evolutionary convergence in mammals, birds and reptiles. Nature Ecology & Evolution 1: 0041
]]>
The highest genetic diversity in humans is found in Africa, in line with Africa being the cradle of humanity. While the three articles we discussed previously during this tutorial (1,2,3) mainly focused on determining the most parsimonious “out-of-Africa” scenarios based on genetic diversity data, this article (Skoglund et al. 2017 4) investigates the population structure of Africa prior to the expansion of food producers (i.e. herders and farmers). In order to reconstruct the prehistoric population structure, the authors analyzed the genomes from 16 ancient African individuals who lived up to 8100 years ago (including 15 newly sequenced genomes), as well as SNP genotypes from 584 present-day Africans, and 300 high coverage genomes from 142 worldwide populations. This is the first study to gather and analyze such a high number of ancient genomes, thereby providing an unpreceded insight into the prehistoric human population structure.
RESULTS
An ancient cline of southern and eastern African hunter-gatherers
The authors used principal component analysis (PCA) and automated clustering in order to relate the 16 ancient individuals to present-day sub-Saharan Africans. This reveals that while the two ancient South African individuals share ancestry with present-day South Africans (Khoe-San), 11 of the 12 ancient individuals living in eastern and south-central Africa between ?8100 and ?400 years BP form a gradient of relatedness to the eastern African Hadza on one extremity and to Khoe-San on the other. This genetic cline is also correlated with geography along a North-South axis. Another pattern which emerged from this analysis is the lack of heterogeneity between the seven ancient individuals from Malawi, indicating a long-standing and distinctive population in ancient Malawi which persisted for at least 5000 years but which is extinct today.
Subsequently, the authors built a model where ancient and present-day African population trace their ancestry to a putative set of nine ancestral populations. They then used data from both ancient and present-day populations showing substantial ancestry to major lineages present in Africa today as proxies for these ancestral populations. These proxy populations consisted of three ancient Near Eastern populations representative of Anatolia, the Levant and Irak, respectively, and six African populations representative of different components of ancestry (western African, southern African before agriculture, northeastern African before agriculture, central African rainforest hunter-gatherer, eastern African early pastoralist context and distinctive ancestry found in Nilotic speakers today). By using qpAdm (a generalization of f4 symmetry statistics), they tested for 1-, 2- or 3-source models and admixture proportions for all other ancient and present-day African populations, with a set of 10 non-African populations as outgroups. We note that the f4 statistics are poorly explained in this article, making it hard for a non-initiated reader to grasp its meaning and the relevance of the results. The main finding from this analysis is that ancestry closely related to the ancient southern Africans was present much farther north and east in the past than is apparent today.
Displacement of forager populations in eastern Africa
Unsupervised clustering and formal ancestry estimation both indicate that present-day Hadza in Tanzania can be modeled as deriving all their ancestry from a lineage related to ancient eastern Africans such as Ethiopia_4500BP. However the contribution of this lineage to present-day Bantu speakers in eastern Africans is small, who instead trace their ancestry to a lineage related to present-day western Africans and additional ancestry components. In present-day Malawians, population replacement by incoming food producers seems to have been almost complete as witnessed by a near absence of ancestry from the ancient individuals sampled, and by most of their ancestry coming from the Bantu expansion of western African origin.
Importantly, of all ancient individuals analyzed, only a 600 BP individual from Zanzibar has a genetic profile similar to present-day Bantu speakers, with even more western African ancestry. Using linkage disequilibrium, the authors estimate that the admixture between western- and eastern-African-related lineages occurred 800-400 years ago. This indicates that there was genetic isolation between early farmers and previously established foragers during the Bantu expansion into eastern Africa, and that this barrier disappeared over time as mixture occurred. However this delayed admixture did not occur in all African populations, as shown in present-day Malawians who display no signs of admixture from previously established hunter-gatherers.
Early Levantine farmer-related admixture in a ?3100-year-old pastoralist from Tanzania
The authors compared estimated the ancestry component from a 3100 BP individual from Tanzania and found that 38% of her ancestry was related to the pre-pottery farmers of the Levant (10000 BP), indicating a critical contribution of Levant-Neolithic-related populations to present-day eastern Africans. The best fitting ancestry component model in Somali indicates that they have ancestry from the 3100 BP Tanzanian individual but also Dinka-related ancestry as well as 16% ancestry related to Iranian-Neolithic-related ancestry. This suggests that ancestry related to the Iranian Neolithic appeared in eastern Africa after an earlier gene flow related to Levant Neolithic populations.
Direct evidence of migration bringing pastoralism to eastern and southern Africa
All three ancient southern Africans show affinities to the ancestry predominant in present-day Tuu speakers in the southern Kalahari. Among them, the 1200 BP sample from western Cape found in a pastoralist context has a similar ancestry composition as present-day pastoralists like the Nama, with affinity to three groups: Khoe-San, western Eurasians and eastern Africans. This is in line with the hypothesis of a non-Bantu-related population transporting eastern African and Levantine ancestry to southern Africa by at least 1200 BP. Using their model to determine the proportions of different ancestries present in western cape 1200 BP, they find mainly a mixture of non-southern African population. This is consistent with the hypothesis that the Savanna Pastoral Neolithic archaeological tradition in eastern Africa is a possible source for the spread of herding to southern Africa.
The earliest divergences among modern human populations
Previous studies indicate that the primary ancestry in the San population (southern Africa) comes from a lineage that separated from all other lineages present in modern humans, before separation of the different modern human lineages. While Skoglund et al. obtain a similar model in absence of admixture, the tree-like representation is a poor fit since ancient southern Africans (2000 BP) were not strictly an outgroup of all other African populations and several examples also show inconsistencies with this model. In order to find models that fit the data, the authors performed admixture graph modeling of the allele frequency correlations and found two parsimonious models. In the first one, present-day western Africans have ancestry from a basal African lineage that contributed more to the Mende than in did to the Yoruba, with the other source of western African ancestry being related to eastern Africans and non-Africans. In the second model, gene flow over long periods of time and over long distances has connected southern and eastern Africa to other groups in western Africa.
A selective sweep targeting a taste receptor locus in southern Africa
The authors then searched for the genomic signature of natural selection in ancient genomes, by searching for regions of greater allele frequency differentiation between ancient and present-day populations than predicted by the genome-wide background. To do this, the researchers compared the two ancient southern African genomes (2000 BP) to six present-day San genomes with minimal recent mixture. Since the small number of ancient genomes does not allow to infer changing allele frequencies at single loci, a scan for high allele frequency differentiation was conducted in 500 kb windows using 10kb steps. This led to the identification of the most differentiated locus which overlapped a cluster of eight taste-receptor genes. Although it is reported that taste receptors have already been identified as targets of natural selection as they affect the ability to detect poisonous compounds in plants, we must be wary that any analysis is bound to find something with such huge datasets, and that the biological interpretation of such finding may not be as straight-forward.
Polygenic adaptation
Skoglund et al. tested for evidence of selection on specific functional gene categories between present-day San and the two ancient genomes from southern Africa using allele frequency differentiation estimation. The functional category with the most extreme allele frequency differentiation between present-day San and the ancient southern Africans corresponded to “response to radiation”. In order to control that this was not a general inflated allele frequency differentiation, the same statistics were used using the Mbuti central African rainforest hunter-gatherer for which no enrichment for “response to radiation” was found. Instead, the top category for Mbutis was “response to growth”. Based on this, the authors speculate that the small stature of hunter-gatherer populations may be an acquired adaptation.
CONCLUSION
This study brings a first and unique view on the genetic makeup of prehistoric Africans. It is indeed a feat realized by 44 authors from institutions in 11 countries, which take advantage of 15 newly sequenced ancient genomes in addition to the only one that was previously available. The results indicate that an ancient lineage related to the San had a wider distribution in the past, depict two plausible scenarios of gene flow that led to the earliest divergences among modern populations and give new insights into the spread of herding and farming within Africa . As a side note, we noticed that all ancient individuals come from eastern or southern Africa, probably because this is where conditions were most favorable for the conservation of these ancient remains, although this could also introduce some biases, it seems to be the only possible way to go.
REFERENCES
Since the first genome of Bacteriophage MS21 was completely sequenced, in 1976, until 2001 when the first draft of human genome2 was released, a lot of work was done to improve and to make accessible different methods to get inside of the genetics of various organisms. For human genome, this step was a very important one and the Human Genome Project was declared complete in 20033. During the last years, more and more projects are involved in deciphering the human wanderlust. To all of previous studies, we can add The Simons Genome Diversity Project, that brought us more information by sequencing 300 new genomes from 142 diverse populations. One of the aim was to chose populations that differ in genetics, language and culture. The study shows that some of the populations separated 100000 years ago and reveals more information about the ancestors of Australian, New Guinean and Andamanese people.
Results
One of the most important thing in discovering the real human peopling of the Earth is to sequence as many as possible genomes, but from individuals coming from diverse populations, that could differ in many aspects. In this study, the 300 samples were prepared by using PCR-free library, through Illumina Ltd. method and the median coverage it was 42-fold (Figure S1.1; Supplementary Data Table 1). The method is using an improved genome coverage to identify the greatest number of variants with some of them previously reported. The single-sample genotypes was made by using the reference-bias free modification of GATK, but the some preprocessing was conducted for eliminating some adapter sequences. For increasing the data accuracy, it was used a filtering system, highly specific to the SGDP dataset. The levels are from 0 to 9 for each sample as a single character and the first level is the best for having a good balance between sensitivity and low error rate, but level 9 is good to be used when there is needed to low the errors rates (Figure S2.1).
The first part of the study is offering us more information about the time needed for the worldwide populations subjected to the study to get separated. The pairwise sequential Markovian coalescent (PSMC) and multiple sequentially Markovian coalescent (MSMC) was used to interpret the changes in size of the populations and the split time, the phased haplotypes of split time estimation were made by using the SHAPEIT and IMPUTE2. The filter used was the level 1. From the Figure 2a we can see evidence about the ancestors of some present populations that were isolated by at least 100kya, that could be an obstacle of certain mutations across the ancestors of all populations. The gene flow continued until around 50kya among the great majority of ancestral populations. The graphs show the moments when the substructure of different populations starts: in the Figure 2a, we can see that the substructure between french and africans start around 200 kya. In the next ones there is a comparison between only africans (the Yoruba separated from KhoeSan 87kya, from Mbuti 56kya and from the Dinka 19kya) or only non-Africans (the oldest substructure is from 50kya, taking part during or shortly after the deepest part of the shared non-African bottleneck 40-60kya). For the Figure 2d-f, it was used the PSMC and PS1 that show the effective population sizes inferred and the cross-coalescence rates inferred.
By using the neighbours-joining tree (pairwise divergence per nucleotide) and FST, Mallick et. al could reconfirm the previous studies regarding the fact that the deepest splits happened among the Africans. Previous studies showed that all non-Africans today possess Neanderthal ancestry and Figure 1c shows that the higher proportion of Neanderthal ancestry we can find it in East Asians. If we compare the EuroAsians between them, the South Asians have highest Denisovan ancestry (heatmap from Figure 1d). Another result is that there are more Denisovan ancestry in eastern than in western EuroAsians. If we take Australia, New Guinea or Oceania we can see that the results from other studies are confirmed by having more ancestry than in mainland Eurasians. In the Figure 3 the deeper the split is, the more divergent is the early dispersal ancestry. By using the cross-population coalescence pattern and allele frequency correlations, the best model is that the Australian, New-Guinean and Andamanese history doesn’t involve ancestry from an early- diverging source. In this study there is no archeological data taken in consideration regarding southern Asia or Australia. So, by using only the data from this study, it is released that the Australians, New Guineans and Andamanese are lacking in an analogous deep ancestry component. All the data referring to Australians seems to be consistent with descending
from a common homogeneous population since separation from New Guineans. Also, New Guineans, Australians and Andamanese appear as part of an eastern clade together with mainland EastAsians.
The 3P-CLR was used to scan the genome for positive selection. In the end, 38 of the largest peaks emerged for selection in the common ancestors of all modern humans. These peaks are the sweeps at the time that the archeological data shows an accelerated evidence of behavioral modernity. This data does not search for the sweeps on chromosome X or in repetitive or difficult-to-analyze sections of the genome.
For the rate of mutation accumulation between the non-Africans (grouped in America, CentralAsiaSiberia, EastAsia, WestEurasia, Oceania) and sub-Saharan Africans (grouped in Pygmy, Khoesan and Africa) it was supposed to be quite equal, but this study revealed an significant average of 0,5% difference. For this part, they used a highly restriction to the samples, by choosing only the samples processed in the same way and the highest level of filtering, pooling the samples from the same regions together. The one strength of this experiment is the fact that they avoid the bias due to different heterozygosity level in different populations (the heterozygosity is higher in Africans), by using only the chromosome X for males. Although, they map everything to chimpanzee, which is equally distant to all present populations. There are differences in observations related to other studies, by having a different rate of CCT>CTT mutation, that is close to Africans in Europeans, but not in East Asians. This could be explained by the decrease in generation interval in non-Africans since separation. Previous studies5 showed a higher X-to-autosome heterozygosity ratio in sub-Saharan Africans than in non-Africans. Mallick et al. confirmed this results by adding more populations to be analyzed: Khoesan for sud-Saharan Africans and New Guineans, Australians, Native Americans, Near Easternes and indigenous Siberians for the non-Africans. The only one exception, that showed a lower X-to-autosome heterozygosity ratio in sud-Saharan African than in non-Africans is in Pygmies (eastern Mbuti and western Biaka). In the Figure 1b through a scatterplot we can observe the two primary clusters: sud-Saharan Africans and all other populations, but without a big difference among the groups, except of the Pygmies with a high autosomal heterozygosity. If we compare the two Pygmies populations with a lower X-to-autosome ratio, we can see that the Mbuti are closer to non-Africans than to Africans, even if in the Neighbour-joining tree based on pairwise divergence, they are integrated to the Africans. The reduction of the X-to-autosome ratio in the non-African compared to African populations could be explained by the repeated waves of male mixture in already mixed population, but in the Pygmy populations, the strongest argument is the sex-biased gene flow supported by the anthropological data.
In the last part, Mallick et al. shows that the non-Africans are presenting a higher accumulation of mutations. This can be explained in two ways: the rate of mutations in non-Africans is increasing by acceleration of it or by a deceleration within Africans. The Extended Data Table 1, shows that none of the populations with strong signals of non-Africans could be in fact a deceleration of Africans. The acceleration in non-Africans could be caused by many possibilities: the life history traits (eg. generation interval) could change after the dispersal of modern humans outside of Africa, increasing the latitudes conquered by the humans or the colder climates, the gene conversion (GC to A or T alleles) was more effective in Africans or a Neanderthal admixture into the ancestors of non-Africans, that could accumulate more mutations than in the modern humans after separation (but there are not clear evidence about this fact).
Conclusion
The Simons Genome Diversity Project is bringing more information by studying 300 new genomes, from 142 diverse populations, that shows an acceleration of accumulation of mutations in non-Africans compared to Africans. Also, the Pygmies seem to be the only African group with a low X-to-Autosome diversity ratio. Regarding the ancestors, the highest proportion of Neanderthal it was present in EastAsians and an excess of Denisovan in some SouthAsians compared to other Euroasians.
Introduction
In the past two decades, considerable research effort has been made to sequence the human genome and subsequently trying to unveil the demographic history underlying the genetic patterns of diversity we observe today across the globe. Here we discuss a recent research article by Pagani et al. 1 that addresses genomic diversity and historic migration patterns of human populations in Eurasia. The first human genome was sequenced in 2003 by the Human Genome Project2 and larger projects rapidly followed, such as HAPMAP3 and the 1000 Genomes Project4, largely due to the considerable technological improvement of sequencing technologies. Despite being extremely useful tools for a number of studies, these genome databases have some important sampling caveats that limit their use to address some particular topics. Indeed, HAPMAP sampled a reduced number of populations whereas the 1000 Genomes sampled a large number of populations but did not attempt to sample individuals of “pure” ancestry. For instance, the sampling in North America focused considerably on city-based individuals that were found to have a very diverse recent ancestry thus blurring the signal of ancient colonisation history. Importantly, in the studied paper, a considerable effort was made on sampling a broad panel of 447 unrelated individuals of pure ancestry from 148 distinct populations, particularly including previously unstudied regions like Siberia and western Asia.
One of the main topics of the demographic history of humans that has long been of interest to researchers is the Out of Africa (OoA) of Anatomically Modern Humans (AMH) – a turning point in which humans dispersed from Africa and colonised Eurasia and ultimately Oceania and the Americas. Among other aspects, the number of OoA events has been the focus of discussion from which two major hypotheses emerged. The first, arguably the most wide-accepted, advocates for a single OoA event estimated at around 40 to 80 kya which gave origin to all extant non-african populations. The second hypothesis, dubbed the multiple-dispersal model5, considers multiple migration waves, more or less successful in settling in new continents, and possible admixture events between them at various points in time, which appears to be supported by previously described fossil evidence6,7. Interestingly, Tucci & Akey8 argue that these theories are not necessarily mutually exclusive but rather complementary as there could have been several failed or low-success OoA events followed by a major one that effectively colonised and subsisted in most continents.
In this study, Pagani et al. argue in favour of a multiple-dispersal scenario based on small remaining genetic contributions in the genomes of extant Papuans from an extinct lineage of AMH OoA earlier than the main OoA 75 kya.
Genetic structure and barriers across space
To obtain the first insight on the genetic structure among the sampled genomes, Pagani et al. employed two different approaches: first, treating SNP as independent markers (with ADMIXTURE9) and second taking into account linkage blocks (with fineSTRUCTURE10). Both strategies identified the major biogeographic groups of populations despite differences in resolution, defining 14 main genetic clusters across the globe (Extended Data Figure 1C). The detailed output from fineSTRUCTURE was interestingly used for a range of analyses from spatial patterns of genetic differentiation (Figure 1), co-ancestry (Extended Data Figure 3) and demographic history reconstruction (Extended Data Figure 7).
Taking advantage of their detailed sampling from Eurasia to Sahul, the authors employed a spatially explicit framework to study genetic differences and gene flow between populations as well as their association with environmental/geographic features at a large scale. Figure 1 illustrates this by representing the magnitude of the gradient of allele frequencies from SNPs across space, allowing to pinpoint the regions of major genetic gradients, i.e. potential barriers to gene flow, specifically mountain ranges, deserts and large water masses. These were consistent in broad strokes among the different analyses with the fineSTRUCTURE output (Figure S2.2.2-I) as well as the complementary migration-based EEMS (Estimating Effective Migration Surfaces; Extended Data Figure 5H). Importantly, the authors tested whether the geographic gaps in their sampling could bias the interpolation of barriers and showed their model remained robust in the face of new gaps (Extended Data Figure 5E-G).
On a second stage, Pagani et al. measured the association between the gradients of allele frequencies (termed as SNPs in Figure 1) and fineSTRUCTURE with three environmental barriers – elevation, temperature and precipitation – to determine the relative importance of the role each played in shaping the genetic patterns observed today. As one can see in the inset of Figure 1, SNPs indicated that elevation and precipitation had a strong spatial correlation with genetic differences whereas fineSTRUCTURE gave higher support to precipitation and temperature. This dissimilarity is likely due to the fact that the latter, as explained above, is dependent on linkage patterns. Linkage blocks are physical associations of loci that recombination renders temporary, unless they are specifically maintained by selection. Thus, current neutral linkage patterns reflect relatively recent demographic history, whereas the bulk of raw allelic frequencies reveal older patterns that influenced the majority of the genome. In the same sense, when taking into account only the rare variants (i.e. more recent), the association of SNPs with elevation was reduced (Figure S2.2.2-II).
The authors conclude these observations by suggesting that elevation contributed to shaping old migration routes (as confirmed by patterns of isolation by distance; Extended Data Figure 5A-C) but has not recently impeded the persistence of human populations. On the other hand, precipitation seems to be of paramount importance as populations continue to this day to avoid inhabiting low-precipitation regions such as deserts.
Despite the credibility of the conclusions, we raised some important questions on the analysis that could bias the interpretation. First, the authors did not address the innate correlation between the environmental variables (ex.: elevation and temperature) nor how or whether it was taken into account. Additionally, it is unclear which time period was used for temperature and precipitation as the study spans 120 thousand years of demographic history. Both these points could change the relative importance of a given variable, and should therefore have been specified clearly in the main text.
Selection screening
The authors scanned the genomes for evidence of purifying and positive selection through a series of different approaches and identified multiple candidate loci, some of which had been identified as targets of positive selection in previous studies. Additionally, the authors highlighted different levels of inter-population purifying selection, such as on olfactory receptor genes in Asians. Interestingly, they identified significantly stronger purifying selection in pigmentation and immune response genes in Africans than in the remaining populations, with the single exception of Papuans for the pigmentation genes (Extended Data Figure 6B). However, the authors did not discuss the possible factors behind such selective forces nor how this section on selection contributed to the main storyline and conclusions of the study.
Demographic history of Papuans
The results of fineSTRUCTURE were summarised with ChromoPainter and revealed very interesting patterns of haplotype co-ancestry and length as well as proportion of shared genome between populations. Leading is the observation that African populations display the highest co-ancestry (Extended Data Figure 3) and the shortest haplotypes (Figure S2.2.1-III), confirming their status as the oldest and most diverse populations. Short haplotypes reflect multiple recombination events through time indicating older ancestry. Thus, the most surprising observation was that Papuans have the shortest average haplotype length of all non-African populations (Figure S2.2.1-III), as well as the shortest African-inherited haplotypes (Extended Data Figure 7), which suggests an older ancestry with Africans than that of the remaining populations.
To investigate this further, the authors used multiple sequential Markovian coalescent (MSMC) to determine mean split times between genomes of Papuans and other populations, and it is represented in Figure 2A. This figure depicts the proportion of genome coalescing between populations over time (in logarithmic axis). However, it is important to take into account that for these calculations they used a generation time of 30 years, whereas the selection scans were done with a 25 years’ generation time. The latter is the most commonly used in the literature and no justification is given for this change. This analysis revealed an old split between the Papuan and African at about 90 kya (represented as Koinanbe in Figure 2A, red line), predating the split between Eurasian and African estimated at 75 kya (black line) and between Papuan and Eurasian at 40 kya (blue line). Despite the possible fluctuation in the absolute split times due to the chosen generation time, the relative differences between them is in line with Papuans harboring high amounts of short haplotypes, all suggesting an older population split than previously thought.
To explain the demographic history behind the observed patterns, the authors propose that a previously unknown admixture event took place in Sahul with either an archaic non-AMH (different from Denisovan and Neanderthal) or with a AMH resulted from an extinct OoA (xOoA). The latter hypothesis, which fits into the multiple-dispersal model explained earlier, would have taken place after the split of AMH with Neanderthal but before the main OoA.
Using coalescent simulations, the authors tried to replicate the split times by adding varying amounts of admixture with a non-AMH or with an AMH from a xOoA. There was no plausible scenario simulated of archaic admixture with non-AMH that could mirror the observed data. On the other hand, including in Papuans a genomic component that diverged from the main human lineage prior to the main OoA, replicated somewhat similar population split times. It is noteworthy that the main text indicates the “observed shift in the African-Papuan MSMC split curve can be qualitatively reproduced” under these conditions. In detail, it obtained a 3ky difference between the Papuan-African and Papuan-Eurasian splits (Figure S2.2.8-III) whereas the observed time-gap between the two is actually 15 kya (Figure 2A). The authors suggest that they may not be able to reach a comparable gap due to higher complexities of the demographic model that were not simulated within this study, such as population expansion and bottlenecks. Although this explanation appears reasonable, we believe it ought to have been made clear in the main text of the article.
To discern the weight of admixture with non-AMH, the authors masked putatively introgressed Denisovan haplotypes in Papuan genomes, which did not change the split times estimated between Papuans and the other populations (dashed lines in Figure 2A). Furthermore, the authors confirmed that MSMC behaved linearly through multiple events of admixture by studying populations with known admixture proportions in time (African Americans and Central and East Asians; Extended Data Figure 8), which allowed the calculation that the hypothesized xOoA would have split from most Africans around 120 kya (Supplementary Information 2.2.4).
On a supplementary line of examination, Pagani et al. looked at the age of African haplotypes in Papuans not present in other Eurasian populations by accessing the density of non-African alleles (nAAs) within them. The rationale behind this lies on the assumption that the rate of accumulation of nAAs, i.e. alleles not found among African genomes, within a haplotype of determined African origin in a non-African genome is proportional to the split date of that given population with Africans. First, this analysis revealed that Papuans had an overall higher amount of nAAs within African haplotypes along the genome than Eurasians (Figure 2B), indicating an older coalescent time with the Africans. Further, the proportions of nAAs within African haplotypes in Papuans were modeled under demographic scenarios of single and multiple-dispersal. The results showed that a xOoA of AMH that split around 120 kya from Africans was necessary to explain the constant elevated proportions of nAAs in Papuans (Figure 2D).
Combining results from the different approaches, the authors support an xOoA that split from Africans around 120 kya, and conclude by estimating it contributes to approximately 2% in contemporary Papuan genomes.
Conclusion
In this wide-ranging study, Pagani et al. discussed three main topics of human evolutionary biology in Eurasia using their extensive sampling: i) detect main geographic barriers to gene flow, ii) identify loci and ultimately pathways under selective pressure and iii) propose an extinct Out of Africa event earlier than 75 kya.
The latter was arguably the most important finding of this study with, as described above, the description of a 2% contribution in the genome of Papuans from an early xOoA. The authors provided multiple lines of compelling evidence pointing to an extinct Out of Africa expansion around 120 kya from Africans that admixed with the main OoA later in Sahul. The complete scenario is described in Extended Data Figure 10.
Nevertheless, the results presented in this paper and their associated methods are consistently poorly detailed and/or not self-explanatory. Such a paper covering a trendy topic in a high impact journal should be less indigestible for neophytes or even to fellow evolutionary biologists. Furthermore, the connection between the three main sets of analyses of the study (geographic barriers to gene flow, selection screening and the possibility of an xOoA) seems to be lacking as there is no global discussion bringing all points together.
Studied papers
Tucci & Akey 2016 Population genetics: A map of human wanderlust. Nature 538: 179–180
Pagani et al 2016 Genomic analyses inform on migration events during the peopling of Eurasia. Nature 538: 238–242
Reference
Exploration of variability of human genomes represents a key step in the holy grail of human genetics – to link genotypes with phenotypes, it also provides insights to human evolution and history. For this purpose Exome Aggregation Consortium (ExAC) have been founded; to capture variability of human exomes using next-generation sequencing. The first ExAC dataset of 63,358 individuals was released 20th of October 2014. Recently, a paper describing updated version of the dataset was published : Analysis of protein-coding genetic variation in 60,706 humans.
Authors made a great work on the reproductibility of the downstream analyses they have performed and generally on the availability of data. All the code is well documented in blogpost and available in GitHub repository. All figures in this blogpost I plotted by my own!
ExAC is composed of almost ten fold more individuals and previous dataset of the similar kind Fig 1a. 91,000 individuals were sequenced, of which 60,706 have been kept after quality filtering. Finnish population was excluded from European due to bottleneck they have gone though.
ExAC was targeting individuals with various genetic background. Principal component analysis have shown very strong geographical pattern in the dataset (Fig 1b). I expected a continuum of haplotypes in the environment without strong geographic obstacle (like European-Latino continuum). The gaps between South Asian samples and the rest Europen samples on the PCA plot is most likely caused by the absence of samples from Middle-East Asia. Middle-East Asian samples have just a colour, but no data points. Central Asians do not even have a colour.
Figure 1: Size and diversity of ExAC dataset a, ExAC dataset is almost ten fold bigger than datasets of similar kind: 1000 Genomes project and Exome Sequencing Project (ESP), but more importantly, it captures a far greater diversity of human populations compared to ESP and 1000 Genomes. b, The geographic signal of populations visualized using Principal component analysis (PCA). The first principal component get all the variability of African samples and it does not tells much about the rest of the dataset (Extended Data Figure 5 in the paper), therefore the second and third principal component has been show.
A 45 million nucleotide positions with sufficient coverage (>10x in at least 80% of individuals) are present in ExAC. These positions correspond to 18 million possible synonymous variants (in theory) of which ExAC is capturing 1.4 million (7.5%).
…mutational reoccurence: 43% of synonymous de Novo variants identified in previous studies were also identified in ExAC, which is a first direct evidence of mutational reoocuarence.
…multiple allels: 7.9% of high quality polymorphic sites are multiallelic, which is fairly close to Poisson expectation (whatever it means…)
…a LOT of variants after all the filtering, 7,404,909 high-quality variants were identified of which 317,381 indels. The density of variant is on the average one over eight bases. 99% of the variants had frequency bellow 1% and 54% of the variants are singletons (i.e. only one individual carries the variant).
…a selection effects The proportion of singletons among polymorphisms can serve as a measure of purifying selection acting on the polymorphisms of given size. The Figure 2 shows that indels that are not affecting open reading frame (ORF) have significantly less singleton variants than indels that actually affect ORF. There is also significant difference between indels of different sizes that are affecting ORF, but we (our topic group) have not found any possible explanation for this pattern.
…saturation of alleles in CpG sites: CpG sites have very high rate of transitions, therefore capturing all possible variants is substantially easier than for other sites. A subset of 20,000 individuals of ExAC dataset shows saturation of alleles – all non-lethal possible synonymous CpG transition variants are present. ExAC is the first dataset showing a saturation of human variation.
Figure 2: Indel frequencies with respect to the size a, Frequency of deletions is higher and smaller indels are more probable than greater. If we take into account the greater probability of smaller indels, frequency of indels that not shifting open reading frame is bit higher than frequency of indels than are not. b, Proportion of singletons in total number of indels (as proxy for strength of selection) is significantly and consistently lower in all indels that are not shifting open reading frame (-6, -3, +3, +6).
Authors introduce a mutability adjusted proportion singleton (MAPS) metric as a measure of selection. This metric is correcting on biases caused by the different mutational rates allowing comparisons of categories with various mutational speed. Comparison across different functional classes have shown at Figure 3. MAPS shows higher values for categories predicted to be deleterious by conservation-based methods.
Figure 3: MAPS values of different functional classes. MAPS is highest for nonense substiturions and it also consistent with PolyPhen and Combined Annotation Dependent Depletion (CADD) classification.
Average ExAC individual carries ~54 variants reported as Mendelian disease causing. Approximately 41 of these alleles were identified with frequency greater than one, therefore it is not expected to be caused by problem is variant calling, but in miss-classification of variants in the database. Evidence of 192 previously variants were manually curated of those only 9 had sufficient evidence in disease association. High allele frequencies were identified mainly in previously underrepresented categories Latino and South Asian.
ExAC have shown importance of matching reference population in identification disease-causing variant. An example is recessive disease North American Indian childhood cirrhosis previously reported to be caused by CIRH1A p.R565W. This variant was identified in homozygotic state in four individuals in Latino population, none of them having a record of liver problems during childhood.
ExAC shows the importance of diversity of sampled population in capturing the real link between genotype and phenotype. Even ExAC provides a lot of new insights, there are still populations that are underrepresented or not represented at all.
Given the richness of ExAC and the effort of authors in data sharing and availability, I guess that it will be a great resource for various analyses in the future for a lot of researchers around the globe.
Lek M, Karczewski KJ, Minikel EV, Samocha KE, Banks E, Fennell T, O’Donnell-Luria AH, Ware JS, Hill AJ, Cummings BB, Tukiainen T, Birnbaum DP, Kosmicki JA, Duncan LE, Estrada K, Zhao F, Zou J, Pierce-Hoffman E, Berghout J, Cooper DN, Deflaux N, DePristo M, Do R, Flannick J, Fromer M, Gauthier L, Goldstein J, Gupta N, Howrigan D, Kiezun A, Kurki MI, Moonshine AL, Natarajan P, Orozco L, Peloso GM, Poplin R, Rivas MA, Ruano-Rubio V, Rose SA, Ruderfer DM, Shakir K, Stenson PD, Stevens C, Thomas BP, Tiao G, Tusie-Luna MT, Weisburd B, Won HH, Yu D, Altshuler DM, Ardissino D, Boehnke M, Danesh J, Donnelly S, Elosua R, Florez JC, Gabriel SB, Getz G, Glatt SJ, Hultman CM, Kathiresan S, Laakso M, McCarroll S, McCarthy MI, McGovern D, McPherson R, Neale BM, Palotie A, Purcell SM, Saleheen D, Scharf JM, Sklar P, Sullivan PF, Tuomilehto J, Tsuang MT, Watkins HC, Wilson JG, Daly MJ, MacArthur DG, & Exome Aggregation Consortium. (2016). Analysis of protein-coding genetic variation in 60,706 humans. Nature, 536 (7616), 285-91 PMID: 27535533
]]>High throughput genotyping and sequencing has led to the discovery of numerous sequence variants associated to human traits and diseases. An important type of variants involved are Loss of Function (LoF) mutations (frameshift indels, stop-gain and essential sites variants), which are predicted to completely disrupt the function of protein-coding genes. In case of Mendelian recessive diseases, for the condition to occur, the LoF variants must be biallelic, i.e. affecting both copies of a gene. The affected gene is then defined as “knockout”.
By studying the Icelandic population, authors aim to identify rare LoF mutations (Minor Allele Frequency, MAF < 2%) present in individuals participating in various disease projects. They then investigate at which frequency in the population these LoF mutations are homozygous (i.e. knockout) in the germline genome.
The Icelandic population Iceland is well-suited for genetic studies for three main reasons. The island was colonized by human population around the 9th century by 8-20 thousand settlers. Since then the population grew to around 320’000 inhabitants today. The initial founder effect and rare genetic admixture make the Icelandic population a genetic isolate. In addition to an unusual genetic isolation, Iceland’s population benefits of a genealogical database containing family histories reaching centuries back in time, as well as a broad access to nationwide healthcare information.
These characteristics led to the development of large-scale genomic studies of Icelanders by deCODE Genetics. This biopharmaceutical company has published various studies, including this paper, related to genetic variants and diseases in Icelanders.
Loss of function mutation and rare complete knockouts Authors sequenced the whole genome of 2’626 Icelanders participating in various disease projects and identified variants in protein coding genes. These variants were annotated with the predicted impact that they have on the gene: LoF, moderate or low impact. A total of 6’795 LoF mutations in 4’924 genes were identified, with most of these variants (6’285) being rare (MAF < 2%).
The identified LoF variants were imputed into an additional 101’584 chip-genotyped and phased Icelanders, allowing the identification of the number of knockout genes in the population. Authors found that 1’485 previously identified LoF mutations (MAF <2%) are contributing to the knockout of 1’171 genes and that 8’041 individuals possess at least 1 of these knockout genes. Out of these 1’171 genes, 88 had been already linked by previous studies to conditions through a recessive mode of inheritance.
Double transmission deficit of LoF variants Because knockout genes should be deleterious for an organisms, we expect a deficit of homozygous for these genes in the population due to embryonic/fetal, perinatal or juvenile lethality. To investigate whether such a deficit was present, authors calculated the transmission probability of LoF variants from parents to their offspring.
Under Mendelian inheritance, the expected percent of transmission of the LoF mutated gene from heterozygous parents to their offspring (i.e. double transmission) is of 25%. However, results show a statistically significant deficit in double transmission, the observed double transmission probability being of 23.6%.
The rare LoF mutations were ranked according to the Residual Variation Intolerance Score (RVIS) percentiles and essentiality score percentiles. Both measures attempt to classify genes according to their tolerance to functional variation, with the lowest rank corresponding to genes being more sensitive to mutations. As expected, the lowest double transmission rate was found for the most sensitive genes (first percentile), suggesting that a homozygous state of LoF mutation in these genes is deleterious.
Tissue specific expression of knockout genes Authors investigated if genes were more likely to be knockout when expressed in specific tissues. By retrieving the information from previous studies of the number of genes that are highly expressed in 1 or more – but not all – 27 tissues, they calculated the fraction of these genes that were knockout in each tissue. They found that the brain and placenta were the tissue with the lowest fraction of knockout genes (3.1% and 3.9%, respectively), and that in testis, small intestine and duodenum were observed the highest fraction of biallelic LoF mutations (5.8%, 6.4%, and 6.9% respectively).
Conclusion and Comments The characteristics of Icelandic population and the incredibly large sample size (~ 1/3 of the total population) allowed authors to identify a large number of new and rare LoF mutations. Part of these mutations was shown to contribute to the knockout of an unexpected large number of genes in an unexpected large number of people. This study is the first to shed a light on the astonishing number of knockout present in human populations. In addition, by investigating the transmission probability, a deficit in homozygous loss-of function offspring was identified, especially when LoF mutations affected essential genes. This result was expected because of the predicted deleterious effect of biallelic LoF mutations.
Besides the aforementioned interesting results of the paper, some aspects were slightly disappointing. First, I was expecting authors to focus more on the genotype-phenotype aspects. Even if they pinpoint a deficit in double transmission, suggesting deleterious consequences for the organism, authors did not discuss the function of the identified knockout genes and their effect on the phenotype. Second, the paper was not an easy read. Many results were only mentioned without additional information on the methods or data used, and it was sometimes difficult to link them with the main aim of the study. Additionally, figures were sometimes misleading because of different axis scales or incomplete legends.
Finally, authors suggested that important tissues, such as the brain, have a lesser number of knockout compared to other tissues, writing that “genes that are highly expressed in the brain are less often completely knocked out than other genes”. However, this result is questionable as we do not have any measure of the number of knockout genes that we expect to be expressed only by chance in the tissues. In other words, the brain could have a lower number of knockout genes expressed compared to other tissues only because the total number of expressed genes in the brain is lower. Therefore we do not know if the lower number of knockout genes in the brain is due to chance or to biological reasons.
Nevertheless, this study opens the door to understanding how many knockout genes occur without phenotypic consequences in humans, what are the genes function and essentiality, and the role of the environment in the buildup of phenotype. The classical search for genetic variants associated to a phenotype, as in GWAS studies, could be reversed by first identifying individuals with the same genetic variants and then precisely phenotyping them.
Sulem, P., Helgason, H., Oddson, A., Stefansson, H., Gudjonsson, S., Zink, F., Hjartarson, E., Sigurdsson, G., Jonasdottir, A., Jonasdottir, A., Sigurdsson, A., Magnusson, O., Kong, A., Helgason, A., Holm, H., Thorsteinsdottir, U., Masson, G., Gudbjartsson, D., & Stefansson, K. (2015). Identification of a large set of rare complete human knockouts Nature Genetics, 47 (5), 448-452 DOI: 10.1038/ng.3243
]]>
Cindy Dupuis, Xinji Li, Casper van der Kooi
The development of new molecular mechanisms and next generation sequencing techniques have advanced our knowledge on the genetic basis underlying phenotypic polymorphism. Over the coarse of recent years, scientific studies have documented large genomic regions with drastic phenotypic effects, the so-called supergenes. A supergene is a set of genes on the same chromosome that exhibit close genetic linkage and thus inherits as one unit.
The evolution of a supergene requires that multiple loci with complementary effects become linked (i.e. they are genetically clustered and recombination between the loci is suppressed) and that optimal alleles at the linked loci are combined. Genetic clustering of different loci can occur when, via mutation, an adaptive interaction between two closely placed loci is created. In addition, gene duplications or translocations that generate a series of (novel) complementary genes can give rise to supergenes. The probability of a recombination event occurring in between loci depends on various factors. The chance of a recombination event occurring in between two loci will be small when the loci are located closely together, as the chance of a recombination event in between two loci generally decreases with physical distance between the loci. Given the large size of supergenes, additional mechanisms seem, nonetheless, important. This can, for instance, be maintained via structural differences, such as inversions, between the supergene and their homologous chromosomal region.
An interesting example of a supergene in an invertebrate is the case documented by Purcell et al. (2014). They documented a large, nonrecombining region that is association with social organisation in an ant species. The nonrecombining region was found to largely constitute one chromosome and was hence aptly called the ‘social chromosome’. They find a structurally similar region with similar effects in another ant species, however the regions exhibit no homology, suggesting parallel evolution of the social chromosome. Examples of vertebrates social systems determined by supergenes are, to our knowledge, unknown.
Two recent articles (Küpper et al., 2016; Lamichhancy et al., 2016) revealed a single supergene controlling alternative male mating tactics in the ruff (Philomachus pugnax). The studies were carried out independently by two research groups, but reach almost the same conclusions. The ruff (Philomachus pugnax) is a lekking wader known for the great diversity in the male plumage color and behavioral polymorphism. Three types of males can be distinguished; these types are characterized by differences in territoriality and behavior that are highly correlated with differences in nuptial plumage and body size. Predominantly dark-colored Independent males are most common (80-95% of males), these males defend small territories on a lek. Smaller, lighter colored Satellite males (5-20%) are non-territorial and less strict to a particular lek. Satellite males make use of – and are largely tolerated by – the residences of Independent males. The third type are the Faeder males, which are very rare (<1% of males). Faeder males lack male display, are small and resemble the unornamented females; however, they have disproportionately large testes.
Previous studies using pedigrees of large, captive populations showed that reproductive polymorphism follows a single-locus autosomal pattern of inheritance (Lank et al., 1995; Lank et al., 2013). The dominant Faeder allele controls development into Faeder males, whereas the Satelllite allele (that is dominant to Independent) controls development into Satellite or Independent males. Ekblom et al. (2012) studied the nucleotide sequence variation and gene expression in ornamental feathers from 5 Independent and 6 Satellites males using transcriptome sequencing. No significant expression divergence of pre-identified coloration candidate genes was found, but many genetic markers showed nucleotide differentiation between the two morphs. Later, Farrell et al. (2013) used linkage analysis and comparative mapping to locate the Faeder locus, and found linkage to microsatellite markers on avian chromosome 11 that included the Melanocortin-1 receptor (MC1R) gene, a strong candidate in alternative male morph determination, because it is considered to be important in plumage coloration.
Using the captive population that was previously phenotyped, Küpper et al. now set out to determine the genomic structure of the existing morph divergence in P. pugnax. The first step in their analysis was to generate and annotate the full genome for one Independent male. Followingly, the authors identified SNPs in the population using RAD sequencing. More than one million SNPs could be distinguished, and Faeder and Satellites could be mapped to a genetic map based on 3’948 SNPs. Interestingly, both morphs mapped to the same region on chromosome 11, but exhibited clear structural differences. This was corroborated by a GWAS analysis on 41 unrelated Satellite, Independant and Faeder males from a natural population.
In order to characterize the genomic region more precisely, they conducted a whole genome sequencing of a small set of Independent, Satellite and Faeder males. They showed that the region on chromosome 11 was highly differentiated between Satellite and Faeder morphs and that this region contained a greater nucleotide variation compared to the adjacent regions. Using the reads orientation, they found clear evidence for an inversion of the chromosomal regions between the different morphs. Interestingly, they found that one breakpoint occurs within an essential gene, CENPN (encoding centromere protein N, recessive lethal), which implies that individuals homozygous for the inversion are not viable – an observation that is confirmed by breeding experiments. The authors also suggested a recombination event or gene conversion to have occurred between the Satellites and Independent alleles.
By comparing gene sequences among morphs, the authors discovered that 78% of the gene sequences were different between morphs, and that those differences had the potential to change the encoded protein. Among the divergent genes, some where found to be involved in hormonal production, like HSD17B2, an enzyme inactivating testosterone and estradiol. Varying specifically depending on the morph, this enzyme may alter steroid metabolism and explain partly why plumage patterns and behavior is different between morphs. The MC1R gene was also found within the altered genomic region. This gene is considered an important locus controlling color polymorphism, which could be at the source of the reduced melanin levels in satellites. The PLCG2 gene, which has been rearranged in Faeders, was found to be a candidate gene for the rather feminine appearance and non-aggressive behavior in Faeders. Presumably, this gene is part of a cascade leading to the development of the usual impressive plumage of other males morphs.
In a second article, Lamichhancy et al., 2016 studied a natural ruff population using whole-genome sequencing. They first established a high-quality reference genome assembly from an Independent male and conducted functional annotation based on both evidence data and de novo gene predictions. Then, whole-genome resequencing and SNP calling were performed for 15 Independent, 9 Satellite and 1 Faeder males. Their genome-wide screen for genetic divergence estimates (FST) between different male morphs identified a 4.5-Mb region, based on which Independents and Satellites could be phylogenetically clustered as distinct groups. Screening for structural variants identified a 4.5-Mb inversion in Satellites that perfectly overlapped with the differentiated region. In addition, PCR-based sequencing confirmed the positions of proximal and distal breakpoints and identified a 2,108-bp insertion of a repetitive sequence at the distal breakpoint. Diagnostic tests showed that Satellite males were heterozygous (S/I), while most Independent males were homozygous (I/I). They suggested the Independent allele to represent the ancestral state, which is consistent with the conserved synteny among birds.
The comparison between Faeder and Independent males showed that the genetic differentiation was equally strong across the same region, creating a mirror image of the differentiation pattern between Satellites and Independents. Accordingly, the region could be subdivided into two parts: region A where Satellite and Faeder chromosomes were closely related and less closely related to Independent, and region B where the Satellite and Independent loci were closer related and divergent from Faeder. Since an inversion is expected to reduce the amount of recombination within the region between the wild-type (I) and mutant alleles (either S or F), the disruption of the differentiation pattern might be considered the result of one or two recombination events between an Independent and a Faeder-like chromosome. The divergence time between the Independent allele and Satellite or Faeder alleles was estimated to be approximately 4 million years, using the nucleotide divergence and estimated mutation rates for birds. The last recombination event was estimated to occur 520,000 ± 20,000 years ago.
To better understand the genetic consequences of the inversion and relate it to the phenotypic variantion in male ruffs, the authors searched for candidate mutations amongst the genes in the inverted region. Mutations in several genes with important functions were found on Satellite and Faeder chromosomes, including the abovementioned CENPN, HSD17B2 and MC1R genes as well as and SDR42E1 (the latter one is important for the metabolism of sex hormones). Missense mutations in derived MC1R were found to be associated to the Satellite and Faeder alleles, hinting at a potential mechanism explaining the male plumage polymorphism during breeding season.
In conclusion, these two studies demonstrated presence of a genomic inversion that led to the evolution of a supergene. This supergene determines the complex phenotypic variation in male ruffs. These two papers contribute to our understanding of supergenes, complex phenotypes and social organization.
Küpper C, Stocks M, Risse JE, Dos Remedios N, Farrell LL, McRae SB, Morgan TC, Karlionova N, Pinchuk P, Verkuil YI, Kitaysky AS, Wingfield JC, Piersma T, Zeng K, Slate J, Blaxter M, Lank DB, & Burke T (2016). A supergene determines highly divergent male reproductive morphs in the ruff. Nature genetics, 48 (1), 79-83 PMID: 26569125
]]>Understanding the evolutionary history of our own species, how migration and mixture of ancestral populations have shaped modern human populations is a key question in evolutionary biology. Here we present three articles related to this topic, the first two dealing with India and the third one focusing on a single Ethiopian group :
1) Moorjani et al 2013 Genetic Evidence for Recent Population Mixture in India AJHG 93,: 422–438
2) Basu et al 2016 Genomic reconstruction of the history of extant populations of India reveals five distinct ancestral components and a complex structure PNAS online before print
3) Van Dorp et al 2016 Evidence for a Common Origin of Blacksmiths and Cultivators in the Ethiopian Ari within the Last 4500 Years: Lessons for Clustering-Based Inference PLOS Genetics 11(8): e1005397
All of them use genome wide data from micro array. After a brief abstract of each paper, showing their similarities and differences, we discuss their methodological approaches.
Ancestral populations of India
The aim of the first two articles is to understand the history of the populations of the Indian subcontinent. The first one (Moorjani et al 2013) reports data from 73 groups living in India for more than 570 individuals sampled. The authors filtered out the data by removing all individuals with evidence of recent admixture or recent ancestry from out of India. The populations that were included in the analysis can be classified into two linguistic categories: the ones speaking Indo-European languages and the ones speaking Dravidian languages.
Previous genetic evidence indicates that most of the groups of India descend from a mixture of two distinct ancestral populations: Ancestral North Indians (ANI) and Ancestral South Indians (ASI). Three different hypothesis exist for the date of mixture of these two populations:
1) arrival of ANI is due to migration prior to agriculture about 30,000-40,000 years ago
2) ANI arrived with the spread of agriculture who probably began around 8,000 and 9,000 years ago
3) ANI arrived very recently (3,000-4,000 years ago) when the Indo-European languages presumably began to be spoken in India.
To prove the admixed origin of Indian groups and estimate the proportion of each ancestry in each population they use a PCA and a statistic called F4 ratio that infers the mixture proportion measuring the correlation in allele frequencies between each pair of groups. They demonstrated that all populations are admixed and lie along an “Indian cline”, that is a gradient going from 17% of ANI ancestry to 71%. These results correlate well with geography and language, with the northern Indo-European populations having more ANI ancestry than the southern Dravidian ones. Then they use linkage disequilibrium (LD) to estimate the dates of admixture : LD blocs are longer if the admixture is younger. By fitting an exponential function to the decay of LD (that is expected from a sudden cessation of admixture) they could estimate that admixture occurred between 1,856 and 4,176 years ago, supporting the third hypothesis. These results correspond with demographic and cultural changes observed in India with the establishment of the caste system leading to strong endogamy that stopped the admixture rapidly. Moreover they found that Indo-Europeans groups have more recent admixture dates, which could be explained by multiple waves of mixture in these populations. Another finding of this paper is that aboriginal Andaman Islanders (Onge) belong to a sister group of ASI.
The second article (Basu et al 2016) has the same focus region and use the same basic dataset, except that the authors kept the all populations in the analyses, including the austro asiatic (AA) and tibeto burman (TB) speakers. They first ran ADMIXTURE on all populations and showed that islanders and mainland populations have distinct ancestral components (islanders share ancestry with oceanic peoples like Papuans). In a second time they ran the same analysis on mainland populations only (thus excluding population from the Andaman and Nicobar islands). The best model was composed of four ancestral components, the ANI, the ASI as well as the ancestral AA and TB and they found that several present day populations are almost pure representatives of these ancestral components (figure 2).
They further estimated the time and extent of admixture using the degree of fragmentation (due to recombination) of haplotypes blocs originating from a donor population into the recipient population. In each population, the distribution fitted again with an exponential curve. They showed that admixture abruptly came to an end about 1575 years ago in upper-caste populations, most likely due to the establishment of endogamy, while tribal populations seemed to have admixed until 1500-1000 years ago.
In short, although they share a common topic, these two papers propose divergent versions of the history of Indian population : while the first considers a priori that austro asiatic and tibeto burman speakers are not component of the ancestral populations of India and only focuses on the mixture between the ANI and ASI components, the second paper claims that the genetic structure of Indian population is the result of admixture events between four ancestral components. However the two views converge on the idea that admixture was a common phenomenon in India that ceased rapidly with the establishment of the caste systems that enforced endogamy.
Common origin of two subgroups of Ari people
The 3rd paper investigates the history of human populations at a smaller scale, focusing on a single ethnic group, the Ari people of Ethiopia. The Ari are composed of two socially and genetically distinct subgroups : the cultivators (Aric) and the blacksmiths (Arib). Anthropologists have proposed two alternatives hypothesis to explain the division of the Ari : under the remnant hypothesis (RN), the blacksmiths are the remnants of an indigenous group that was assimilated by the more recently arrived cultivators, whereas the marginalization (MA) hypothesis proposes that the two groups share a common ancestry but the blacksmith were recently marginalized due to their activity. While anthropologists traditionally favour the MA hypothesis, recent genetic studies have provided support for the RN hypothesis. In this article the authors use a new methodology on the same genetic dataset to bring evidence for the MA hypothesis. They show that when ADMIXTURE, fineSTRUCTURE or CHROMOPAINTER analysis are run on a complete dataset of 237 samples of 12 Ethiopian and neighbouring populations, the Arib are grouped into a single homogeneous cluster. But when the patterns of haplotype sharing are inferred by composing the Ari as a genetic mixture of all other groups, except themselves, the genetic differences between Arib and Aric disappear. In fact, their analyses reveal that the two Ari groups have the same mixture events with non Ari populations (figure 3).
To explain this pattern they propose that the genetic differentiation of the blacksmith is due to a bottleneck effect. Their hypothesis is supported by the fact that identity-by-descent (IBD) is stronger in blacksmiths than cultivators which is consistent with reduced genetic diversity in the blacksmiths. Using the D-statistic, they also show that the Arib and Aric are more closely related to each other than they are to any other Ethiopian group. Therefore they conclude that the observed genetic differentiation between the Arib and Aric does not represent separate ancestry but is rather the result of strong genetic drift due to a bottleneck effect induced by the social marginalization of the blacksmiths.
Methodological discussion
What stands out from reading these three articles is that selection of a proper methodology is crucial within an hypothesis testing framework. While the two articles on Indian populations use the same initial dataset, the way they filter and analyse it results in very different conclusions. The inclusion or exclusion of some populations from an admixture analysis or outgroup selection for an f4 ratio estimation directly impact the output of these analysis and can lead the authors to tell very different stories. Before disclaiming or putting forward one hypothesis, it is important to be aware of the limitations of the method that is used to produce the results. For example the authors of the second paper on India’s ancestral populations, claim to demonstrate a more complex history than shown in the first paper but their result is solely based on a clustering analyse (implemented in various softwares such as STRUCTURE or ADMIXTURE).
The basic principle of those STRUCTURE/ADMIXTURE like programs is to take the K most different groups of the dataset, consider them as the pure ancestral groups and force the others to be a combination of those. This means that the results depend on the populations and the number of clusters K that are input in the program. There are different methods to determine which K provide the best fit to the data (cross-validation error, delta K …) but in numerous cases the inferred mixture proportions are wrong. Only in very simple cases, like the African American genetic history (well explained in Daniel Falush’s blog) that involves three clearly defined and very differentiated ancestral populations (West Africans, Europeans and Native Americans) we can be confident in the results of the clustering analyse.
But in many cases the history is more complex and no current population actually corresponds to a pure ancestral population because of multiple waves of admixtures. In this case the most differentiated groups correspond only to the most extreme groups but it does not mean that these groups are pure or ancestral. This is well explained in Razib Khan’s blog using the simple example of Uygurs and Europeans : it is known that the Uygurs are a recently mixed group (between European and Asian) but if K is fixed to 2 with Uygurs and Europeans, STRUCTURE will form two different clusters at 100% levels, one with the Uygurs and one with Europeans. This is why, in the 2nd paper, the apparently pure AAA, ATB, ASI and ANI populations and all the clustering implications are probably meaningless. In fact, when using the f4 ratio (as in the first paper) all groups are found to be admixed to a certain extent (with the smallest rate of admixture being 17%).
This critic of clustering analysis is a key element of the study on the Ari people where the authors point out that results from such methods should not be taken for granted but interpreted with caution. Indeed this kind of method cannot discriminate between alternative scenarios of recent mixture of separate populations or shared ancestry followed by population divergence. Therefore support for one of these hypotheses should rely on additional tests. Instead of directly accepting the story suggested by a clustering analysis, a more reasonable work-flow would be to use other methods in order to address the specific implications of one hypothesis. This is exactly what is done in the third article where, as we previously explained, the authors constrain the analysis of mixture by forbidding self ancestry in the two groups of interest which remove the confounding effect of recent bottleneck. In such complex cases, associating PCA and STRUCTURE-like analyses with F-statistics and simulations allow to draw a more robust conclusion. Indeed statistics such as Fst or Dxy that estimate the genetic differentiation between two populations can be simulated under alternative scenarios, representing competing hypothesis (figure 5). These simulated statistics can be subsequently compared with the ones estimated from real data to favour one hypothesis over the other. Simulations can also give an idea of how difficult it is to discriminate between the different hypothesis, which avoid over interpretation of the results. In the second paper, where the authors put forward an new hypothesis, radically different from the classical hypothesis of anthropology and other genetic studies, additional tests like these seem necessary to strengthen their conclusions.
Although it was not mentioned in any of the articles, the quality of the data and the way to obtain them, i.e. the kind of sequencing methodology, should also be a matter of precaution. Indeed, they all use micro arrays designed from European populations. These micro arrays consist of thousands of DNA spots containing a predefined sequence, known to be polymorphic in Europeans and only the complementary sequence can fix to this spot and be sequenced. So using these micro arrays to study the history of non european populations may be problematic as only SNPs that are variable for europeans will be targeted, probably leading to the exclusion of meaningful information for non European populations. Today, with New Generation Sequencing (NGS) there are many alternatives, such as RAD sequencing or Whole Genome Sequencing, that allow to sequence tens of thousands non-predefined SNPs.
Conclusion
To conclude, the take home messages from these three articles are :
– Social systems leading to endogamy can influence and modify rapidly and dramatically the genetic structure and patterns of humans populations.
– It is difficult to reconstruct the ancestry of human populations, especially when they involve a complex process with multiple waves of admixture.
– Clustering methods are designed to find a structure in a genetic dataset but they do not necessarily reflect real shared ancestry. Further test using other methods are required to robustly support one hypothesis.
]]>Introduction
Darwin’s finches from Galapagos and Cocos Island are classic example of young adaptive radiation, entirely intact because none of the species having become extinct as a result of human activity. They have diversified in beak sizes and shapes, feeding habits and diets in adapting to different food resources. Although traditional taxonomy of Darwin’s is based on morphology and has been largely supported by observations of breeding birds finches, in this paper, authors showed the results of whole-genome re-sequencing of 120 individuals representing all of the Darwin’s finch species inhabiting Galapagos archipelago (Fig. 1a) and two close relatives, trying to analyse patterns of intra-and interspecific genome diversity and phylogenetic relationships among the species.
Figure 1a. Sample location of Darwin’s finches
Summary and comments of the paper
The authors analyzed location and phylogeny of Darwin’s finches and found widespread evidence of interspecific gene flow that may have enhanced evolutionary diversification throughout phylogeny. They also reported discovery of a locus with the major effect on beak shape. They generated 10x sequence coverage per individual bird and using 2×100 base-pair (bp) paired-end reads and found evidence of introgression from three sources: ABBA-BABA tests, discrepancies between phylogenetic trees based on autosomal and sex linked loc, and mtDNA. Extensive sharing of genetic variation among populations was evident, particularly among ground and tree finches, with almost no fixed differences between species in each group. Their maximum-likelihood phylogenetic tree based on autosomal genome sequences is generally consistent with current taxonomy showing several interesting deviations (Fig. 1b).
Figure 1b. Phylogeny of Darwin’s finches
Revised and dated phylogeny of Darwin’s finches shows that the adaptive radiation took place in the past million years, with a rapid accumulation of species recently. Genomic characterization of the entire radiation revealed a striking connection between past and present evolution. Evidence of introgressive hybridization is found throughout the radiation, showing that hybridization always gives rise to species of mixed ancestry, which is explained in detail (species and location) in this paper. The most obvious morphological difference among Darwin’s finches concerns beak shape. The authors performed a genome wide scan on the basis of populations that are closely related but show different beak morphology. In this study, they indicated a polygenic basis for beak diversity, discovering 15 regions with strong genetic differentiation between groups of finches with blunt or pointed beaks. Their analysis revealed that ALX homeobox 1 is an excellent candidate for variation in beak morphology, because it encodes a paired-type homeodomain protein (transcription factor), that plays a crucial role in development of structures derived from craniofacial mesenchyme, the first branchial arch and the limb bud, and have influence on migration of cranial neural crest cells, highly relevant to beak development. They observed single nucleotide polymorphisms (SNPs) in ALX1 gene of various finch species and concluded that blunt haplotype has a long evolutionary history because it’s origin predates the radiation of vegetarian, tree and ground finches. The haplotype might have evolved by accumulating both coding and regulatory changes affecting ALX1 function. Natural selection and introgression affecting this locus have contributed to the diversification of beak shapes among Darwin’s finches and hence to their expanded utilization of food resources on different Galapagos islands.
Lamichhaney, S., Berglund, J., Almén, M., Maqbool, K., Grabherr, M., Martinez-Barrio, A., Promerová, M., Rubin, C., Wang, C., Zamani, N., Grant, B., Grant, P., Webster, M., & Andersson, L. (2015). Evolution of Darwin’s finches and their beaks revealed by genome sequencing Nature, 518 (7539), 371-375 DOI: 10.1038/nature14181
]]>