An oligogenic architecture underlying ecological and reproductive divergence in sympatric populations

Abstract

The evolutionary trajectories and genetic architectures underlying ecological divergence with gene flow are poorly understood. Sympatric timing types of the intertidal insect Clunio marinus (Diptera) from Roscoff (France) differ in lunar reproductive timing. One type reproduces at full moon, the other at new moon, controlled by a circalunar clock of yet unknown molecular nature. Lunar reproductive timing is a magic trait for a sympatric speciation process, as it is both ecologically relevant and entails assortative mating. Here, we show that the difference in reproductive timing is controlled by at least four quantitative trait loci (QTL) on three different chromosomes. They are partly associated with complex inversions, but differentiation of the inversion haplotypes cannot explain the different phenotypes. The most differentiated locus in the entire genome, with QTL support, is the period locus, implying that this gene could not only be involved in circadian timing but also in lunar timing. Our data indicate that magic traits can be based on an oligogenic architecture and can be maintained by selection on several unlinked loci.

Keywords: Clunio marinus; circalunar clock; ecological speciation; genetics; genomics; magic trait; sympatric speciation.

Figure 1.. Population structure of the FM and NM types in Roscoff based on 721,000 genetic variants.

(A, B) Principal component analysis (PCA; A) and admixture analysis (B) identify one migrant in time (blue; pure NM genotype caught at full moon), many potential F1 hybrids (yellow) and four individuals in the FM strain that appear genetically distict from all other samples (red). (C) Global genetic differentiation is limited, but there is a block of strong differentiation on chromosome 1.

Figure 2.. A complex inversion system on chromosome 1.

(A) Chromosome 1 harbors a block of markedly elevated genetic differentiation. (B) This genomic block coincides with two windows of elevated long-range LD in the FM and NM strains. (C) Structural variant (SV) calling from long read sequencing data supports that the larger LD block is due to an inversion. (D–G) Principal component analysis (PCA) for chromosomal sub-windows corresponding to suggested inversions separates the individuals into clusters corresponding to standard haplotype homozygotes (SS), inversion homozygotes (I₁I₁, I₂I₂) and inversion heterozygotes (SI₁, SI₂, SI₃, I₁I₂). The four individuals which were already found distinct in whole genome analysis (red question mark) cannot be assessed. (H–K) Observed heterozygosity is clearly elevated in inversion heterozygotes, underpinning substantial genetic differentiation between the inversions. (L) A schematic overview of the sequence of inversion events and the resulting haplotypes. The frequencies of the haplotypes differ markedly between the FM and NM types.

Figure 2—figure supplement 1.. Long-range linkage disequilibrium along chromosome 1.

Blocks of elevated genetic differentiation (A) correspond to blocks of elevated long-range LD in the FM type (B) and the NM type (C). Elevated long-range LD suggests that the inversion is polymorphic in the respective population.

Figure 2—figure supplement 2.. Local admixture plots for the genomic windows corresponding to different inversions.

Gray bars and black bars represent the two classes of inversion homozygotes. Bars which are half black, half gray represent inversion heterozygotes. Support for In(1c) is not fully congruent with support from heterozygosity and PCA analyses. As for In(1c) only one group of homozygotes was detected in PCA, the fully gray bars may indicate inversion heterozygotes. Similarly, for In(2L) and In(3L) only one class of homozygotes was detected in PCA, suggesting that here the fully black bars represent heterozygotes.

Figure 2—figure supplement 3.. Genetic differentiation between homozygotes of the standard haplotype (SS) and inversion haplotypes (I₁I₁, I₂I₂).

Figure 2—figure supplement 4.. Windowed PCA along chromosome 1, part 1 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 5.. Windowed PCA along chromosome 1, part 2 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 6.. Windowed PCA along chromosome 1, part 3 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 7.. Windowed PCA along chromosome 1, part 4 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 8.. Windowed PCA along chromosome 1, part 5 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 9.. Windowed PCA along chromosome 1, part 6 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 10.. Windowed PCA along chromosome 1, part 7 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 11.. Windowed PCA along chromosome 1, part 8 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 12.. Windowed PCA along chromosome 1, part 9 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 13.. Windowed PCA along chromosome 1, part 10 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 14.. Windowed PCA along chromosome 1, part 11 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 15.. Windowed PCA along chromosome 1, part 12 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 16.. Windowed PCA along chromosome 1, part 13 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 17.. Windowed PCA along chromosome 1, part 14 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 18.. Windowed PCA along chromosome 1, part 15 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 19.. Windowed PCA along chromosome 1, part 16 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 20.. Windowed PCA along chromosome 1, part 17 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 21.. Windowed PCA along chromosome 1, part 18 (500kb sliding windows with 100 kb steps).

Figure 2—figure supplement 22.. Windowed analysis for observed heterozygosity (A) and Admixture (B) along chromosome 1.

(A) A number of individuals (red arrows) show elevated heterozygosity beyond the end of In(1a). In three of them (individuals 38, 41, and 43), the ends of that region coincide (red lines). These individuals are considered SI₃ heterozygotes, carrying In(1c). In the other two individuals, the region is shorter or longer, possibly pointing to yet other inversions. (B) Admixture scores from windowed analysis support the same individuals to be genetically distinct in the right arm of In(1c).

