1. Evolutionary Biology
  2. Genetics and Genomics

Intra-species differences in population size shape life history and genome evolution

  1. David Willemsen
  2. Rongfeng Cui
  3. Martin Reichard
  4. Dario Riccardo Valenzano  Is a corresponding author
  1. Max Planck Institute for Biology of Ageing, Germany
  2. Czech Academy of Sciences, Institute of Vertebrate Biology, Czech Republic
  3. Department of Botany and Zoology, Faculty of Science, Masaryk University, Czech Republic
  4. CECAD, University of Cologne, Germany
https://doi.org/10.7554/eLife.55794
Share this article
Cite this article
  1. David Willemsen
  2. Rongfeng Cui
  3. Martin Reichard
  4. Dario Riccardo Valenzano
(2020)
Intra-species differences in population size shape life history and genome evolution
eLife 9:e55794.
https://doi.org/10.7554/eLife.55794
  2. Download BibTeX
  3. Download .RIS

Abstract

The evolutionary forces shaping life history divergence within species are largely unknown. Turquoise killifish display differences in lifespan among wild populations, representing an ideal natural experiment in evolution and diversification of life history. By combining genome sequencing and population genetics, we investigate the evolutionary forces shaping lifespan among wild turquoise killifish populations. We generate an improved reference genome assembly and identify genes under positive and purifying selection, as well as those evolving neutrally. Short-lived populations from the outer margin of the species range have small population size and accumulate deleterious mutations in genes significantly enriched in the WNT signaling pathway, neurodegeneration, cancer and the mTOR pathway. We propose that limited population size due to habitat fragmentation and repeated population bottlenecks, by increasing the genome-wide mutation load, exacerbates the effects of mutation accumulation and cumulatively contribute to the short adult lifespan.

Introduction

The extent to which drift and selection shape life history trait evolution across species in nature is a fundamental question in evolutionary biology. Variations in population size among natural populations is expected to affect the rate of accumulation of advantageous and slightly deleterious gene variants, hence impacting the relative contribution of selection and drift to genetic polymorphisms (Lanfear et al., 2014). Populations living in fragmented habitats, subjected to continuous and severe bottlenecks, are expected to undergo dramatic population size reduction and drift, which can significantly impact the accumulation of genetic polymorphisms in genes affecting important life history traits (Nonaka et al., 2019). The two main evolutionary theories of aging explain aging as the consequence of two fundamentally different processes. The mutation accumulation theory of aging (MA) attributes the evolution of aging to germline-encoded genetic variants accumulating in populations due to the age-dependent weakening of purifying selection, which becomes less efficient to remove from the gene pool gene variants that negatively impact fitness in late life (Charlesworth, 2000). The antagonistic pleiotropy (AP) theory of aging, instead, states that positive selection could favor gene variants that, while overall beneficial for individual fitness, may have detrimental effects in late life (Charlesworth, 2000; Williams, 1957). Although the two theories are not mutually exclusive and can both in principle explain the evolution of aging-related variants across species, their genetic traces in the genome should be distinguishable. In fact, while aging-determining gene variants occurring due to mutation accumulation evolve as nearly-neutral variants, aging-determining gene variants emerging via antagonistic pleiotropism evolve as positive selected variants.

Among vertebrates, killifish represent a unique system, as they repeatedly and independently colonized highly fragmented habitats, characterized by cycles of rainfalls and drought (Furness, 2016). While on the one hand intermittent precipitation and periodic drought pose strong selective pressures leading to the evolution of embryonic diapause, an adaptation that enables killifish to survive in absence of water (Cellerino et al., 2016; Hu and Brunet, 2018), on the other hand they cause habitat and population fragmentation, promoting inbreeding and genetic drift. The co-occurrence of strong selective pressure for early-life on the one hand and population size decline leading to genetic drift on the other hand characterizes life history evolution in African annual killifishes (Cui et al., 2019).

The turquoise killifish (Nothobranchius furzeri) is the shortest-lived vertebrate with a thoroughly documented post-embryonic life, which, in the shortest-lived strains, amounts to four months (Cellerino et al., 2016; Hu and Brunet, 2018; Kim et al., 2016; Blazek, 2017). Turquoise killifish has recently emerged as a powerful new laboratory model to study experimental biology of aging due to its short lifespan and to its wide range of aging-related changes, which include neoplasias (Di Cicco et al., 2011), decreased regenerative capacity (Wendler et al., 2015), cellular senescence (Ahuja et al., 2019; Valenzano et al., 2006), and loss of microbial diversity (Smith et al., 2017). At the same time, while sharing physiological adaptations that enable embryonic diapause and rapid sexual maturation, different wild turquoise killifish populations display differences in lifespan, both in the wild and in captivity (Terzibasi et al., 2008; Valenzano et al., 2015; Vrtílek et al., 2018), making this species an ideal evolutionary model to study the genetic basis underlying life history trait divergence within species.

Characterization of life history traits in wild-derived laboratory strains of turquoise killifish revealed that while different populations have similar rates of sexual maturation (Blazek, 2017), populations from arid regions exhibit the shortest lifespans, while populations from more semi-arid regions exhibit longer lifespans (Blazek, 2017; Terzibasi et al., 2008). Hence, speed of sexual maturation and adult lifespan appear to be independent in turquoise killifish populations. The evolutionary mechanisms responsible for the lifespan differences among turquoise killifish populations are not yet clearly understood. Mapping genetic loci associated with lifespan differences among turquoise killifish populations showed that adult survival has a complex genetic architecture (Valenzano et al., 2015; Kirschner et al., 2012). Here, combining genome sequencing and population genetics, we investigate to what extent genomic divergence in natural turquoise killifish populations that differ in lifespan is driven by adaptive or neutral evolution, compatible with either the antagonistic pleiotropy (AP) theory of aging or with the mutation accumulation (MA) theory of aging, respectively.

Results

Genome assembly improvement and gene annotation

To identify the genomic mechanism that led to the evolution of differences in lifespan between natural populations of the turquoise killifish (Nothobranchius furzeri), we combined the currently available reference genomes (Valenzano et al., 2015; Reichwald et al., 2015) into an improved reference turquoise killifish genome assembly. Due to the high repeat content, genome assembly from short reads required a highly integrated and multi-platform approach. We ran Allpaths-LG with all the available pair-end sequences, producing a combined assembly with a contig N50 of 7.8 kb, corresponding to a ~ 2 kb improvement from the previous versions. Two newly obtained 10X Genomics linked read libraries were used to correct and link scaffolds, resulting in a scaffold N50 of 1.5 Mb, that is a three-fold improvement from the best previous assembly. With the improved continuity, we assigned 92.2% of assembled bases to the 19 linkage groups using two RAD-tag maps (Valenzano et al., 2015). Gene content assessment using the BUSCO method improved “complete” BUSCOs from 91.43% (Valenzano et al., 2015) and 94.59% (Reichwald et al., 2015) to 95.20%. We mapped Genbank N. furzeri RefSeq RNA to the new assembly to predict gene models. The predicted gene model set is 96.1% for ‘complete’ BUSCOs. The overall size of repeated regions (masked regions) is 1.003 Gb, accounting for 66% of the entire genome, that is 20% higher than a previous estimate (Reichwald et al., 2009).

Population genetics of natural turquoise killifish populations

Natural populations of turquoise killifish occur along an aridity gradient in Zimbabwe and Mozambique and populations from more arid regions are associated with shorter captive lifespan (Blazek, 2017; Terzibasi et al., 2008). A QTL study performed between short-lived and long-lived turquoise killifish populations showed a complex genetic architecture of lifespan (measured as age at death), with several genome-wide loci associated with lifespan differences among long-lived and short-lived populations (Valenzano et al., 2015). To further investigate the evolutionary forces shaping genetic differentiation in the loci associated with lifespan among wild turquoise killifish populations, we performed pooled whole-genome-sequencing (WGS) of killifish collected from four sampling sites within the natural turquoise killifish species distribution, which vary in altitude, annual precipitation and aridity (Figure 1—figure supplement 1, Supplementary file 1A). Population GNP is located within the Gonarezhou National Park at high altitude and in an arid climate (Koeppen-Geiger classification ‘BWh’, Figure 1—figure supplement 1), in a region at the outer edge of the turquoise killifish distribution (Figure 1—figure supplement 1; Dorn et al., 2011; Bartáková et al., 2013; Bartáková et al., 2015), which corresponds to the place of origin of the ‘GRZ’ laboratory strain, which has the shortest lifespan of all laboratory strains of turquoise killifish (Terzibasi et al., 2008; Valenzano et al., 2015). Population NF414 (MZCS 414) is located in an arid area in the center of the Chefu river drainage in Mozambique (‘BWh’, Figure 1—figure supplement 1; Dorn et al., 2011; Bartáková et al., 2013; Bartáková et al., 2015), and population NF303 (MZCS 303) is located in a semi-arid area in transition to more humid climate zones in the center of the Limpopo river drainage system (Koeppen-Geiger classification ‘BSh’, Figure 1—figure supplement 1; Dorn et al., 2011; Bartáková et al., 2013; Bartáková et al., 2015). Altitude among localities ranges from 344 m (GNP) to 68 m (NF303, Figure 1—figure supplement 1a and Supplementary file 1A). The temporary habitat of turquoise killifish populations differs in terms of altitude and aridity, as the ephemeral pools at higher altitude are drained earlier and persist for shorter time, while water bodies in habitats at lower altitude last longer (Terzibasi et al., 2008). Population GNP is therefore named ‘dry’, population NF414 is named ‘intermediate’ and population NF303 ‘wet’ throughout the manuscript. The populations used in this study are from localities that belong to the same drainage system as those used in the previous QTL study and their relative position is included in Figure 1—figure supplement 2.

High genetic differentiation and contrasting population demography in dry and wet populations

We asked whether populations from dry, intermediate and wet areas, corresponding to shorter and progressively longer lifespan, differ in genetic variability. We calculated genome-wide estimates of average pairwise difference (π) and genetic diversity (θWatterson) based on 50kb-non-overlapping sliding windows using PoPoolation (Kofler et al., 2011a). We found that π and θWatterson decrease from wet to dry population (θWatterson GNP: 0.0011, θWatterson NF414: 0.0036, θWatterson NF303: 0.0072; πGNP: 0.0009, πNF414: 0.0031, and πNF303: 0.0054). Hence, dry populations have overall smaller genetic diversity than populations from less dry regions. To infer the genetic distance between the populations, we computed the genome-wide pairwise genetic differentiation between populations using FST (Kofler et al., 2011b). Overall, the genetic differentiation between populations ranged between 0.14 and 0.26 and was the highest between the more geographically distant population GNP (dry) and population NF303 (wet) (Figure 1a).

Figure 1 with 2 supplements see all
Download asset Open asset
Demography and natural occurrence of turquoise killifish populations.

(a) Inferred ancestral effective population size (Ne) (using PSMC’) on y-axis and past generations on x-axis in GNP (red, orange), NF414 (black, grey) and NF303 (blue). Inset: unrooted neighbor joining tree based on pairwise genetic differentiation (FST) values. (b) Geographical locations of sampled natural population of turquoise killifish (Nothobranchius furzeri). The area of the colored circles represents the estimated effective population size (Ne) based on θWatterson. (c) Natural environment of turquoise killifish and schematic of the annual life cycle. Figure 1 was partly made with Biorender.

Next, we inferred the demographic history of the populations using pairwise sequentially Markovian coalescent (PSMC) by resequencing at high-coverage single individuals for each population (Schiffels and Durbin, 2014). The population GNP (dry) experienced a strong population decline starting approximately 150 k generations ago, a result consistent for both the sequenced individuals from the two sampling sites (GNP-G1-3 and GNP-G4, Figure 1a). In contrast to the demographic history in GNP, we found indications for recent population expansions in populations from the center of the Chefu and Limpopo basins clades. Analysis of population NF414 (intermediate) (Figure 1a, NF414-Y and NF414-R) and NF303 (wet) (Figure 1a, blue line) shows population expansion until recent time (~50 k generations ago). To infer the effective population size (Ne) of the populations, we used the published mutational rate of 2.6321e−9 per base pair per generation for Nothobranchius, computed via dated phylogeny and θWatterson (Cui et al., 2019). In line with the decrease in genetic diversity from wet to dry population, we found a decrease in Ne estimates (107221.8, 338849.48 and 683693.25 for GNP, NF414 and NF303, respectively; Figure 1b). Hence, our findings show that dry populations from the outer edge of the species distribution show lower genetic diversity and smaller effective population size compared to population from intermediate and more wet regions.

Genetic differentiation among turquoise killifish populations

To test whether regions underlying longevity QTL in turquoise killifish (Valenzano et al., 2015; Kirschner et al., 2012) display a genetic signature for positive or purifying selection in these wild populations, we took advantage of the improved turquoise killifish genome assembly and the newly sequenced wild turquoise killifish populations (Figure 2). The strongest QTL for lifespan differences among long-lived and short-lived populations mapped on the sex chromosome (Valenzano et al., 2015; Kirschner et al., 2012), in proximity to the sex determining locus (Valenzano et al., 2015).

Figure 2
Download asset Open asset
Genomic regions of high and low genetic divergence between pairs of turquoise killifish populations.

Left) Genomic regions with high or low genetic differentiation between turquoise killifish populations identified with an FST outlier approach. Z-transformed FST values of all pairwise comparisons in solid lines, with ‘NF303vsNF414’ in yellow, ‘NF303vsGNP’ in blue, and ‘NF414vsGNP’ in green. The significance thresholds of upper and lower 5‰ are shown as dotted lines with same color coding. Center) Circos plot of Z-transformed FST values between all pairwise comparisons with ‘NF303vsNF414’ in the inner circle (yellow), ‘NF414vsGNP’ in the middle circle (green), and ‘NF303vsGNP’ in the outer circle (blue). Right) Pairwise genetic differentiation based on FST in the four main clusters associated with lifespan (QTL from Valenzano et al., 2015).

To identify a genomic signature of strong selection, we performed an outlier approach based on the pairwise genetic differentiation index (FST). To find highly differentiated regions that may underlie positive selection in natural turquoise killifish populations, we scanned for regions with elevated genetic differentiation between pairs of populations, that is exceeding the 0.995 quantile of Z-transformed non-overlapping 50 kb sliding windows of FST. To find regions under purifying selection, we scanned for regions with lowered genetic differentiation among populations, that is below the 0.005 quantile of Z-transformed non-overlapping 50 kb sliding windows of FST (Supplementary file 1G).

The outlier approach did not reveal clear signatures of positive or purifying selection based on genetic differentiation in the four main chromosomal clusters associated with lifespan in experimental strains of turquoise killifish (Figure 2).

We then analyzed genomic regions carrying signatures of positive and purifying selection in the natural turquoise killifish populations irrespective of the QTL regions (Figure 2). The FST outlier approach led to the identification of several potential regions under strong selection between populations, in particular between the intermediate and wet populations (Supplementary file 1D) and only two between the dry and wet populations (Supplementary file 1E). Genes significantly different and within regions of larger genetic differentiation based on Z-transformed non-overlapping sliding windows of FST were located on chromosomes 6 and 10. The region on chromosome six includes the gene slc8a1, which contains mutations with significant difference in allele frequencies between the wet and intermediate population (Fisher’s exact test implemented in PoPoolation; adjusted p value < 0.001). The region on chromosome 10 contains four genes: XM_015941868, XM_015941869, lss and hibch. All genes under the major FST peak on chromosome 10 showed significant difference in allele frequencies between the intermediate and wet population (Fisher’s exact test; adjusted p value < 0.001) and additionally, hibch had significantly different allele frequencies between the dry and wet population (Fisher’s exact test; adjusted p value < 0.001).

Age-specific changes in genes with sequence divergence between populations

Genes under FST peaks between populations that differ in lifespan, are not necessarily causally involved in lifespan differences between populations, as sequence differences could segregate in populations due to population structure and drift. However, to test whether the genes located in genomic regions that are significantly divergent between populations could be functionally involved in age-related phenotypes, we investigated whether gene expression in these genes varied as a function of age. Analyzing available turquoise killifish longitudinal RNA-Seq datasets generated in liver, brain and skin (Baumgart et al., 2017), we found that hibch, lss and slc8a1 are differentially expressed between adult and old killifish (Supplementary file 1J, adjusted p value < 0.01). hibch, lss and slc8a1 are involved in amino acid metabolism (Ferdinandusse et al., 2013), biosynthesis of cholesterol (Huff and Telford, 2005), and proton-mediated accelerated aging (Osanai et al., 2018), respectively. Gene XM_015956265 (ZBTB14) is the only gene that is an FST outlier and that is differentially expressed in adult vs. old individuals between at least two populations in all tissues (liver, brain and skin). XM_015956265 encodes a transcriptional modulator with ubiquitous functions, ranging from activation of dopamine transporter to repression of myc, fmr1 and thymidine kinase promoters (Orlov et al., 2007). However, although genomic regions that have sequence divergence between turquoise killifish populations contain genes that are differentially expressed during aging in different tissues, whether any of these genes are causally involved in modulating aging-related changes between turquoise killifish wild populations still remains to be assessed. We could not find enrichment of significant differentially expressed genes within the FST outlier regions (Fisher’s exact test p value > 0.05).

Genomic regions of low genetic differentiation among populations

Based on the outlier approach, we found two genomic regions with low genetic differentiation between all pairs of populations, suggesting strong purifying selection. The first region is located on the sex chromosome and contains the putative sex determining gene gdf6 (Reichwald et al., 2015), which is hence conserved among these populations. This same region also contains sybu, a maternal-effect gene associated with the establishment of embryo polarity (Nojima et al., 2010). The second region under low genetic differentiation is located on chromosome nine and harbors the genes XM_015965812 (abi2-like), cnot11 and lcp1, which are involved in phagocytosis (Ulvila et al., 2011), mRNA degradation (Mauxion et al., 2013) and cell motility (Kell et al., 2018), respectively. Signatures of low and high genetic differentiation between populations can be the result of purifying or positive selection. However, balancing selection, a mechanism that could maintain polymorphism above the expected genetic diversity, could also in part result in genetic differentiation (Brandt et al., 2018). To account for balancing selection, we compared the pairwise genetic diversity (π) among populations and we could not find signatures of elevated genetic diversity within the investigated regions under strong selection.

Hence, we could not find a clear evidence of positive or purifying selection in correspondence with the survival QTL previously identified, suggesting that genomic regions associated with natural lifespan differences may have not evolved due to positive selection or have being maintained under purifying selection. However, we cannot exclude that we could not detect positive selection at the QTL regions due to statistical power or that the populations used in this study and those used for the QTL analysis had a different genetic architecture of lifespan.

Evolutionary origin of the sex chromosome

Since we found reduced genetic differentiation among populations in the chromosomal region containing the putative sex-determining gene in the sex chromosome, we used synteny analysis and the new genome assembly to investigate the genomic events that led to evolution of this chromosomal region (Figure 3). We found that the structure of the turquoise killifish sex chromosome is compatible with a chromosomal translocation within an ancestral chromosome and a fusion event between two chromosomes. The translocation event within an ancestral chromosome corresponding to medaka´s chromosome 16 and platyfish´s linkage group three led to a repositioning of a chromosomal region containing the putative sex-determining gene gdf6 (Figure 3b). The fusion of the translocated chromosome with a chromosome corresponding to medaka´s chromosome eight and platyfish´s linkage group 16, possibly led to the origin of turquoise killifish's sex chromosome. We could hence reconstruct a model for the origin of the turquoise killifish sex chromosome (Figure 3c), which parsimoniously places a translocation event before a fusion event. The occurrence of two major chromosomal rearrangements could have then contributed to suppressing recombination around the sex-determining region (Valenzano et al., 2015; Valenzano et al., 2009).

Figure 3
Download asset Open asset
Synteny and sex chromosome evolution in turquoise killifish.

(a) Synteny circos plots based on 1-to-1 orthologous gene location between the new turquoise killifish assembly (black chromosomes) and platyfish (Xiphophorus maculatus, colored chromosomes, left circos plot) and between the new turquoise killifish assembly (black chromosomes) and medaka (Oryzias latipes, colored chromosomes, right circos plot). Orthologous genes in concordant order are visualized as one syntenic block. Synteny regions are connected via color-coded ribbons, based on their chromosomal location in platyfish or medaka. If the direction of the syntenic sequence is inverted compared to the compared species, the ribbon is twisted. Outer data plot shows –log(q-value) of survival quantitative trait loci (QTL, ordinate value between 0 and 3.5, every value above 3.5 is visualized at 3.5 [Valenzano et al., 2015]) and the inner data plot shows –log(q-value) of the sex QTL (ordinate value between 0 and 3.5, every value above 3.5 is visualized at 3.5). Boxes between the two circos plots show genes within the peak regions of the four highest –log(q-value) of survival QTL on independent chromosomes (red box) and the highest association to sex (black box). (b) High resolution synteny map between the sex-chromosome of the turquoise killifish (Chr3) with platyfish chromosome 16 and 3 in the upper plot, and between the turquoise killifish and medaka chromosome 8 and 16 (lower plot). The middle plot shows the QTLs for survival and sex along the turquoise killifish sex chromosome. (c) Model of sex chromosome evolution in the turquoise killifish. A translocation event within one ancestral autosome led to the emergence of a chromosomal region harboring a new sex-determining-gene (SDG). The fusion of a second autosome led to the formation of the current structure of the turquoise killifish sex chromosome.

Relaxed selection in turquoise killifish populations

Since we could not identify specific signatures of genetic differentiation in the genomic regions associated with longevity from previous QTL mapping, we asked whether other evolutionary forces than directional selection may underlie differences in survival among wild turquoise killifish populations. The difference in the recent and past demography between populations (Figure 1) led us to ask whether demography could have led to evolutionary changes on genome-wide scale between natural populations. For each population, we calculated the fraction of substitutions driven to fixation by positive selection since divergence from the outgroup species Nothobranchius orthonotus (NOR) using the asymptotic McDonald-Kreitman α (Messer and Petrov, 2013). The original McDonald-Kreitman α (which ranges from – ∞ to 1) was designed to calculated the rate of adaptation by comparing the polymorphisms (within species) and divergence (between species) at neutral and functional sites (McDonald and Kreitman, 1991). While McDonald-Kreitman α = 0 indicates neutrality, larger and positive values of α mean that a given population has an elevated proportion of genetic variants driven by natural selection, while negative values of α can be an indication of deleterious variants. The asymptotic McDonald-Kreitman α accounts for a range of derived allele frequencies, enabling to identify slightly deleterious mutations as those segregating at lower derived allele frequencies (Messer and Petrov, 2013). Variants at low derived allele frequency are either neutral (if they were beneficial they would have higher frequency) or are slightly deleterious. Hence, negative McDonald-Kreitman α values at low derived allele frequency bins likely reflect slightly deleterious gene variants. Additionally, negative α values at intermediate to higher frequency bins may indicate drift of deleterious variants. We set out to adopt this method to assess the genetic variants for each population, compared to two different outgroup species.

Using Nothobranchius orthonotus (NOR) as an outgroup, we inferred the fraction of positive selection by pooling all coding sites (Figure 4a). SNPs were called with the program SNAPE (Raineri et al., 2012), which specifically deals with pooled sequencing. We only included SNPs with a derived frequency between 0.05–0.95 and performed stringent filtering. The asymptotic McDonald-Kreitman α ranged from −0.21 to −0.01 in comparison to the very closely related sister species N. orthonotus, confirming limited genome-wide positive selection since divergence from N. orthonotus (Figure 4a). The population GNP, located in an arid region at higher altitude and associated with the shortest recorded lifespan, shows the lowest asymptotic McDonald-Kreitman α, as well as lower McDonald-Kreitman α values throughout all derived frequency bins, potentially suggesting a higher load of slightly deleterious mutations segregating in this population (Figure 4a). Using as an outgroup species another annual killifish species, Nothobranchius rachovii (NRC), we confirmed the lowest asymptotic McDonald-Kreitman α value in the dry population GNP (Figure 4b). Additionally, using Nothobranchius rachovii (NRC) as outgroup species, the asymptotic McDonald-Kreitman α ranged from −0.06 to 0.23 among populations, indicating that more alleles were driven to fixation by positive selection in the ancestral lineage leading to Nothobranchius furzeri and Nothobranchius orthonotus. In particular, the wet population NF303 had the highest asymptotic McDonald-Kreitman α value (Figure 4b). Using both N. orthonotus and N. rachovii as outgroups, we found that the dry GNP population had the lowest McDonald-Kreitman α values at the low derived frequency bins, potentially consistent with a genome-wide accumulation of slightly deleterious mutations in these isolated populations.

Figure 4 with 2 supplements see all
Download asset Open asset
Genome-wide signatures of natural and relaxed selection in turquoise killifish populations.

Asymptotic McDonald-Kreitman alpha (MK α) analysis based on derived frequency bins using as outgroups (a) Nothobranchius orthonotus and (b) Nothobranchius rachovii. Population GNP is shown in red, NF414 in black, and NF303 in blue. (c) Proportion of non-synonymous SNPs binned in allele frequencies of non-reference (alternative) alleles for GNP (red), NF414 (black) and NF303 (blue). (d) Negative distribution of fitness effects of populations GNP (red), NF414 (black) and NF303 (blue) with cumulative proportion of deleterious SNPs on y-axis and the compound measure of 4Nes on x-axis. (e) Proportion of different effect types of SNPs in coding sequences of all populations. The effect on amino acid sequence for each genetic variant is represented by colors (legend). Significance is based on ratio between synonymous effects to non-synonymous effects (significance based on Chi-square test).

Estimating the distribution of fitness effect across populations

To directly estimate the fitness effect of gene variants associated with each population, we analyzed population-specific genetic polymorphisms to assign mutations as beneficial, neutral or detrimental, and determine the distribution of fitness effect (DFE) (Tataru et al., 2017) of new mutations. Consistently with the overall lower McDonald-Kreitman α values throughout all derived frequency bins, we found more new mutations assigned as the slightly deleterious category in the dry GNP population, compared to the other two populations (indicated by the higher number of deleterious SNPs in proximity to 4NeS ~ 0 in the GNP population, Figure 4d, Supplementary file 1H). To independently validate our findings, we ran a simulation using SLiM3, which recapitulated the population divergence from an ancestral population, followed with diverging population size as inferred from the PSMC’ analysis (Figure 4—figure supplement 1). Analyzing the distribution of fitness effect in these simulated populations, we confirmed that populations with smaller effective population size have a higher proportion of new slightly deleterious variants, compared to larger populations, which have relatively more newly arising gene variants that are highly deleterious, indicating that purifying selection in the large population is more efficient in removing mutations with deleterious effects. To further infer the effect of the putative deleterious mutations on protein function, we used the new turquoise killifish genome assembly as a reference and adopted an approach that, by analyzing sequence polymorphism among populations, predicts functional consequences at the protein level (Cingolani et al., 2012). We found that the proportion of mutations causing a change in protein function is significantly larger in the GNP population compared to populations NF414 and NF303 (Chi-square test: PGNP-NF303 <1.87e-119, PGNP-NF414 <4.96e-57, PNF303-NF414< 3.51e-35, Figure 4e). Additionally, the mutations with predicted deleterious effects on protein function reached also higher frequencies in the dry population GNP (Figure 4c).

Distribution of mutations at conserved sites

To further investigate the impact of mutations on protein function, we calculated the Consurf (Pupko et al., 2002; Mayrose et al., 2004; Glaser et al., 2003; Ashkenazy et al., 2016) score, which determines the evolutionary constraint on an amino acid, based on sequence conservation. Mutations at amino acid positions with high Consurf score (i.e. otherwise highly conserved) are considered to be more deleterious. We found that the dry population GNP had a significantly higher mean Consurf score for mutations at non-synonymous sites in frequency bins from 5–20% up to 40–60%, compared to populations NF414 (intermediate) and NF303 (wet) (Figure 4—figure supplement 2). The mutations in the dry GNP population had significantly higher Consurf scores than the other populations using both outgroup species N. orthonotus and N. rachovii (Figure S3). Upon exclusion of potential mutations at neighboring sites (CMD: codons with multiple differences), CpG hypermutation and genes containing mutations with highly detrimental effect on protein function based on SnpEFF analysis, the dry population GNP had higher mean Consurf score at the low frequency bin (Figure 4—figure supplement 2, Supplementary file 1K-L). To note, we also found a significantly higher average Consurf score at synonymous sites in GNP at low derived frequencies (Figure 4—figure supplement 1, Supplementary file 1K-L), possibly suggestive of an overall higher mutational rate in GNP.

Relaxation of selection in age-related disease pathways

To further test whether populations from dry environments accumulated a higher load of deleterious gene variants, we computed the gene-wise direction of selection (DoS) (Stoletzki and Eyre-Walker, 2011) index, which measures the strength of selection based on the count of mutations in non-synonymous and synonymous sites. Indeed, we found support to the hypothesis that the dry, short-lived population GNP has significantly more slightly deleterious mutations segregating in the population, compared to the populations NF414 and NF303 (Figure 5a, Median NOR: GNP: −0.17, NF414: −0.02, NF303: −0.01; Median NRC: GNP: −0.14, NF414: 0.00, NF303: 0.00; Wilcoxon rank sum test: NOR: PGNPNF303 <2.21e-105, PGNP-NF414 <1.19e-76, PNF303-NF414< 1.39e-06; NRC: PGNP-NF303 <4.61e-179, PGNP-NF414 <1.42e-100, PNF303-NF414< 5.96e-22), indicating that purifying selection is relaxed in GNP. We calculated DoS in all populations using independently as outgroup species N. orthonotus and N. rachovii (Figure 5a).

Figure 5
Download asset Open asset
Pathway enrichment in genes under adaptive and neutral evolution in turquoise killifish populations.