Figure 3.. Inversions in chromosomes 2 and 3.

(A) Chromosome arm 2L harbors a block of mildly elevated genetic differentiation. (B) Several blocks of long-range linkage disequilibrium (LD) point to additional inversions in the FM and NM strains, one on each chromosome arm. (C) Structural variant (SV) calling from long read sequencing data supports the inversion on chromosome arm 2R. (D–G) Principal component analysis (PCA) for chromosomal sub-windows corresponding to suggested inversions separates the individuals into clusters corresponding to standard haplotype homozygotes (SS), inversion homozygotes (II) and inversion heterozygotes (SI). The four individuals which were already found distinct in whole genome analysis (red question mark) cannot be assessed. (H–K) Observed heterozygosity is clearly elevated in inversion heterozygotes, underpinning substantial genetic differentiation between the inversions. (L–O) The frequencies of the haplotypes do not differ much between the FM and NM types.

Figure 3—figure supplement 1.. Long-range linkage disequilibrium along chromosomes 2 and 3.

The block of slightly elevated genetic differentiation in chromosome arm 2L (A) correspond to a block of elevated long-range LD in the NM type (C), but not the full moon type (B). Blocks of elevated LD are also found in all other chromosome arms. Elevated long-range LD suggets that the inversion is polymorphic in the repective population, congruent with the genotyping results presented in Figure 2 (the higher the allele frequency of the inverted haplotype, the stronger is LD).

Figure 4.. Lunar emergence time is heritable.

In a cross between FM type (A) and NM type (B), the F1 hybrids emerges intermediate between the parents (C). In the F2 (D), the phenotypes spread out again, but do not segregate completely, suggesting that more than one locus controls lunar emergence time. The peaks around day 20 (gray shading) are direct responses to the artificial moonlight treatment (see Kaiser et al., 2021). They do not occur in the field and were not considered in our analyses. Yellow shading indicates the days with artificial moonlight treatment in the laboratory cultures.

Figure 5.. Quantitative trait loci (QTL) and genomic regions underlying FM vs NM emergence.

(A) Multiple QTL Mapping (MQM) identified four significant QTL controlling lunar emergence time. Confidence intervals were determined in MQM and Composite Interval Mapping (CIM). The marker order supports In(1a) and the nested In(1b) on chromosome 1, as well as In(2L) on chromosome 2. (B) QTL intervals were colour-coded, transferred to the reference sequence and assessed for loci which are strongly differentiated (F_ST) between the FM and NM types. Inversions are represented by blue bars below the plot. The QTL in In(1a) overlaps with the most differentiated locus in the genome (red asterisk) – the period locus. Other markedly differentiated loci within the QTL intervals are the plum/cask locus (red circle) and the stat1 locus (red cross). (C) Genetic differentiation is strongest in the first intron and the intergenic region just upstream of the period gene. (D) The period locus was re-genotyped based on an insertion-deletion mutation in the first intron (green arrows in (C)). The long (L) and short (S) period alleles are only loosely associated with inversion haplotypes.

Figure 5—figure supplement 1.. Scan for interacting QTL.

All chromosmes were scanned for interacting QTLs with the scantwo function of R/qtl. (A) For each possible QTL combination three models are tested, allowing for either additve and epistatic effects (full model; Mf), only additive effects (Ma) or assuming that one of the two QTL does not exist (M1). (B) The LOD score for the full model (Mf) and the additive-only model (Ma) are similar across the genome. (C) Congruently, the LOD scores are much lower for the difference between Mf and Ma (upper triangle). The additive model Ma does much better in explaining the data than assuming only one QTL exists (Ma-M1; lower triangle). (D,E) Statistical tests for the detected QTL combinations. Only in one case the full model is significantly better than the additive-only model (see p(Mf-Ma)), suggesting that epistatic effects are limited to an interaction between the middle of chromosome 1 and the end of chromosome 2 (D). The additive model is significantly better than assuming there is only one QTL (E).

Figure 5—figure supplement 2.. LOD profiles from Composite Interval Mapping (CIM).

(A) Backward selection models. (B) Forward selection models. CM = size of the exclusion window in centimorgan; COV = number of covariates. CIM generally identifies the same four QTL as MQM, but sometimes identifies several adjacent QTL in the same genomic region.

Figure 5—figure supplement 3.. Inversion genotypes of the parents of the FMxNM cross.

Inversion genotypes were determined by local PCA for the inversions overlapping with detected QTL. For the FM type mother the inversion genotype of In(2R) is not clear. Inversion genotypes suggest partial recombination suppression in the F1 for all inversions.

Figure 5—figure supplement 4.. Genotyping of an insertion-deletion in the period locus.

The NM strain is dominated by the long allele (L), the FM strain by the short allele (S).

Figure 5—figure supplement 5.. Gene models and genetic differentiation at the plum locus.

Dots are individual SNPs, the line is F_ST in 1 kb sliding windows with 200 bp steps. In the gene models, gray boxes are untranslated regions (UTRs), black boxes are coding sequence (CDS).

Figure 5—figure supplement 6.. Gene models and genetic differentiation at the stat1 locus.

Dots are individual SNPs, the line is F_ST in 1 kb sliding windows with 200 bp steps. In the gene models, gray boxes are untranslated regions (UTRs), black boxes are coding sequence (CDS).