(a) Distribution of direction of selection (DoS) represented with median of distribution for population GNP (red), NF414 (grey) and NF303 (blue). Left panel shows DoS distribution computed using Nothobranchius orthonotus as outgroup and right panel shows DoS distribution computed using Nothobranchius rachovii as outgroup. Significance based on Wilcoxon-Rank-Sum test. (b) Pathway over-representation analysis of genes below the 2.5% level of gene-wise DoS values are shown with red background and above the 97.5% level of gene-wise DoS values are shown with green background. Only pathway terms with significance level of FDR corrected q-value <0.05 are shown (in -log(q-value)). Terms enriched in population GNP have red dots, enriched in population NF414 have black dots, and enriched in population NF303 have blue dots, respectively.

To assess whether specific biological pathways were significantly more impacted by the accumulation of slightly deleterious mutations, we performed pathway overrepresentation analysis. We found a significant overrepresentation in the lower 2.5th DoS quantile (i.e. genes under relaxation of selection) in the GNP population for pathways associated with age-related diseases, including gastric cancer, breast cancer, neurodegenerative disease, mTOR signaling and WNT signaling (q-value <0.05, Figure 5b, Supplementary file 1I). Overall, relaxed selection in the dry GNP population affected accumulation of deleterious mutations in age-related and in the WNT pathway. Analyzing the pathways affected by genes within the upper 2.5th DoS values – corresponding to genes undergoing adaptive evolution – we found a significant enrichment for mitochondrial pathways – potentially compensatory (Cui et al., 2019) – in population NF303 (Figure 5b, Supplementary file 1I). Overall, our results show that differences in effective population size among wild turquoise killifish are associated with an extensive relaxation of purifying selection, significantly affecting genes involved in age-related diseases, and which could have cumulatively contributed to reducing individual survival.

Discussion

The turquoise killifish (Nothobranchius furzeri) is the shortest-lived known vertebrate and while its natural populations show similar timing for sexual maturation, exhibit differences in lifespan along a cline of altitude and aridity in south-eastern Africa (Blazek, 2017; Terzibasi et al., 2008). Here we generate an improved genome assembly (NFZ v2.0) in turquoise killifish (Nothobranchius furzeri) and study the evolutionary forces shaping genome evolution among natural populations.

Using the new turquoise killifish genome assembly and synteny analysis with medaka and platyfish, we reconstructed the origin of the turquoise killifish sex chromosome, which appears to have evolved through two independent chromosomal events, that is a translocation and a fusion event.

Using the new genome assembly and pooled sequencing of natural turquoise killifish populations, we found that genetic differentiation among populations of the short-lived turquoise killifish is consistent with differences in demographic constraints. While we found that strong purifying selection maintains low genetic diversity among populations at genomic regions underlying key species-specific traits, such as in proximity to the sex-determining region, demography and genetic drift largely shape genome evolution, leading to relaxation of selection and the accumulation of deleterious mutations. We showed that isolated populations from an arid region, dwelling at higher altitude and characterized by shorter lifespan, experienced extensive population bottlenecking and a sharp decline in effective population size. Populations from dryer regions at higher altitudes experience genetic isolation and possibly steady decline in population size due to limited incoming gene flow and possibly more severe bottlenecks due to recent founder effect. However, populations from more wet regions likely undergo extensive gene flow, maintaining larger population size. We found that relaxation of selection in more drifted populations significantly affected the accumulation of deleterious gene variants in pathways associated with neurodegenerative diseases and WNT-signaling (Figure 5). While simple traits, such as male tail color and sex have a simple genetic architecture among turquoise killifish populations (Valenzano et al., 2015; Valenzano et al., 2009), we find that the complex genetic architecture of lifespan differences among killifish populations (Valenzano et al., 2015) is entirely compatible with genome-wide relaxation of selection. Additionally, the absence of genomic signature of positive selection in genomic regions underlying survival QTL in killifish suggest that, rather than directional selection, the neutral accumulation of deleterious mutations may be the evolutionary mechanism underlying survival differences among turquoise killifish populations, in line with the mutation accumulation theory of aging. The antagonistic pleiotropy theory of aging states that positive selection could lead to the fixation of gene variants that, while overall beneficial for fitness, could reduce survival and reproductive capacity in late life (Williams, 1957). However, the lack of genomic signatures of positive selection at the genomic regions underlying survival QTL in turquoise killifish rather suggests that the accumulation of deleterious mutations due to neutral drift may have played a key role in shaping genome and phenotype differences among natural turquoise killifish populations. One of the deductions of the antagonistic pleiotropy theory is that a reduction in speed of maturation should be associated with increased lifespan (Williams, 1957). However, different wild populations of turquoise killifish have similar time to sexual maturation and yet different lifespan (Blazek, 2017). Hence, the uncoupling of age of sexual maturation from adult lifespan in different turquoise killifish populations is more compatible with the mutation accumulation theory of aging. However, although we did not find evidence for it, our results do not exclude a priori the possibility that genes under strong selection may in part contribute to lifespan differences among different turquoise killifish populations, hence acting compatibly with the antagonistic pleiotropy theory of aging. However, historical fluctuations in the size of natural turquoise killifish populations, especially in isolated populations living in more arid and elevated habitats, weakened the strength of natural selection, ultimately contributing to increased load of deleterious gene variants, preferentially in genes associated with aging-related diseases and in the WNT pathway. We hypothesize that small effective population size leads to the accumulation of aging-causing mutations that together contribute to the genetic architecture of lifespan. Overall, our findings highlight the role of demographic constraints in shaping life history within species.

Materials and methods

Merging and improvement of the turquoise killifish genome assembly (BioProject ID: PRJNA599375)

10x genomics read clouds

Request a detailed protocol

A single GRZ male individual was sacrificed with MS222 (Sigma-Aldrich, Steinheim, Germany). Blood was drawn from the heart and high molecular weight DNA was isolated with Qiagen MagAttract kit following manufacturer’s instructions. Gemcode v2 DNA library generation was performed by Novogene (Beijing, China). Briefly, a proportion of the sample was run on a pulse field agarose gel to confirm high molecularity >100 kb. Based on a genome size estimate of 1.54 Gb (half of human genome), 0.6 ng of DNA was used to construct 2 Gemcode libraries, sequenced on two HiSeq X lanes to obtain a raw coverage of approximately 60X each. The reported input molecular length by SuperNova (Weisenfeld et al., 2017) was 118 kb for library 1 and 60.73 kb for library 2. Both libraries were used to correct and scaffold the Allpath-LG assembly (see below), and library one was also de novo assembled with the SuperNova assembler v.2 with default parameters. The SuperNova assembly totaled 802.6 Mb, with a contig N50 of 19.65 kb, scaffolded into 6.78 thousand scaffolds with an N50 of 3.83 Mb. Despite high continuity, however, the BUSCO (Simão et al., 2015) metrics are much lower than the Allpath-LG assemblies.

Nanopore long reads

Request a detailed protocol

DNA was extracted from a single GRZ male individual’s muscle tissue by grinding in liquid nitrogen followed by phenol-chloroform extraction (Sigma). The rapid sequencing kit (SQK-RAD004) and the ligation kit (SQK-LSK108) were sued to prepare six libraries and were sequenced on 6 MinION flow cells (R9.4.1). These runs yielded a total of 3.3 Gb of sequences after trimming and correction by HALC (Bao and Lan, 2017). For correction, Allpath-LG contigs (see below) and short reads from the 10X genomic run were used.

Allpath-LG assembly

Request a detailed protocol

Two independent short read datasets were previous collected for the GRZ strain of Nothobranchius furzeri. Allpath-LG (Gnerre et al., 2011) was used on the pooled datasets. Together, 4 Illumina short read pair-end libraries with a fragment size distribution from 158 bp to 179 bp were used to construct the contigs (sequence coverage 191.9X, physical coverage 153.5X), and 22 pair-end and mate pair libraries distributed at 92 bp, 135 bp, 141 bp, 176 bp, 267 bp, 2 kb, 3 kb 5 kb and 10 kb were used for the scaffolding step (sequence coverage 135.7X, physical coverage 453.8X). The published BAC library ends (Reichwald et al., 2015) with an insert size of 112 kb were also included in the ALLPaths-LG run (physical coverage 0.6X). The resulting assembly has a total contig length of 823,583,106 bp distributed in 151,307 contigs > 1 kb, with an N50 of 7.8 kb. The total scaffold length is 943,793,727 bp distributed in 7830 scaffolds with an N50 of 421 kb (with gaps). The resulting assembly was further scaffolded by ARCS v1.0 (Coombe et al., 2018) + LINKS v1.8.5 (Warren et al., 2015) with the following parameters: arcs -e 50000 c 3 r 0.05 s 98 and LINKS -m -d 4000 k 20 -e 0.1 l 3 -a 0.3 t 2 -o 0 -z 500 r -p 0.001 -x 0. This increased the scaffold N50 to 1.527 Mb. Next, scaffolds were assigned to the RAD-tag linkage map (Valenzano et al., 2015) collected from a previous study with Allmaps (Tang et al., 2015), using equal weight for the two independent mapping crosses. This procedure assigned 90.6% of the assembled bases in 1131 scaffolds to 19 linkage groups, in which 76.6% can be oriented. Misassemblies were corrected with the 10X genomic read cloud. Read clouds were mapped to the preliminary assembly with longranger v2.1.6 using default parameters, and a custom script was used to scan for sudden drops in barcode shares along the assembled linkage groups. The scaffolds were broken at the nearest gap of the drop in 10x barcodes. The same ARCS + LINKS pipeline was again run on the broken scaffolds, increasing the scaffold N50 to 1.823 Mb. Next, BESST_RNA (https://github.com/ksahlin/BESST_RNA) was used to further scaffold the assembly with RNASeq libraries, Allmaps was again used to assign the fixed scaffolds back to linkage groups, increasing the assignable bases to 92.2% (879Mb) with 80.3% (765Mb) with determined orientation. The assembly was again broken with longranger and reassigned to LG with Allmaps, and the scaffolds were further partitioned to linkage groups due to linkage of some left-over scaffolds with an assigned scaffold. Each partitioned scaffold groups were subjected to the ARCS + LINKS pipeline again, to constraint the previously unassigned scaffolds onto the same linkage group. Allmaps was run again on the improved scaffolds, resulting in 94.5% (903.4Mb) of bases assigned and 89.1% (852Mb) of bases oriented. Longranger was run again, visually checked and compared with the RADtag markers. Eleven mis-oriented positions were identified and corrected. Gaps were further patched by GMCloser (Kosugi et al., 2015) with ~2X of nanopore long reads corrected by HALC using BGI500 short PE reads with the following parameters: gmcloser --blast --long_read --lr_cov 2 l 100 -i 466 -d 13 min_subcon 1 min_gap_size 10 --iterate 2 --mq 1 c. The corrected long reads not mapped by GMCloser were assembled by CANU (Koren et al., 2017) into 7.9 Mb of sequences, which are likely unassigned repeats.

Meta assembly

Request a detailed protocol

Five assemblies were integrated by MetAssembler (Wences and Schatz, 2015) in the following order (ranked by BUSCO scores) using a 20 kb mate pair library: 1) The improved Allpaths-LG assembly assigned to linkage groups produced in this study, 2) A previously published assembly with Allpaths-LG and optical map (Reichwald et al., 2015) 3) A previously published assembly using SGA (Valenzano et al., 2015), 4) The SuperNova assembly with only 10x Genomic reads and 5) Unassigned nanopore contigs from CANU. The final assembly NFZ v2.0 has 911.5 Mb of scaffolds assigned to linkage groups. Unassigned scaffolds summed up to 142.2 Mb, yielding a total assembly length of 1053.7 Mb, approximately 2/3 of the total genome size of 1.53 Gb. The final assembly has 95.2% complete and 2.24% missing BUSCOs.

Mapping of NCBI Genbank gene annotations

Request a detailed protocol

RefSeq mRNAs for the GRZ strain (PRJNA314891, PRJEB5837) were downloaded from GenBank (Sayers et al., 2019), and aligned to the assembly with Exonerate (Slater and Birney, 2005). The RefSeq mRNAs have a BUSCO score of 98.0% complete, 0.9% missing. The mapped gene models resulted in a BUSCO score of 96.1% complete, 2.1% missing.

Pseudogenome assembly generation

Request a detailed protocol

The pseudogenomes for Nothobranchius orthonotus and Nothobranchius rachovii were generated from sequencing data and the same method used in Cui et al., 2019. Briefly, the sequencing data were mapped to the NFZ v2.0 reference genome by BWA-mem v0.7.12 in PE mode (Li, 2013; Li and Durbin, 2010). PCR duplicates were marked with MarkDuplicates tool in the Picard (version 1.119, http://broadinstitute.github.io/picard/) package. Reads were realigned around INDELs with the IndelRealigner tool in GATK v3.4.46 (McKenna et al., 2010). Variants were called with SAMTOOLS v1.2 (Li et al., 2009) mpileup command, requiring a minimal mapping quality of 20 and a minimal base quality of 25. A pseudogenome assembly was generated by substituting reference bases with the alternative base in the reads. Uncovered regions, INDELs and sites with >2 alleles were masked as unknown ‘N’. The allele with more supporting reads was chosen at biallelic sites.

Mapping of longevity and sex quantitative trait loci

Request a detailed protocol

The quantitative trait loci (QTL) markers published in Valenzano et al., 2015 were directly provided by Dario Riccardo Valenzano. In order to map the markers associated with longevity and sex, a reference database was created using BLAST (Altschul et al., 1990). The nucleotide database was created with the new reference genome of N. furzeri (NFZ v2.0). Subsequently, the QTL marker sequences were mapped to the database. Only markers with full support for the total length of 95 bp were considered as QTL markers.

Synteny analysis

Request a detailed protocol

Synteny analysis was performed using orthologous information from Cui et al., 2019 determined by the UPhO pipeline (Ballesteros and Hormiga, 2016). For this, the 1-to-1 orthologous gene positions of the new turquoise killifish reference genome (NFZ v2.0) were compared to two closely related teleost species, Xiphophorus maculatus and Oryzias latipes. Result were visualized using Circos (Krzywinski et al., 2009) for the genome-wide comparison and the genoPlotR package (Guy et al., 2010) in R for the sex chromosome synteny analysis. Synteny plots for orthologous chromosomes of Xiphophorus maculatus and Oryzias latipes were generated with Synteny DB (http://syntenydb.uoregon.edu) (Catchen et al., 2009).

Koeppen-Geiger index and bioclimatic variables

Request a detailed protocol

The Koeppen-Geiger classification data was taken from Peel et al., 2007 and the altitude, precipitation per month, and the bioclimatic variables were obtained from the Worldclim database (v2.0 [Fick and Hijmans, 2017]). The monthly evapotranspiration was obtained from Trabucco and Zomer, 2019. Aridity index was calculated based on the sum of monthly precipitation divided by sum of monthly evapotranspiration. Maps in Figure 1—figure supplement 1 were generated with QGIS version 2.18.20 combined with GRASS version 7.4 (Neteler et al., 2012), the Koeppen-Geiger raster file, data from Natural Earth, and the river systems database from Lehner and Verdin, 2006.

DNA isolation and pooled population sequencing

Request a detailed protocol

The ethanol preserved fin tissue was washed with 1X PBS before extraction. Fin tissue was digested with 10 µg/mL Proteinase K (Thermo Fisher) in 10 mM TRIS pH 8; 10 mM EDTA; 0.5 SDS at 50°C overnight. DNA was extracted with phenol-chloroform-isoamylalcohol (Sigma) followed by a washing step with chloroform (Sigma). Next, DNA was precipitated by adding 2.5 vol of chilled 100% ethanol and 0.26 vol of 7.5M Ammonium Acetate (Sigma) at −20°C overnight. DNA was collected via centrifugation at 4°C at 12000 rpm for 20 min. After a final washing step with 70% ice-cold ethanol and air drying, DNA was eluted in 30 µl of nuclease-free water. DNA quality was checked on one agarose gels stained with RotiSafe (Roth) and a UV-VIS spectrometer (Nanodrop2000c, Thermo Scientific). DNA concentration was measured with Qubit fluorometer (BR dsDNA Assay Kit, Invitrogen). For each population, the DNA of the individuals were pooled at equimolar contribution (GNP_G1_3, GNP_G4 N = 29; NF414, NF303 N = 30). DNA pools were given to the Cologne Center of Genomic (CCG, Cologne, Germany) for library preparation. The total amount of DNA provided to the sequencing facility was 3.2 μg per pooled population sample. Libraries were sequenced with 150 bp x two paired-ends on the HiSeq4000. Sequencing of pooled samples resulted in a range of 419–517 million paired-end reads for each population (Supplementary file 1B).

Mapping of pooled sequencing reads

Request a detailed protocol

Raw sequencing reads were trimmed using Trimmomatic-0.32 (ILLUMINACLIP:illumina-adaptors.fa:3:7:7:1:true, LEADING:20, TRAILING:20, SLIDINGWINDOW:4:20, MINLEN:50 [Bolger et al., 2014]). Data files were inspected with FastQC (version 0.11.22, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Trimmed reads were subsequently mapped to the reference genome with BWA-MEM v0.7.12 (Li, 2013; Li and Durbin, 2010). The SAM output was converted into BAM format, sorted, and indexed via SAMTOOLS v1.3.1 (Li et al., 2009). Filtering and realignment was conducted with PICARD v1.119 (http://broadinstitute.github.io/picard/) and GATK (McKenna et al., 2010). Briefly, the reads were relabeled, sorted, and indexed with AddOrReplaceReadGroups. Duplicated reads were marked with the PICARD feature MarkDuplicates and reads were realigned with first creating a target list with RealignerTargetCreator, second by IndelRealigner from the GATK suite. Resulting reads were again sorted and indexed with SAMTOOLS. For population genetic bioinformatics analyses the BAM files of the pooled populations were converted into the required MPILEUP format via the SAMTOOLS mpileup command. Low quality reads were excluded by setting a minimum mapping quality of 20 and a minimum base quality of 20. Further, possible insertion and deletions (INDELs) were identified with identifygenomic-indel-regions.pl script from the PoPoolation package (Kofler et al., 2011a) and were subsequently removed via the filter-pileup-by-gtf.pl script (Kofler et al., 2011a). Coding sequence positions that were identified to be putative ambiguous were removed by providing the filter-pileup-by-gtf.pl script a custom modified GTF file with the corresponding coordinates. After adapter and quality filtering, mapping to the newly assembled reference genome resulted in mean genome coverage of 35x, 39x, and 47x for the population NF303, NF414, and GNP, respectively (Supplementary file 1B).

Merging sequencing reads of populations from the Gonarezhou National Park

Request a detailed protocol

Population GNP consists of two sampling sites (GNP-G1_3, GNP-G4) with very low genetic differentiation (Figure 1—figure supplement 1c, Supplementary file 1C). Sequencing reads of the two populations from the Gonarezhou National Park (GNP) were combined used the SAMTOOLS ‘merge’ command. The populations GNP-G1-3 and GNP-G4 were merged together and this population was subsequently denoted as GNP.

Estimating genetic diversity

Request a detailed protocol

Genetic diversity in the populations was estimated by calculating the nucleotide diversity π (Nei and Li, 1979) and Wattersons’s estimator θ (Watterson, 1975). Calculation of π and θ was done with a sliding window approach by using the Variance-sliding.pl script from the PoPoolation program (Kofler et al., 2011a). Non-overlapping windows with a length of 50 kb with a minimum count of two per SNP, minimum quality of 20 and the population specific haploid pool size were used (GNP = 116; NF414 = 60; NF303 = 60). Low covered regions that fall below half the mean coverage of each population were excluded (GNP = 23; NF414 = 19; NF303 = 18), as well as regions that exceed a two times higher coverage than the mean coverage (GNP = 94; NF414 = 77; NF303 = 70). The upper threshold is set to avoid regions with possible wrong assemblies. Mean coverage was estimated on filtered MPILEUP files. Each window had to be at least covered to 30% to be included in the estimation.

Estimation of effective population size

Request a detailed protocol

Wattersons’s estimator of θ (Watterson, 1975) is referred to as the population mutation rate. The estimate is a compound parameter that is calculated as the product of the effective population size (Ne), the ploidy (2 p, with p is ploidy) and the mutational rate μ (θ = 2pNeμ). Therefore, Ne can be obtained when θ, the ploidy and the mutational rate μ are known. The turquoise killifish is a diploid organism with a mutational rate of 2.6321e−nine per base pair per generation (assuming one generation per year in killifish Cui et al., 2019 and θ estimates were obtained with PoPoolation (see previous Section) (Kofler et al., 2011a).

Estimating population differentiation index FST

Request a detailed protocol

The filtered and realigned BAM files of each population were merged into a single pileup file with SAMTOOLS mpileup, with a minimum mapping quality and a minimum base quality of 20. The pileup was synchronized using the mpileup2sync.jar script from the PoPoolation2 program (Kofler et al., 2011b). Insertions and deletions were identified and removed with the identify-indel-regions.pl and filter-sync-by-gtf.pl scripts of PoPoolation2 (Kofler et al., 2011b). Again, coding sequence positions that were identified to be putative ambiguous were removed by providing the filter-pileup-by-gtf.pl script a custom modified GTF file with the corresponding coordinates. Further a synchronized pileup file for genes only were generated by providing a GTF file with genes coordinates to the create-genewise-sync.pl from PoPoolation2 (Kofler et al., 2011b). FST was calculated for each pairwise comparison (GNP vs NF303, GNP vs NF414, NF414 vs NF303) in a genome-wide approach using non-overlapping sliding windows of 50 kb with a minimum count of four per SNP, a minimum coverage of 20, a maximum coverage of 94 for GNP, 77 for NF414, and 70 for NF303 and the corresponding pool size of each population (N = 116; 60; 60). Each sliding window had to be at least covered to 30% to be included in the estimation. The same thresholds, except the minimum covered fraction, with different sliding window sizes were used to calculate the gene-wise FST for the complete gene body (window-size of 2000000, step-size of 2000000) and single SNPs within genes (window-size of 1, step-size of 1). The non-informative positions were excluded from the output. Significance of allele differences per base-pair within the gene-coordinates were calculated with the fisher´s exact test implemented in the fisher-test.pl script of PoPoolation2 (Kofler et al., 2011b). Calculation of unrooted neighbor joining tree based on the genome-wide pairwise FST averages was performed with the ape package in R (Paradis et al., 2004).

Detecting signatures of selection based on FST outliers

Request a detailed protocol

For FST-outlier detection, the pairwise 50kb-window FST-values for each comparison were Z-transformed (ZFST). Next, regions potentially under strong selection were identified by applying an outlier approach. Outliers were identified as non-overlapping windows of 50 kb within the 0.5% of lowest and highest genetic differentiation per comparison. To reduce the number of false-positive results, the outlier threshold was chosen at 0.5% highest and lowest percentile of each pairwise genetic differentiation (Pruisscher et al., 2018; Guo et al., 2016). To find candidate genes within windows of highest differentiation, a total of three selection criteria were used. First, the window-based ZFST value had to be above the 99.5th percentile of pairwise genetic differentiation. Second, the gene FST value had to be above the 99.5th percentile of pairwise genetic differentiation and last, the gene needed to include at least one SNP with significant differentiation based on Fisher’s exact test (calculated with PoPoolation2 [Kofler et al., 2011b]; p<0.001, Benjamini-Hochberg corrected P-values [Benjamini and Hochberg, 1995]).

Identifying polymorphic sites

Request a detailed protocol

SNP calling was performed with Snape (Raineri et al., 2012). The program requires information of the prior nucleotide diversity θ. Hence, the initial values of nucleotide diversity obtained with PoPoolation were used. Snape was run with folded spectrum and prior type informative. As Snape requires the MPILEUP format, the previously generated MPILEUP files were used. SNP calling was separately performed on coding and non-coding parts of the genome. Therefore, each population MPILEUP file was filtered by coding sequence position with the filter-pileup-by-gtf.pl script of PoPoolation. For coding sequences, the --keep-mode was set to retain all coding sequences. The non-coding sequences were obtained by using the default option and thus discarding the coding sequences from the MPILEUP file. Snape produces a posterior probability of segregation for each position. The posterior probability of segregation was used to filter low-confidence SNPs and indicated in the specific section.

Divergence and polymorphisms in 0-fold and 4-fold sites

Request a detailed protocol

Polarization of synonymous sites (four-fold degenerated sites) and non-synonymous sites (zero-fold degenerated sites) was done using the pseudogenomes of outgroups Nothobranchius orthonotus and Nothobranchius rachovii. For each population the genomic information of the respective pseudogenome was extracted with bedtools getfasta command (Quinlan, 2014; Quinlan and Hall, 2010) and the derived allele frequency of every position was inferred with a custom R script. Briefly, only sites with the bases A, G, T or C in the outgroup pseudogenome were included and checked whether the position has an alternative allele in each of the investigated populations. Positions with an alternative allele present in the population data were treated as possible divergent or polymorphic sites. The derived frequency was determined as frequency of the allele not present in the outgroup. Occasions with an alternate allele present in the population data were treated as possible divergent or polymorphic sites. Divergent sites are positions in the genome were the outgroup allele is different from the allele present in the population. Polymorphic sites are sites in the genome that have more than one allele segregating in the population. Only biallelic polymorphic sites were used in this analysis. The DAF was determined as frequency of the allele not shared with the respective outgroups. In general, positions with only one supporting read for an allele were treated as monomorphic sites. SNPs with a DAF < 5% or>95% were treated as fixed mutations. Further filtering was done based on the threshold of the posterior probability of >0.9 calculated with Snape (see previous subsection), combined with a minimum and maximum coverage threshold per population (GNP: 24, 94; NF414:19, 77; NF303: 18, 70).

Asymptotic McDonald-Kreitman α

Request a detailed protocol

The rate of substitutions that were driven to fixation by positive selection was evaluated with an improved method based on the McDonald-Kreitman test (McDonald and Kreitman, 1991). The test assumes that the proportion of non-synonymous mutations that are neutral has the same fixation rate as synonymous mutations. Therefore, under neutrality the ratio between non-synonymous to synonymous substitutions (Dn/Ds) between species is equal to the ratio of non-synonymous to synonymous polymorphisms within species (Pn/Ps). If positive selection takes place, the ratio between non-synonymous to synonymous substitutions between species is larger than the ratio of non-synonymous to synonymous polymorphism within species (McDonald and Kreitman, 1991). The concept behind this is that the selected variant reaches fixation in a shorter time than by random drift. Therefore, the selected variant increases Dn, not Pn. The proportion of non-synonymous substitutions that were fixed by positive selection (α) was estimated with an extension of the McDonald-Kreitman test (Smith and Eyre-Walker, 2002). Due to the presence of slightly deleterious mutations the estimate of α can be underestimated. For this reason, the method used by Messer and Petrov was implemented to calculate α as a function of the derived allele frequency x (Messer and Petrov, 2013; Haller and Messer, 2017). With this method the true value of α can be inferred as the asymptote of the function of α. Additionally, the value of α(x) for low derived frequencies should give an estimate of the number of slightly deleterious mutations that segregate in the population.

Direction of selection (DoS)

Request a detailed protocol

To further investigate the signature of selection, the direction of selection (DoS) index for every gene was calculated (Stoletzki and Eyre-Walker, 2011). DoS standardizes α to a value between −1 and 1. A positive value of DoS indicates adaptive evolution (positive selection) and a negative value indicates the segregation of slightly deleterious alleles, therefore weaker purifying selection (Stoletzki and Eyre-Walker, 2011). This ratio is undefined for genes without any information about polymorphic or substituted sites. Therefore, only genes with at least one polymorphic and one substituted site were included.

Inference of distribution of fitness effects

Request a detailed protocol

The distribution of fitness effects (DFE) was inferred using the program polyDFE2.0 (Tataru et al., 2017). For this analysis the unfolded site frequency spectra (SFS) of non-synonymous (0-fold) and synonymous sites (4-fold) were projected into 10 chromosomes for each population. Information about the fixed derived sites was included in this analysis (using Nothobranchius orthonotus). PolyDFE2.0 estimates either the full DFE, containing deleterious, neutral and beneficial mutations, or only the deleterious DFE. The best model for each population was obtained using a model testing approach with three different models implemented in PolyDFE2.0 (Model A, B, C). Due to possible biases from erroneous polarization or unknown demography, runs accounting for polarization errors and demography (+eps, +r) were included. Initial parameters were automatically estimated with the –e option, as recommended. To ensure that the parameter space is explored thoroughly, the basin hopping option was applied with a maximum of 500 iterations (-b). The best model for each population was chosen based on the Akaike Information Criterion (AIC). Confidence intervals were generated by running 200 bootstrap datasets with the same parameters used to infer the best model.

Simulations of population demography and the distribution of fitness effects

Request a detailed protocol

We performed simulations with SLiM (version 3.3) (Haller and Messer, 2019), using demographic parameters from the PSMC’ analysis. We simulated four models with two populations that diverged from the same ancestral population and are identical in all population genetic parameters. We adapted the simulation from Rousselle et al., 2018 and simulated 1500 coding sequences of 500 bp length, separated by non-coding regions and re-scaled the population genetic parameter to an initial population size of 10 k (mutation rate of 5.468928e-08, recombination rate of 6.604764e-06). The neutral mutations were initialized in ~1/3 of the coding region and deleterious mutations were initialized in ~2/3 of the coding region, resembling the structure of a codon. Deleterious mutations were simulated with a selection coefficient s, a measure of relative fitness, drawn from the same reflected gamma distribution (mean of −2.5 and shape of 0.35) and a dominance coefficient of 0.1. All models start with an initial burn-in phase of 50000 generations to generate genetic diversity and a stable population. The models vary in either having immediate population size changes followed by constant population size (models A, B) or having exponential growth (models C, D). We further distinguished between a split directly after the burn-in phase (models B, D) or at a later timepoint, following the PSMC’ analysis interpretation (models A, C). The scripts are available on https://github.com/valenzano-lab/simulation_DFE. We performed 25 replicates per model and used a sample size of 30 diploid individuals to retrieve the fixed and segregating mutations. The synonymous and non-synonymous site-frequency-spectra were built by projecting the retrieved mutations into 10 chromosomes, comparable to the site-frequency-spectra generated for the observed data of the turquoise killifish populations and subsequently used as the input for PolyDFE2.0. Model choice of the negative distribution of fitness effects was performed similar to the initial analysis as described previously, except we did not include the ancestral misidentification error.

Variant annotation

Request a detailed protocol

Classification of changes in the coding-sequence (CDS) was done with the variant annotator SnpEFF (Cingolani et al., 2012). The new genome of Nothobranchius furzeri (NFZ v2.0) was implemented to the SnpEFF pipeline. Subsequently, a database for variant annotation with the genome NFZ v2.0 FASTA file and the annotation GTF file was generated. For variant annotation the population specific synonymous and non-synonymous sites with a change in respect to the reference genome NFZ v2.0 were used to infer the impact of these sites. The possible annotation impact classes were low, moderate, and high. SNPs with a frequency below 5% or above 95% were excluded for this analysis. To be consistent with the analysis of the distribution of fitness effects, only positions also found to be present in the N. orthonotus pseudogenome were considered. Positions with warnings in the variant annotation were removed.

Consurf analysis

Request a detailed protocol

The Consurf score was calculated accordingly to the method used in Cui et al., 2019. We used the Consurf (Pupko et al., 2002; Mayrose et al., 2004; Glaser et al., 2003; Ashkenazy et al., 2016) package to assign each AA a conservation score based on the evolutionary rate in homologs of other vertebrates. Consurf scores were estimated for 12575 genes of N. furzeri and synonymous and non-synonymous genomic positions were matched with the derived allele frequency of N. orthonotus and N. rachovii, respectively. The derived frequencies were binned in five bins and we used pairwise Wilcoxon rank sum test to assess significance after correcting for multiple testing (Benjamini and Hochberg adjustment) between each subsequent bin per population and matching bins between populations.

Over-representation analysis

Request a detailed protocol

Gene ontology (GO) and pathway overrepresentation analysis was performed with the online tool ConsensusPathDB (http://cpdb.molgen.mpg.de; version34) (Herwig et al., 2016) using ‘KEGG’ and ‘REACTOME’ databases. Briefly, each gene present in the outlier list was provided with an ENSEMBL human gene identifier (Zerbino et al., 2018), if available, and entered as the target list into the user interface. All genes included in the analysis and with available human ENSEMBL identifier were used as the background gene list. ConsensusPathDB maps the entries to the databases and calculates the enrichment score for each entity by comparing the proportion of target genes in the entity over the proportion of background genes in the entity. For each of the enrichment a P-value is calculated based on a hyper geometric model and is corrected for multiple testing using the false discovery rate (FDR). Only GO terms and pathways with more than two genes were included. Overrepresentation analysis was performed on genes falling below the 2.5th percentile or above the 97.5th percentile thresholds. The percentiles for either FST or DoS values were calculated with the quantile() function in R.

Statistical analysis and data processing

Request a detailed protocol

Statistical analyses were performed using R studio version 1.0.136 (R version 3.3.2 [RStudio, 2015]) on a local computer and R studio version 1.1.456 (R version 3.5.1) in a cluster environment at the Max-Planck-Institute for Biology of Ageing (Cologne). Unless otherwise stated, the functions t.test() and wilcox.test() in R have been used to evaluate statistical significance. To generate a pipeline for data processing we used Snakemake (Köster and Rahmann, 2012). Figure style was modified using Inkscape version 0.92.4. For circular visualization of genomic data we used Circos (Krzywinski et al., 2009).

Inference of demographic population history with individual resequencing data

Request a detailed protocol

To infer the demographic history, we performed whole genome re-sequencing of single individuals from all populations resulting in mean genome coverage between 13-21x (Supplementary file 1B). Demographic history was inferred from single individual sequencing data using Pairwise Sequential Markovian Coalescence (PSMC’ mode from MSMC2 [Schiffels and Durbin, 2014]). Re-sequencing of single individuals was performed with the DNA of single individuals extracted for the pooled sequencing for each examined population. The Illumina short-insert library was constructed based on a published protocol (Rowan et al., 2015). Extracted DNA (500 ng) was digested with fragmentase (New England Biolabs) for 20 min at 37°C, followed by end-repair and A-tailing (1.0 µl NEB End-repair buffer, 0.5 µl Klenow fragment, 0.5 µl Taq.Polymerase, 0.2 µl T4 polynucleotide kinase, 10 µl reaction volume, 30 min at 25°C, 30 min at 75°C) and adapter ligation (NEB Quick ligase buffer 12.5 µl, Quick ligase 0.5 µl, 1 µl adapter P1 (D50X), 1 µl adapter P2 (universal), 5 µM each; 20 min at 20°C, 25 µl reaction volume). Next, ligation mix was diluted to 50 µl and used 0.583:1 vol of home-brewed SPRI beads (SPRI binding buffer: 2.5M NaCL, 20 mM PEG 8000, 10 mM Tris-HCL, 1 mM EDTA,ph = 8, 1 mL TE-washed SpeedMag beads GE Healthcare, 65152105050250 per 100 mL buffer) for purification. The ligation products were amplified with 9 PCR cycles using KAPA Hifi kit (Roche, P5 universal primer and P7 indexed primer D7XX). The samples were pooled and sequenced on Hiseq X. Raw sequencing reads were trimmed using Trimmomatic-0.32 (ILLUMINACLIP:illumina-adaptors.fa:3:7:7:1:true, LEADING:20, TRAILING:20, SLIDINGWINDOW:4:20, MINLEN:50) (Bolger et al., 2014). Data files were inspected with FastQC v0.11.22. Trimmed reads were subsequently mapped to the reference genome with BWA-MEM (version 0.7.12). The SAM output was converted into BAM format, sorted, and indexed via SAMTOOLS v1.3.1 (Li et al., 2009). Filtering and realignment was conducted with PICARD v1.119 and GATK (McKenna et al., 2010). Briefly, the reads were relabeled, sorted, and indexed with AddOrReplaceReadGroups. Duplicated reads were marked with the PICARD feature MarkDuplicates and reads were realigned with first creating a target list with RealignerTargetCreator, second by IndelRealigner from the GATK suite. Resulting reads were again sorted and indexed with SAMTOOLS. Next, the guidance for PSMC’ (https://github.com/stschiff/msmc/blob/master/guide.md) was followed; VCF-files and masked files were generated with the bamCaller.py script (MSMC-tools package). This step requires the chromosome coverage information to mask regions with too low or too high coverage. As recommended in the guidelines, the average coverage per chromosome was calculated using SAMTOOLS. In addition, this step was performed using a coverage threshold of 18 as recommended by Nadachowska-Brzyska et al., 2016. Final input data were generated using the generate_multihetsep.py script (MSMC-tools package). Subsequently, for each sample PSMC’ was run independently. Bootstrapping was performed for 30 samples per individual and input files were generated with the multihetsep_bootstrap.py script (MSMCtools package).

Analysis of differential expressed genes with age

Request a detailed protocol

We downloaded the previously published RNAseq data from a longitudinal study of Nothobranchius furzeri (Reichwald et al., 2015). The data set contains five time points (5 w, 12 w, 20 w, 27 w, 39 w) in three different tissues (liver, brain, skin). The raw reads were mapped to the NFZ v2.0 reference genome and subsequently counted using STAR (version 2.6.0 .c) (Dobin et al., 2013) and FeatureCounts (version 1.6.2) (Liao et al., 2014). We performed statistical analysis of differential expression with age using DESeq2 (Love et al., 2014) and age as factor. Genes are classified as upregulated in young (log(FoldChange)<0, adjusted p<0.01), upregulated in old (log(FoldChange)>0,adjusted p<0.01).

Data availability

Genome sequencing data have been deposited on GenBank under the following accession code: BioProject ID: PRJNA599375 The pooled data sequences used for population genetics under the following accession code: Bioproject PRJNA627180.

The following data sets were generated
    1. Cui R
    2. Willemsen D
    3. Valenzano DR
    (2020) NCBI BioProject
    ID PRJNA599375. Nothobranchius furzeri assembly 2.0.
    https://www.ncbi.nlm.nih.gov/bioproject/PRJNA599375/
    1. Willemsen D
    2. Valenzano DR
    (2020) NCBI BioProject
    ID PRJNA627180. Pooled whole genome sequencing and individual resequencing of wild-caught Nothobranchius furzeri.
    https://www.ncbi.nlm.nih.gov/bioproject/PRJNA627180/

References

    1. Charlesworth B
    (2000)
    Medawar, Hamilton and the evolution of aging
    Genetics 156:927–931.
  44. Software
    1. Lehner B
    2. Verdin K
    (2006) HydroSHEDS Technical Documentation. World Wildlife Fund US, Washington
    HydroSHEDS Technical Documentation. World Wildlife Fund US, Washington.
    http://hydrosheds.cr.usgs.gov
  73. Software
    1. RStudio
    (2015)
    RStudio: Integrated Development for R
    RStudio: Integrated Development for R.

Decision letter

  1. Pankaj Kapahi
    Reviewing Editor; Buck Institute for Research on Aging, United States
  2. Jessica K Tyler
    Senior Editor; Weill Cornell Medicine, United States
  3. Nathan Clark
    Reviewer; Utah, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

Your study on using turquoise killifish to determine the evolutionary forces shaping life history divergence within species is of great interest to the field of aging. The use of genome sequencing and population genetics allows to test some of the major evolutionary theories of aging. Your findings that small population size and genetic drift determine the accumulation of deleterious gene variants in specific genes leading to short lifespan will have great impact on the field bringing those who think of evolutionary and molecular mechanisms of aging closer.

Decision letter after peer review:

Thank you for submitting your article "Intra-Species differences in population size shape life history and genome evolution" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Jessica Tyler as the Senior Editor.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

As the editors have judged that your manuscript is of interest, but as described below that additional experiments are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is 'in revision at eLife'. Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)

In addition to a point by point rebuttal, in particular please respond to the following points:

1) The major conclusions are based on too narrow an interpretation of the senescence literature. They have not considered mutation accumulation as a basis for their results. They need to do so and to adjust the rest of the paper accordingly or at least need to explain to me why they did not include it.

2) The authors should explain the inferred DFEs and show that the DFE inference method to produce reasonable results on simulated data with a similar demographic history to what they infer from the killifish data.

“There is also concern about the DFE study, since, if I understand Figure 4D correctly, the wetter habitat populations have a higher proportion of strongly deleterious SNPs than the GNP (dry) population. Does this go against the argument presented in the paper?”

Reviewer #1:

In their manuscript, the authors investigate the evolutionary forces shaping lifespan in killifish populations. They produced an improved genome assembly and annotation and generate population genetic data for 3 populations of the species. By contrasting the relative wetness of the habitats and using an array of population genetic analyses, they make statements about their effective population sizes, loci with signatures of positive selection (extreme differentiation), and contrasts in their distributions of fitness effects (DFEs). They conclude that the QTL loci previously linked to differences in lifespan did not show signs of positive selection and that genetic drift has allowed a relatively large proportion of slightly deleterious polymorphisms to persist in the smaller population with the most arid habitat.

Generally, I found the study interesting and the choice of system an important one since this model is being used for longevity research. The conclusion that lifespan differences between populations is driven mainly by drift has implications for how results from this model are interpreted in the future. There are a few major concerns which I believe are related to this conclusion. They are detailed below.

1) The study ultimately cannot comment on the initial ancestral evolution of shorter time to reproductive maturity or lifespan because it is analyzing differences between populations, all of which already have short lifespans and time to maturity. This is not a critical concern however, and the authors seems to be aware of this, but it should be stated upfront so readers appreciate the distinction. It remains a possibility that the species-wide trait of short lifespan and time to reproductive maturity arose through selection in response to the extreme seasonality of their habitat.

2) It is not entirely clear whether the populations studied here are the same as those in which the lifespan QTL were discovered. What are the genetic and geographic distances between these populations and those studied by Valenzano et al.

If the populations are substantially different, then the finding of no evidence of positive selection at the QTL is limited, since lifespan associated QTL in these populations could be entirely different. This is an important concern to address.

3) The reader would benefit from explanations about the expectations for the α statistic, since many will not have encountered it before. What are expectations under relevant scenarios? What do values of 0, 1, and negative really mean? How robust are such inferences to demographics and other confounders?

Because conclusions of that paper are somewhat reliant on these values and related tests about the DFEs and DoS, it is important to spell this out and to understand their interpretation, strengths and weaknesses.

4) Each of the population genetic tests could be affected by demographic changes. The authors should address this possibility head on and explain how the method or its interpretation is or is not affected by the scenarios they inferred in Figure 1A.

5) The DFE analysis in Figure 4D seems to suggest a very high proportion of highly deleterious SNPs. ~60% of SNPs with a value of 4Nes < -100. This value seems high. Wouldn't we expect most segregating SNPs to be nearly neutral?

Could the authors explain? Is this high proportion of highly deleterious SNPs due to SNPs at very low frequency?

6) Given point 5 above, the wet population has a much higher proportion of highly deleterious SNPs than the arid population. Couldn't it be argued that the wet population is suffering poorer genetic health, since it has a higher proportion of highly damaging SNPs?

The Consurf analysis seems more convincing to me than the DFE analysis.

Reviewer #2:

The authors present a diversity of analyses of patterns of genetic variation among populations of the annual killifish Nothobranchus fuerzi. These fish live in ponds that are dry for large portions of every year. They are able to persist because they lay eggs that are desiccation resistant, but also because the adults have very rapid development, enabling them to mature quickly and reproduce during the brief times that the ponds hold water. The adults also age rapidly and have short life expectancies. The different populations represented in the study come from ponds distributed along a gradient of annual rainfall, so some ponds are destined to retain water for much shorter intervals of time than others. The fish from ponds which are expected to have short durations do not have faster development in the laboratory, which is surprising, yet they do have faster aging and shorter life expectancies. This paper is one of a series in which the authors argue that the real reason these fish have shorter life expectancies is that these ponds of smaller effective population sizes and hence are more susceptible to genetic drift.

In offering this review, I must also acknowledge my background and limitations. I am familiar with the biology of these fish, with the empirical study of aging and with the literature on the evolution of aging. I also have a good background in evolutionary biology. I am not a practitioner of the sorts of genetic methods employed in this study and am less able to judge them.

I have two categories of concern about this paper. The first is that I feel it is poorly written and inappropriately written for the general audience targeted by a journal like eLife. Second, I feel that there are critical deficiencies in the science.

Presentation: This work is actually a test of alternative hypotheses, yet the authors never clearly state this. They instead present an argument in favor of one of the two hypotheses while the contrast between the two is implicit. The two hypotheses are: 1) the differences in life expectancy among the various populations are a consequence of how lifespan evolved in each of them, 2) the shorter lifespans of fish from ponds of shorter duration are a by-product of genetic drift.

Hypothesis 1: There are diverse hypotheses for the evolution of lifespan and aging, but I will only address three of them here. The two that date to the 1950's are "antagonistic pleiotropy" (Williams, 1957) and "mutation accumulation" (Medawar, 1952). Medawar's hypothesis is relevant here, since he predicts that species/populations that experience shorter life expectancy will more readily accumulate deleterious mutations that act late in life because natural selection (purifying selection) will weaken with age as a consequence of the small number of individuals that survive to advanced ages.

In supporting their argument for the importance of genetic drift, the authors show an accumulation of genetic variation in the population that presumably has reduced life expectancy plus they present evidence that this variation is likely to be deleterious. A problem I see with their argument is that mutation accumulation makes the same prediction, save for the restriction that the deleterious mutations be ones that are age-specific in their effects. The fact is that they do not know whether or not the variation they characterize is age-specific in its expression. I cannot see how they can discriminate between deleterious mutations that accumulate because of genetic drift versus those that accumulate because of mutation accumulation. The only compelling argument I do see is their evidence for effective population size, but the low end of the estimated population size for the more rapidly aging population is in the vicinity of 20,000, which is not that small and not by itself likely to generate severe genetic drift. Furthermore, this population has been declining in size over time. My understanding of the expected trend for long term declines in population size (based on the application of this idea to whales by Steve Palumbi) is that there will be more genetic variation than expected, given the current population size.

A second implicit argument for genetic drift is the results of their 2017 paper in Evolution in which they quantify other aspects of the life histories of multiple populations of this species. This paper is a lab study of the life histories of multiple populations of a number of different species of annual killifish, one of which is N. fuerzi. The study is good in terms of the spectrum of life history variables they assess. A noteworthy feature of the results for N. fuerzi is the differences among populations in lifespan in a controlled laboratory environment (populations from ponds with shorter expected duration also tend to have shorter lifespans) without any clear evidence of a tradeoff in other aspects of the life history (e.g., those populations that are shorter lived also mature at an earlier age and/or invest more in offspring). In fact, the short-lived population did tend to invest more in eggs early in life, but the difference was relatively small. As an aside, they put some effort into addressing a third hypothesis for the evolution of aging, which is condition-dependent selection (Williams and Day, 2003). I think that condition dependence is an unlikely mechanism in this system. It can be effective if short life expectancy is due to predation, for example, because adaptation to predation can indirectly affect aging. When a pool dries, there is no option for adaptation in a fish – they all die.

We can agree that there is not strong support for any form of antagonistic pleiotropy revealed in the 2017 paper. The best shot is a slightly higher level of reproductive investment in the population from a shorter-lived pond, which would balance against reduced lifespan. I would not feel tempted to make this argument, but I would advise not banking too heavily on this study fully eliminating antagonistic pleiotropy as a source of differences among populations in lifespan. I would be inclined to first take another shot at this approach.

In summary, they need to be explicit about the alternative hypotheses and need to address the possible role of mutation accumulation. Note that there is empirical evidence in favor of a contribution of mutation accumulation to aging. In many other regards, their presentation is not one that will readily appeal to a general audience. I found their figures to be particularly difficult to interpret.

Science: I have confounded presentation with science a bit since my argument about mutation accumulation applies to the scientific aspects of the presentation, but there are other issues as well. I summarize some of these in the following comments on specific portions of the manuscript.

Other comments:

Subsection “Population genetics of natural turquoise killifish populations”: If they are working with lab lines, then how can we know whether any genetic drift associated with the different lines represents differences among natural populations vs. drift that occurred during prolonged laboratory culture. The effective size of the laboratory lines and their duration in captivity matter. If the genomes were just from wild-caught fish then this is not an issue, but it could apply to the results reported in the 2017 Evolution paper.

Another issue is to better illustrate the climatic differences among the study populations. I can see that they fall at different places in a rainfall gradient, but it would be better to give a more concrete comparison, like the duration of holding water. Are any such comparisons available?

Subsection “Genetic differentiation among turquoise killifish populations”: Are positive natural selection and purifying selection the only interpretations for the FST values? I had thought that various forms of balancing selection can yield the same results that you attribute to purifying selection (with regard to purifying selection – shouldn't this interpretation be broadened to include balancing selection?).

This section is a bit deceptive in presentation. These regions are all ones that were not associated with any QTL for lifespan, correct? This means that you are looking for age-specific differences in expression in the absence of any affiliation with lifespan, save one exception. I understand the bottom line, which is that you will suggest that the absence of any positive evidence for a genetic basis in lifespan the alternative of genetic drift becomes more plausible. You can do a much better job of making these alternatives more explicit. But, there are very good reasons why there could be portions of the genome under selection that were not revealed by your QTL study. Consider adding a power analysis of these results if you really want to make this argument; e.g., how big could an effect have been yet not show up as significant? What if the evolution of aging and lifespan were polygenic, with each gene making a small contribution? Positive results from QTL can be very informative, but I do not thing the reverse is true.

Delete subsection “Evolutionary origin of the sex chromosome”. It has nothing to do with the central theme of the paper, so it is an unnecessary digression. It also impressed me as a topic that could be expanded on and published on its own.

Subsection “Relaxed selection in turquoise killifish populations”: I understand the goal of this section of the paper, which is to make a compelling case for genetic drift as the cause of differences among populations in lifespan. This is the part of the paper that I am not well qualified to judge so I hope other reviewers can do better. But, I can at least question whether or not mutation accumulation could be a plausible alternative, since I believe it is consistent with the patterns reported here. If it is not, then you need to take it on then explain why not.

Subsection “Relaxation of selection in age-related disease pathways”: My reaction to this section is the same as the last section. I am limited in my ability to evaluate this, but I also question whether mutation accumulation represents a viable alternative explanation.

Discussion paragraph two: Eliminate these lines and save them for a later paper on the sex chromosome.

Discussion: I think that the Discussion must be redone in concert with the Introduction. The Introduction should formally present the alternative explanations for differences among populations in aging – antagonistic pleiotropy, mutation accumulation, drift. I think it is okay to just refer to the 2017 paper to dismiss other alternatives, like condition dependent selection. I do not see how you can discriminate between drift and mutation accumulation. I am particularly concerned by your reporting a reasonably large effective population size for the population from the driest environment, in spite of its being the smallest population. If you really want to go this route, then I think it is incumbent on you to do some simulations, or at least something concrete, to argue that the population size is small enough for drift to have been an issue. These issued would be subjects for the Discussion.

Subsection “Koeppen-Geiger index and bioclimatic variables”: This is an obscure way of reporting differences in local climate and pond longevity. It would make a big difference if you could give me something more concrete than this since the values you report look small and subtle but the implied effects on biology are large. The two do not fit well together.

Reviewer #3:

This paper presents a new killifish genome assembly as well as population sequencing data from several killifish populations living in wet and dry habitats. Using this dataset, the authors show that killifish lifespan is inversely correlated with several measures of population diversity. They find no evidence for positive selection favoring genetic variation that appears to shorten lifespan based on a prior QTL analysis, but instead find support for the idea that increased genetic drift could have caused lifespan-shortening variation to reach high frequency in small, outlying killifish populations. Overall the study is very nice and provides some compelling evidence that killifish lifespan has been shortened by weakly deleterious mutation accumulation.

1) My one comprehensive issue with the manuscript is that it implies too strongly that the dryness of the outlying killifish habitats directly caused those populations to lose diversity and accumulate deleterious mutations. Even if all killifish habitats had identical rainfall and resources, we would still expect the populations at the origin of the range expansion to have higher diversity and less genetic load than smaller, more recently founded, outlying groups. In the absence of any data from control pairs of killifish populations that are geographically distant from each other but have similar levels of rainfall, the paper should clearly state up front that the differences in diversity between populations might not be driven primarily by wetness vs dryness but would be expected to exist whenever one population is more recently founded than another. In addition to addressing this broad issue, it would be helpful for the authors to address the following more specific issues:

2) The manuscript notes that FST outlier genes are differentially expressed between fish of different ages but does not comment on how often random control genes likewise vary in expression as a function of age.

3) For the gene with unusually low FST that is located on the sex chromosome, is its outlier status appropriately calibrated with respect to the overall divergence and diversity levels of the sex chromosomes, which can differ substantially from that of the autosomes?

4) In this discussion of the origin of the killifish sex chromosome, it wasn't clear to me whether the sex determining locus was inherited ancestrally or seems to have acquired a de novo sex determining function.

5) How good is the evidence that the DFE inference method being used is robust to non-equilibrium demographic history?

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "Intra-Species differences in population size shape life history and genome evolution" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Jessica Tyler as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

As the editors have judged that your manuscript is of interest, but as described below that additional experiments are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is "in revision at eLife". Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)

In this study the authors use turquoise killifish to determine the evolutionary forces shaping life history divergence within species. Their results support that mutation accumulation plays an important role in shaping lifespan, but the reviewers have also raised concerns that must be addressed before acceptance of this manuscript.

Essential revisions:

1) The authors still do not state how close these populations are to those studied in the QTL study (Kirschner et al., 2012). Couldn't they just provide coordinates, a map, physical distances, or something to help the reader understand how close they are? As I previously stated, if they are not close, then it raises issues of whether these populations would even have the same loci affecting differences in lifespan. The authors provide the new statements below, but it is not clear how to interpret them objectively. Were these fish collected meters from the previous ones? Kilometers? 100 km?

"The populations used in this study are from localities in proximity to those used in the previous QTL study."

"we collected wild populations from natural localities as close as possible to the locations of origin of the strains used for the QTL study."

The revision also adds:

"…or that the populations used in this study and those used for the QTL analysis had a different genetic architecture of lifespan"

This statement could be unnecessary if exact locations show them to be quite close.

2) The explanation of MK α is still inadequate. Some PopGen aficionados will know, but general eLife readers won't know what these values mean. What are high values, low values… etc. What has changed in the sequences to cause these values to change? The authors must really present these statistics in a didactic way.

3) From the new explanations of the inferred DFE of new mutations, it seems this analysis only reflects inferred population demographics, rather than the actual deleteriousness of polymorphisms in the population. Is it correct to say that it uses inferred demographics to just rescale a statistical distribution chosen to represent the DFE?

If so, why not then just present the inferred demographic history? If the DFE analysis isn't analyzing anything about the actual polymorphisms present in the populations, why make that extra inferential jump? Just report the demographics. If this analysis reflects some other quality of the data that translates into actual fitness effects, then it should be stated. If it reflects demographics alone as they would affect a theoretical DFE, then that should be clearly stated.

4) The authors have done a thorough job with revisions and substantially improved their manuscript. The one point from my review that I'd like to revisit is my point that more recently founded populations tend to have smaller effective size and less diversity than older populations regardless of founding order. I don't agree with the reasoning that this only applies to temporally separate populations-there is extensive literature in humans on a "serial bottleneck model" that explains why population diversity is proportional to the distance of a population from the origin of the human range expansion in Africa (see e.g. Ramachandran et al., 2005). This work was based on genetic sampling of human populations at the same time, as was done with the killifish populations here, but founding order is important because more recently founded populations have had less time to recover from the associated bottlenecking founder effect. I don't think it's necessary for the authors to do more work to address this point, but I do think they should point out in the text that some proportion of the population size difference between the wet and dry populations is probably due to the fact that killifish have existed at large population size in the wet habitats for a long time, whereas the dry population probably expanded from a small founder population more recently.

https://doi.org/10.7554/eLife.55794.sa1

Author response

In addition to a point by point rebuttal, in particular please respond to the following points:

1) The major conclusions are based on too narrow an interpretation of the senescence literature. They have not considered mutation accumulation as a basis for their results. They need to do so and to adjust the rest of the paper accordingly or at least need to explain to me why they did not include it.

We thank the reviewer to raise this important point. Our findings indeed lend full support to mutation accumulation being the main driver of senescence evolution in our empirical example of natural turquoise killifish populations. Indeed, we believe that our work provides a direct link between the classical theories of ageing to empirical findings in ecology and evolution, connecting the mutation accumulation theory with ecological and demographic constraints (habitat aridity and population size, respectively) in wild species. Our main result is that genetic drift due to small effective population size exacerbates the effects of mutation accumulation. We agree that this was not clearly stated in the text and have now explicitly spelled out in Abstract, Introduction and Discussion.

2) The authors should explain the inferred DFEs and show that the DFE inference method to produce reasonable results on simulated data with a similar demographic history to what they infer from the killifish data.

“There is also concern about the DFE study, since, if I understand Figure 4D correctly, the wetter habitat populations have a higher proportion of strongly deleterious SNPs than the GNP (dry) population. Does this go against the argument presented in the paper?”

We fully agree with the editor and with the reviewers that the DFE analysis demands clearer explanation.

We would like to point out that the difference between the DFE spectra between dry and wet populations is along the expected direction. DFE refers to the fitness effects of newly arising mutations. When purifying selection is strong, more new mutations will be strongly deleterious and thus be removed effectively from the population. That is, the strongly deleterious mutations are not observed in the data. However, when purifying selection is weaker, some of the deleterious mutations “move” to the weakly deleterious category and segregate in the population. This is why the population with the smaller effective population size (dry population) has less strongly deleterious new mutations, but relatively more slightly deleterious mutations. Another way to see this point is that mutations with same negative selection coefficient s would be more deleterious (more to the left in the DFE plot) in larger populations than in smaller ones, as 4*Ne*s scales linearly with population size.

To address this point numerically, we performed a new simulation with SLiM (version 3.3, Haller and Messer, 2019), using similar demographic parameters to those from the PSMC’ analysis. We simulated four models with two populations that diverged from the same ancestral population and are identical in all population genetic parameters. Importantly, the selection coefficient s of the deleterious mutations was drawn from the same reflected gamma distribution (mean of -2.5 and shape of 0.35). All models start with an initial burn-in phase of 50000 generations to generate genetic diversity and a stable ancestral population. We re-scaled the population genetic parameters to an initial population size of 10k to increase the speed of the simulation. This should not affect the results and is a common method. The models vary in either having immediate population size changes (Model A, B) or having exponential growth (Model C, D). We further distinguished between a split directly after the burn-in phase (Model B, D) or at a later timepoint, following the PSMC’ analysis interpretation (Model A, C).

The output of the simulation was used to build the site-frequency-spectra that were used as the input for PolyDFE. We found that the PolyDFE results recovered the same conclusions as those we empirically observed between the large and the small population. In fact, although the selection coefficient for the simulations was drawn from the same gamma distribution for large and small populations, the small population evolved fewer strongly deleterious novel mutations, but more slightly deleterious mutations, as expected under weaker purifying selection. The results of this simulation are displayed in Figure 4—figure supplement 1.

Reviewer #1:

[…]

1) The study ultimately cannot comment on the initial ancestral evolution of shorter time to reproductive maturity or lifespan because it is analyzing differences between populations, all of which already have short lifespans and time to maturity. This is not a critical concern however, and the authors seems to be aware of this, but it should be stated upfront so readers appreciate the distinction. It remains a possibility that the species-wide trait of short lifespan and time to reproductive maturity arose through selection in response to the extreme seasonality of their habitat.

As correctly noted by the reviewer, the current study specifically focuses on population differences within the model species Nothobranchius furzeri. Certainly, since all Nothobranchius share common adaptations, including embryonic diapause and rapid developmental time, we cannot use the current dataset to search for genes important for genus-wide adaptations that emerged in the common ancestor of Nothobranchius. To note, different species in the genus Nothobranchius differ in lifespan – with species living up to 2 years – even though they have embryonic diapause and reach sexual maturity roughly at the same time. Similarly, different wild populations of the species Nothobranchius furzeri reach sexual maturation at the same time, even though have different lifespan in nature and in captivity. Hence, early life history adaptations (embryonic diapause and rapid sexual maturation) appear somehow disjointed in the genus Nothobranchius and in the species Nothobranchius furzeri. We believe that African killifish are an example of how early life history traits may be under different selective constraints and evolutionary trajectories than late life history traits. This very aspect was addressed in a previous work conducted by our group and published in 2019 (Cui et al., 2019). The present work is intended to give insights to the mechanisms that resulted in lifespan differences between populations of the turquoise killifish. We now clarify in the text that our study does not address the evolution of rapid sexual maturation and diapause in killifish, but rather provide an evidence of how nearly neutral evolution adds genome-wide deleterious gene variants that predominantly affect late life history traits.

2) It is not entirely clear whether the populations studied here are the same as those in which the lifespan QTL were discovered. What are the genetic and geographic distances between these populations and those studied by Valenzano et al.

If the populations are substantially different, then the finding of no evidence of positive selection at the QTL is limited, since lifespan associated QTL in these populations could be entirely different. This is an important concern to address.

This point is well taken. Although the strains used in Valenzano et al. are derived from wild populations that differ in their lifespan, we cannot assume that the lifespan QTL architecture is shared across any pair of long- vs. short-lived turquoise killifish populations. This very question is what we wanted to explore when we first started this study, and therefore we collected wild populations from natural localities as close as possible to the locations of origin of the strains used for the QTL study. We clarify this point now in results.

3) The reader would benefit from explanations about the expectations for the α statistic, since many will not have encountered it before. What are expectations under relevant scenarios? What do values of 0, 1, and negative really mean? How robust are such inferences to demographics and other confounders?

Because conclusions of that paper are somewhat reliant on these values and related tests about the DFEs and DoS, it is important to spell this out and to understand their interpretation, strengths and weaknesses.

We fully agree with the reviewer and have now added explanations of what each of these statistics means, what are the expectations, and how our results fit such expectations.

4) Each of the population genetic tests could be affected by demographic changes. The authors should address this possibility head on and explain how the method or its interpretation is or is not affected by the scenarios they inferred in Figure 1A.

PolyDFE fully accounts for demographic changes (population size is in the x axis of each plot). Fluctuations in demography are cancelled out when computing DoS and MK α, a desirable feature of these statistics. Indeed, the main argument of our work is that demographic changes, affecting drift, lead to differences in genetic load. Rather than a confounder, we believe that demographic changes are a direct cause of the observed differences among turquoise killifish populations.

5) The DFE analysis in Figure 4D seems to suggest a very high proportion of highly deleterious SNPs. ~60% of SNPs with a value of 4Nes < -100. This value seems high. Wouldn't we expect most segregating SNPs to be nearly neutral?

Could the authors explain? Is this high proportion of highly deleterious SNPs due to SNPs at very low frequency?

We realize the DFE analysis needed more explanation and we are thankful to the reviewer for requesting more clarity on this part. DFE computes the fitness effect of novel mutations. Mutations with similar (e.g. negative) selection coefficient s would have a more deleterious effect in large populations than in small populations, since 4*s*Ne is bigger as Ne gets larger. Hence, DFE interprets how purifying selection purges mutations given their selection coefficient and the population size. We have now explained this part extensively and performed a novel simulation with SLiM3 to model how DFE differs in a small and a large population evolved from the same ancestral population. The results of this simulation support our findings and interpretations and are now presented as Figure 4—figure supplement 1. We thank again the reviewer for this note, as we feel that our finding was strengthened by the simulation.

6) Given point 5 above, the wet population has a much higher proportion of highly deleterious SNPs than the arid population. Couldn't it be argued that the wet population is suffering poorer genetic health, since it has a higher proportion of highly damaging SNPs?

See response to point 5). More novel mutations are strongly deleterious in the wet population, and thus they are removed very effectively by purifying selection.

The Consurf analysis seems more convincing to me than the DFE analysis.

We appreciate that the reviewer judges that the Consurf analysis is convincing, and we hope that now also the expanded explanation of the DFE analysis strengthens our findings.

Reviewer #2:

[…[

I have two categories of concern about this paper. The first is that I feel it is poorly written and inappropriately written for the general audience targeted by a journal like eLife. Second, I feel that there are critical deficiencies in the science.

Presentation: This work is actually a test of alternative hypotheses, yet the authors never clearly state this. They instead present an argument in favor of one of the two hypotheses while the contrast between the two is implicit. The two hypotheses are: 1) the differences in life expectancy among the various populations are a consequence of how lifespan evolved in each of them, 2) the shorter lifespans of fish from ponds of shorter duration are a by-product of genetic drift.

Hypothesis 1: There are diverse hypotheses for the evolution of lifespan and aging, but I will only address three of them here. The two that date to the 1950's are "antagonistic pleiotropy" (Williams, 1957) and "mutation accumulation" (Medawar, 1952). Medawar's hypothesis is relevant here, since he predicts that species/populations that experience shorter life expectancy will more readily accumulate deleterious mutations that act late in life because natural selection (purifying selection) will weaken with age as a consequence of the small number of individuals that survive to advanced ages.

In supporting their argument for the importance of genetic drift, the authors show an accumulation of genetic variation in the population that presumably has reduced life expectancy plus they present evidence that this variation is likely to be deleterious. A problem I see with their argument is that mutation accumulation makes the same prediction, save for the restriction that the deleterious mutations be ones that are age-specific in their effects. The fact is that they do not know whether or not the variation they characterize is age-specific in its expression. I cannot see how they can discriminate between deleterious mutations that accumulate because of genetic drift versus those that accumulate because of mutation accumulation.

We thank the reviewer for emphasizing that our manuscript requires a language that is more geared towards a general audience and that we need to be more explicit about the hypotheses that are tested. We hope we now address both these fundamental aspects.

Indeed, our finding that reduced population size in populations from dry environments leads to genome-wide accumulation of deleterious mutations fully echoes Medawar’s mutation accumulation. Our statement is that genetic drift exacerbates the effects of mutation accumulation and that, based on a population genetics approach, we do not find strong support for antagonistic pleiotropism being the main driving evolutionary mechanism leading to short lifespan in killifish populations from dry (hence more ephemeral) environments. Furthermore, in previous work done in our group (Cui et al., 2019), we showed that in annual killifish, genes with higher expression in early life (e.g. developmental genes) are under stronger purifying selection – i.e. accumulate fewer deleterious alleles – than those expressed in late life. In the present work, studying the shortest-known living killifish (turquoise killifish), which is a unique example of species having populations that naturally differ in wild and captive lifespan, we used population genetics to identify the evolutionary forces that have moulded the genome. Indeed, we find that small population size is strictly associated with the genome-wide expansion of slightly deleterious gene variants. Since our main conclusion was not clearly stated in the text, we have now explicitly added key sentences in Abstract, Introduction and Discussion.

The only compelling argument I do see is their evidence for effective population size, but the low end of the estimated population size for the more rapidly aging population is in the vicinity of 20,000, which is not that small and not by itself likely to generate severe genetic drift. Furthermore, this population has been declining in size over time. My understanding of the expected trend for long term declines in population size (based on the application of this idea to whales by Steve Palumbi) is that there will be more genetic variation than expected, given the current population size.

The “small population size” of the population from drier habitat has to be considered relative to the populations from more wet habitats. On a relative scale, smaller populations are expected to experience higher drift. We refrain from making a statement on the absolute scale of population size and drift.

A second implicit argument for genetic drift is the results of their 2017 paper in Evolution in which they quantify other aspects of the life histories of multiple populations of this species. This paper is a lab study of the life histories of multiple populations of a number of different species of annual killifish, one of which is N. fuerzi. The study is good in terms of the spectrum of life history variables they assess. A noteworthy feature of the results for N. fuerzi is the differences among populations in lifespan in a controlled laboratory environment (populations from ponds with shorter expected duration also tend to have shorter lifespans) without any clear evidence of a tradeoff in other aspects of the life history (e.g., those populations that are shorter lived also mature at an earlier age and/or invest more in offspring). In fact, the short-lived population did tend to invest more in eggs early in life, but the difference was relatively small. As an aside, they put some effort into addressing a third hypothesis for the evolution of aging, which is condition-dependent selection (Williams and Day, 2003). I think that condition dependence is an unlikely mechanism in this system. It can be effective if short life expectancy is due to predation, for example, because adaptation to predation can indirectly affect aging. When a pool dries, there is no option for adaptation in a fish – they all die.

We can agree that there is not strong support for any form of antagonistic pleiotropy revealed in the 2017 paper. The best shot is a slightly higher level of reproductive investment in the population from a shorter-lived pond, which would balance against reduced lifespan. I would not feel tempted to make this argument, but I would advise not banking too heavily on this study fully eliminating antagonistic pleiotropy as a source of differences among populations in lifespan. I would be inclined to first take another shot at this approach.

We thank the reviewer for the input. We now cite these studies in our Discussion. Indeed, we are unable to reject antagonistic pleiotropy, even though we could not find good support for it. This may be partly due to the difficulty of detecting positive selection from genomic data in general.

In summary, they need to be explicit about the alternative hypotheses and need to address the possible role of mutation accumulation. Note that there is empirical evidence in favor of a contribution of mutation accumulation to aging. In many other regards, their presentation is not one that will readily appeal to a general audience. I found their figures to be particularly difficult to interpret.

Science: I have confounded presentation with science a bit since my argument about mutation accumulation applies to the scientific aspects of the presentation, but there are other issues as well. I summarize some of these in the following comments on specific portions of the manuscript.

Other comments:

Subsection “Population genetics of natural turquoise killifish populations”: If they are working with lab lines, then how can we know whether any genetic drift associated with the different lines represents differences among natural populations vs. drift that occurred during prolonged laboratory culture. The effective size of the laboratory lines and their duration in captivity matter. If the genomes were just from wild-caught fish then this is not an issue, but it could apply to the results reported in the 2017 Evolution paper.

We would like to clarify that the samples used in this study were all wild caught.

Another issue is to better illustrate the climatic differences among the study populations. I can see that they fall at different places in a rainfall gradient, but it would be better to give a more concrete comparison, like the duration of holding water. Are any such comparisons available?

We agree that the duration of holding water is an important factor and would be a very valuable resource. Satellite data from NASA indeed showed that the duration of holding water is longer in the populations in Mozambique compared to the population in Zimbabwe, but the raw data is not accessible to us. We are confident that the provided environmental data is sufficient to illustrate the differences in the habitat of the investigated samples. Furthermore, a paper published by one of the authors in 2008 (Terzibasi et al., 2008), used as a proxy for (inverse) duration of holding water in this region, the ratio of monthly evaporation over precipitation.

Subsection “Genetic differentiation among turquoise killifish populations”: Are positive natural selection and purifying selection the only interpretations for the FST values? I had thought that various forms of balancing selection can yield the same results that you attribute to purifying selection (with regard to purifying selection – shouldn't this interpretation be broadened to include balancing selection?).

We thank the reviewer for requesting this clarification. Balancing selection is predicted to cause low genetic differentiation between the populations in most of the cases (e.g. Brandt et al., 2017). Only if none of the alleles are shared by both populations, balancing selection can cause a high FST values. Importantly, balancing selection causes a higher genetic diversity compared to the background genomic diversity. We could not find elevated pairwise nucleotide diversity in the FST outlier regions, which would support balancing selection. We now included this in the text.

This section is a bit deceptive in presentation. These regions are all ones that were not associated with any QTL for lifespan, correct? This means that you are looking for age-specific differences in expression in the absence of any affiliation with lifespan, save one exception. I understand the bottom line, which is that you will suggest that the absence of any positive evidence for a genetic basis in lifespan the alternative of genetic drift becomes more plausible. You can do a much better job of making these alternatives more explicit. But, there are very good reasons why there could be portions of the genome under selection that were not revealed by your QTL study. Consider adding a power analysis of these results if you really want to make this argument; e.g., how big could an effect have been yet not show up as significant? What if the evolution of aging and lifespan were polygenic, with each gene making a small contribution? Positive results from QTL can be very informative, but I do not thing the reverse is true.

We agree that we cannot be certain that positive selected genes would show up in the previous QTL results. However, our population genetic analysis showed overall little evidence for positive selection. We now discuss the possibility of other explanations in subsection “Genetic differentiation among turquoise killifish populations”.

Delete subsection “Evolutionary origin of the sex chromosome”. It has nothing to do with the central theme of the paper, so it is an unnecessary digression. It also impressed me as a topic that could be expanded on and published on its own.

We thank the reviewer for showing general interest in the findings, despite questioning the suitability of this result to this study. We agree that this part does not belong to the relaxed selection part of the study, but as the sex determining region is under low genetic divergence between the populations and we provide a new general genome resource, we believe that this result provides a valuable addition to this paper.

Subsection “Relaxed selection in turquoise killifish populations”: I understand the goal of this section of the paper, which is to make a compelling case for genetic drift as the cause of differences among populations in lifespan. This is the part of the paper that I am not well qualified to judge so I hope other reviewers can do better. But, I can at least question whether or not mutation accumulation could be a plausible alternative, since I believe it is consistent with the patterns reported here. If it is not, then you need to take it on then explain why not.

We thank the reviewer for explicitly stating the main implication of our study, which is that mutation accumulation is fully compatible with our findings. Indeed, our results fully support mutation accumulation as a major force shaping genome-wide distribution of deleterious gene variants. We now state this several times through the revised manuscript.

Subsection “Relaxation of selection in age-related disease pathways”: My reaction to this section is the same as the last section. I am limited in my ability to evaluate this, but I also question whether mutation accumulation represents a viable alternative explanation.

Mutation accumulation is indeed supported.

Discussion paragraph two: Eliminate these lines and save them for a later paper on the sex chromosome.

Discussion: I think that the Discussion must be redone in concert with the Introduction. The Introduction should formally present the alternative explanations for differences among populations in aging – antagonistic pleiotropy, mutation accumulation, drift. I think it is okay to just refer to the 2017 paper to dismiss other alternatives, like condition dependent selection. I do not see how you can discriminate between drift and mutation accumulation.

We amended the Introduction and Discussion and included the hypotheses of evolution of ageing. We think that mutation accumulation and drift work hand in hand. When drift increases, slightly deleterious mutations are more likely to become highly frequent. Our study links genetic drift to the mutation accumulation theory of ageing, combining the nearly neutral theory of molecular evolution with the mutation accumulation theory of ageing. We now make this point clearer throughout the revised manuscript.

I am particularly concerned by your reporting a reasonably large effective population size for the population from the driest environment, in spite of its being the smallest population. If you really want to go this route, then I think it is incumbent on you to do some simulations, or at least something concrete, to argue that the population size is small enough for drift to have been an issue. These issued would be subjects for the Discussion.

The level of drift is always relative to the large population. It is not the absolute level of drift per se that matters, but the proportion of bad alleles that becomes effectively neutral that matters. Our study highlights that the difference in population size between the population could be plausible for the accumulation of deleterious mutations in the dry population. To address this part, we performed several simulations with SLiM and added Figure 4—figure supplement 1.

Subsection “Koeppen-Geiger index and bioclimatic variables”: This is an obscure way of reporting differences in local climate and pond longevity. It would make a big difference if you could give me something more concrete than this since the values you report look small and subtle but the implied effects on biology are large. The two do not fit well together.

The Koeppen-Geiger is the most widely used climate classification system. We believe that monthly precipitation over evapotranspiration is a reasonable way to infer water persistence. However, we understand that these measures are no perfect replacement for measured water persistence.

Reviewer #3:

[…]

1) My one comprehensive issue with the manuscript is that it implies too strongly that the dryness of the outlying killifish habitats directly caused those populations to lose diversity and accumulate deleterious mutations. Even if all killifish habitats had identical rainfall and resources, we would still expect the populations at the origin of the range expansion to have higher diversity and less genetic load than smaller, more recently founded, outlying groups. In the absence of any data from control pairs of killifish populations that are geographically distant from each other but have similar levels of rainfall, the paper should clearly state up front that the differences in diversity between populations might not be driven primarily by wetness vs dryness but would be expected to exist whenever one population is more recently founded than another. In addition to addressing this broad issue, it would be helpful for the authors to address the following more specific issues:

We thank the reviewer for requiring clarification on these important aspects.

In a previous study, we have indeed shown a robust phylogenetic correlation between climatic parameters and annual life cycle across different killifish species (Cui et al., 2019). Furthermore, our previous study showed in two independent dry-wet pairs of Nothobranchius species (N. orthonotus and N. rachovii) that populations from dry habitats have higher genetic load. We hypothesize that dry habitats lead to population fragmentation (smaller ponds), hence creating a constellation of small, isolated populations. We now emphasize this point further in Discussion.

We understand the point about founding order, but we believe that this reasoning better applies to populations temporally separated, i.e. with more recently founded populations recently diverging from the ancestral population. However, in our case we sampled populations from dry and wet habitats at the same time. We consider these as divergent populations, i.e. derived from the same ancestral population. Sampling two populations simultaneously prevents us from treating any of the sampled populations as the ancestor, as their divergence time from the ancestral population is identical. Thus, according to this model, genetic differences can occur on either branch. Because in the analysis we were making a relative comparison between two contemporary populations, we believe that founding order won’t be applicable. However, as discussed above, smaller populations in dry regions at the margin of the habitat distribution of the species may experience extreme genetic drift due to population fragmentation, contributing to weakening of purifying selection and accumulation of deleterious gene variants. To note, populations from dry regions are also at higher altitude than populations from wet regions. Hence, gene flow is expected to happen from dry to wet, and not vice versa. This, in turn, may further contribute to genetic diversity depletion in dry areas and maintenance of higher diversity in wet regions.

2) The manuscript notes that FST outlier genes are differentially expressed between fish of different ages but does not comment on how often random control genes likewise vary in expression as a function of age.

We fully agree with this point. We now performed this additional analysis and we couldn’t find an enrichment of significant differentially expressed genes in the FST outliers. We reported the statistics in the text and still report which genes are significant.

3) For the gene with unusually low FST that is located on the sex chromosome, is its outlier status appropriately calibrated with respect to the overall divergence and diversity levels of the sex chromosomes, which can differ substantially from that of the autosomes?

We thank the reviewer for this important point. We now checked the range of the FST values on all chromosomes. We found that the FST values on the sex chromosome (chromosome 3) were neither the lowest nor the highest values of all chromosomes in all comparisons. Here we attach a boxplot of Z-transformed FST values for every chromosome and every comparison we did.

Author response image 1
Download asset Open asset

4) In this discussion of the origin of the killifish sex chromosome, it wasn't clear to me whether the sex determining locus was inherited ancestrally or seems to have acquired a de novo sex determining function.

With our analysis we want to show that the current sex chromosome of the turquoise killifish emerged from two autosomes. Unfortunately, we cannot say whether the sex determining locus was inherited ancestrally. However, since the sex-determining gene occurs in a region corresponding to an ancestral chromosomal translocation, it is likely that it may have acquired its current function simultaneously with this chromosomal rearrangement. However, we cannot yet date the timing of this event.

5) How good is the evidence that the DFE inference method being used is robust to non-equilibrium demographic history?

The program polyDFE can account for demography and polarization error (i.e. knowledge of the ancestral state). We included runs with these options enabled in the analysis of polyDFE. We also simulated SFS (site frequency spectra) using the demography extracted from the PSMC’ results to validate that demography is accounted for and that polyDFE is precise. The simulated data recovered the reflected gamma distribution and confirmed that polyDFE worked for our model. We now present these results as Figure 4—figure supplement 1.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Essential revisions:

1) The authors still do not state how close these populations are to those studied in the QTL study (Kirschner et al., 2012). Couldn't they just provide coordinates, a map, physical distances, or something to help the reader understand how close they are? As I previously stated, if they are not close, then it raises issues of whether these populations would even have the same loci affecting differences in lifespan. The authors provide the new statements below, but it is not clear how to interpret them objectively. Were these fish collected meters from the previous ones? Kilometers? 100 km?

"The populations used in this study are from localities in proximity to those used in the previous QTL study."

"we collected wild populations from natural localities as close as possible to the locations of origin of the strains used for the QTL study."

The revision also adds:

"…or that the populations used in this study and those used for the QTL analysis had a different genetic architecture of lifespan"

This statement could be unnecessary if exact locations show them to be quite close.

We have now provided the exact locations of the populations used in the present and in the QTL work as supplementary Figure 1—figure supplement 1. Additionally, we have specified that these populations belong to the same drainage system. The map shows the physical distance among localities. We have referred to this map in the text and spelled out the relative distance in Km in the figure caption.

2) The explanation of MK α is still inadequate. Some PopGen aficionados will know, but general eLife readers won't know what these values mean. What are high values, low values… etc. What has changed in the sequences to cause these values to change? The authors must really present these statistics in a didactic way.

We have now expanded the section relative to the McDonald-Kreitman statistics to clarify this metrics to a broader audience. We now explain the implications of this metrics and how the plots should be interpreted. I hope this expanded part clarifies better the utility of using this powerful approach.

3) From the new explanations of the inferred DFE of new mutations, it seems this analysis only reflects inferred population demographics, rather than the actual deleteriousness of polymorphisms in the population. Is it correct to say that it uses inferred demographics to just rescale a statistical distribution chosen to represent the DFE?

If so, why not then just present the inferred demographic history? If the DFE analysis isn't analyzing anything about the actual polymorphisms present in the populations, why make that extra inferential jump? Just report the demographics. If this analysis reflects some other quality of the data that translates into actual fitness effects, then it should be stated. If it reflects demographics alone as they would affect a theoretical DFE, then that should be clearly stated.

DFE indeed accounts for both demography and deleteriousness of the variants. Hence, it is not redundant to the inferred demographic history. When inferring the DFE, demography (which is inferred from "neutral"/synonymous sites in the genome) needs to be statistically accounted for, because the pattern of genetic variation observed in the data changes due to past demographic history. However, the inferred DFE itself reflects the deleterious or beneficial newly arising functional mutations (nonsynonymous sites). Hence, DFE takes into consideration the actual polymorphisms present in the populations at neutral and functionally relevant sites, from which it estimates actual fitness effects. By showing how different populations with different demographic history have a different distribution of fitness effects of newly emerging mutations, we do feel that the DFE analysis adds a significant novelty to our results, connecting demography with fitness. We clearly stated in the text that DFE is measured based on polymorphisms “To directly estimate the fitness effect of gene variants associated with each population, we analyzed population-specific genetic polymorphisms to assign mutations as beneficial, neutral or detrimental, and determine the distribution of fitness effect (DFE)42 of new mutations”, hence we are probably unable to better address this point raised by the reviewer.

4) The authors have done a thorough job with revisions and substantially improved their manuscript. The one point from my review that I'd like to revisit is my point that more recently founded populations tend to have smaller effective size and less diversity than older populations regardless of founding order. I don't agree with the reasoning that this only applies to temporally separate populations-there is extensive literature in humans on a "serial bottleneck model" that explains why population diversity is proportional to the distance of a population from the origin of the human range expansion in Africa (see e.g. Ramachandran et al., 2005). This work was based on genetic sampling of human populations at the same time, as was done with the killifish populations here, but founding order is important because more recently founded populations have had less time to recover from the associated bottlenecking founder effect. I don't think it's necessary for the authors to do more work to address this point, but I do think they should point out in the text that some proportion of the population size difference between the wet and dry populations is probably due to the fact that killifish have existed at large population size in the wet habitats for a long time, whereas the dry population probably expanded from a small founder population more recently.

We thank the reviewer for the feedback on our revised manuscript. We also agree with the statement that founding order and population isolation may play a key role in explaining the differences in population size that we observe among populations. Indeed, in the Introduction we mention the importance of population fragmentation in dry populations and we now have added in Discussion a mention to the possible isolation of populations from drier habitats. However, we still know too little about the temporal sequence of population separation. Our favorite hypothesis though is that indeed populations from more coastal areas may be more ancestral. We have now added in Discussion a part referring to how limited gene flow in isolated populations may be also caused by recent founder effect.

https://doi.org/10.7554/eLife.55794.sa2

Article and author information

Author details

  1. David Willemsen

    Max Planck Institute for Biology of Ageing, Cologne, Germany
    Contribution
    Resources, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared

  2. Rongfeng Cui

    Max Planck Institute for Biology of Ageing, Cologne, Germany
    Contribution
    Formal analysis
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8146-6958

  3. Martin Reichard

    1. Czech Academy of Sciences, Institute of Vertebrate Biology, Brno, Czech Republic
    2. Department of Botany and Zoology, Faculty of Science, Masaryk University, Brno, Czech Republic
    Contribution
    Resources
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9306-0074

  4. Dario Riccardo Valenzano

    1. Max Planck Institute for Biology of Ageing, Cologne, Germany
    2. CECAD, University of Cologne, Cologne, Germany
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    dvalenzano@age.mpg.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8761-8289

Funding

Max Planck Society (Valenzano Group core funding)

  • David Willemsen
  • Rongfeng Cui
  • Dario Riccardo Valenzano

Czech Science Foundation (19-01789S)

  • Martin Reichard

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to thank Patience and Edson Gandiwa for their administrative support, Tamuka Nhiwatiwa for helping with logistics and samples handling; Evious Mpofu, Hugo and Elsabe van der Westhuizen and all the rangers of the Gonarezhou National Park for their support in the field. We are thankful to Zimbabwe National Parks for allowing our team to conduct research in the Gonarezhou National Park; Itamar Harel, Matej Polacik and Radim Blazek for hands-on contribution with the field work. We further thank all members of the Valenzano lab for their continuous scientific input and support. The Czech Science Foundation provided financial support to MR for sampling Mozambican populations (19–01789S). This project was funded by the Max Planck Institute for Biology of Ageing, the Max Planck Society and the CECAD at the University of Cologne.

Senior Editor

  1. Jessica K Tyler, Weill Cornell Medicine, United States

Reviewing Editor

  1. Pankaj Kapahi, Buck Institute for Research on Aging, United States

Reviewer

  1. Nathan Clark, Utah, United States

Version history

  1. Received: February 6, 2020
  2. Accepted: August 14, 2020
  3. Version of Record published: September 1, 2020 (version 1)
  4. Version of Record updated: November 20, 2020 (version 2)

Copyright

© 2020, Willemsen et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,937
    Page views
  • 358
    Downloads
  • 20
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. David Willemsen
  2. Rongfeng Cui
  3. Martin Reichard
  4. Dario Riccardo Valenzano
(2020)
Intra-species differences in population size shape life history and genome evolution
eLife 9:e55794.
https://doi.org/10.7554/eLife.55794

Categories and tags

Further reading

    1. Chromosomes and Gene Expression
    2. Cell Biology

    Aging, Geroscience and Longevity: A Special Issue

    Edited by Matt Kaeberlein et al.
    Collection

    eLife is pleased to present a Special Issue to highlight recent advances in the mechanistic understanding of aging and interventions that extend longevity.

    1. Developmental Biology
    2. Evolutionary Biology

    Hierarchical morphogenesis of swallowtail butterfly wing scale nanostructures

    Kwi Shan Seah, Vinodkumar Saranathan
    Research Article

    The study of color patterns in the animal integument is a fundamental question in biology, with many lepidopteran species being exemplary models in this endeavor due to their relative simplicity and elegance. While significant advances have been made in unraveling the cellular and molecular basis of lepidopteran pigmentary coloration, the morphogenesis of wing scale nanostructures involved in structural color production is not well understood. Contemporary research on this topic largely focuses on a few nymphalid model taxa (e.g., Bicyclus, Heliconius), despite an overwhelming diversity in the hierarchical nanostructural organization of lepidopteran wing scales. Here, we present a time-resolved, comparative developmental study of hierarchical scale nanostructures in Parides eurimedes and five other papilionid species. Our results uphold the putative conserved role of F-actin bundles in acting as spacers between developing ridges, as previously documented in several nymphalid species. Interestingly, while ridges are developing in P. eurimedes, plasma membrane manifests irregular mesh-like crossribs characteristic of Papilionidae, which delineate the accretion of cuticle into rows of planar disks in between ridges. Once the ridges have grown, disintegrating F-actin bundles appear to reorganize into a network that supports the invagination of plasma membrane underlying the disks, subsequently forming an extruded honeycomb lattice. Our results uncover a previously undocumented role for F-actin in the morphogenesis of complex wing scale nanostructures, likely specific to Papilionidae.

    1. Evolutionary Biology

    Coevolution of the CDCA7-HELLS ICF-related nucleosome remodeling complex and DNA methyltransferases

    Hironori Funabiki, Isabel E Wassing ... Thomas Carroll
    Research Article

    5-Methylcytosine (5mC) and DNA methyltransferases (DNMTs) are broadly conserved in eukaryotes but are also frequently lost during evolution. The mammalian SNF2 family ATPase HELLS and its plant ortholog DDM1 are critical for maintaining 5mC. Mutations in HELLS, its activator CDCA7, and the de novo DNA methyltransferase DNMT3B, cause immunodeficiency-centromeric instability-facial anomalies (ICF) syndrome, a genetic disorder associated with the loss of DNA methylation. We here examine the coevolution of CDCA7, HELLS and DNMTs. While DNMT3, the maintenance DNA methyltransferase DNMT1, HELLS, and CDCA7 are all highly conserved in vertebrates and green plants, they are frequently co-lost in other evolutionary clades. The presence-absence patterns of these genes are not random; almost all CDCA7 harboring eukaryote species also have HELLS and DNMT1 (or another maintenance methyltransferase, DNMT5). Coevolution of presence-absence patterns (CoPAP) analysis in Ecdysozoa further indicates coevolutionary linkages among CDCA7, HELLS, DNMT1 and its activator UHRF1. We hypothesize that CDCA7 becomes dispensable in species that lost HELLS or DNA methylation, and/or the loss of CDCA7 triggers the replacement of DNA methylation by other chromatin regulation mechanisms. Our study suggests that a unique specialized role of CDCA7 in HELLS-dependent DNA methylation maintenance is broadly inherited from the last eukaryotic common ancestor.