Genomic divergence landscape in recurrently hybridizing Chironomus sister taxa suggests stable steady state between mutual gene flow and isolation

Abstract Divergence is mostly viewed as a progressive process often initiated by selection targeting individual loci, ultimately resulting in ever increasing genomic isolation due to linkage. However, recent studies show that this process may stall at intermediate stable equilibrium states without achieving complete genomic isolation. We tested the extent of genomic isolation between two recurrently hybridizing nonbiting midge sister taxa, Chironomus riparius and Chironomus piger, by analyzing the divergence landscape. Using a principal component‐based method, we estimated that only about 28.44% of the genomes were mutually isolated, whereas the rest was still exchanged. The divergence landscape was fragmented into isolated regions of on average 30 kb, distributed throughout the genome. Selection and divergence time strongly influenced lengths of isolated regions, whereas local recombination rate only had minor impact. Comparison of divergence time distributions obtained from several coalescence‐simulated divergence scenarios with the observed divergence time estimates in an approximate Bayesian computation framework favored a short and concluded divergence event in the past. Most divergence happened during a short time span about 4.5 million generations ago, followed by a stable equilibrium between mutual gene flow through ongoing hybridization for the larger part of the genome and isolation in some regions due to rapid purifying selection of introgression, supported by high effective population sizes and recombination rates.

the divergence process proceeds in time, how widespread the phenomenon is, and whether it always and inevitably leads to complete isolation. Comparing the genomes of individuals of two regularly hybridizing sister taxa of nonbiting midges, we could show that they diverged during a short period millions of generations ago. Their divergence process apparently ceased before the entire genome was mutually isolated. The taxa remain distinct since, even though they share most of their genome. Our findings thus extend our view of the nature If selection acts simultaneously in disparate directions, favoring different phenotypes, divergence processes are initiated (Rundle and Nosil 2005). However, such divergence processes may not need to affect the genome as a whole but rather lead to heterogeneous genomic differentiation (Martin et al. 2013;Feulner et al. 2015;Duranton et al. 2018;Ravinet et al. 2018). When selection targets individual loci underlying the favorable phenotype, only these as well as linked loci will diverge (Via and West 2008). These diverging genomic regions under selection are shielded from gene flow by locally reduced recombination rates (Via 2012), whereas the predominant part of the genome remains freely exchangeable among emerging taxa (Wu 2001;Feder et al. 2012b). This process is called divergence hitchhiking (DH; Via and West 2008) and results in heterogeneous genomic differentiation patterns with small genomic regions of elevated differentiation referred to as "islands of divergence" (Turner et al. 2005). Such patterns of genomic divergence are already well documented across a wide range of taxa: Studies on Anopheles gambiae (Turner et al. 2005) as well as Ficedula species (Ellegren et al. 2012) revealed "islands of divergence" often located near telomeres and centromeres, regions with reduced recombination rates (Chen et al. 2008). In Heliconius species, divergence peaks were found in presumed regulatory regions associated with wing patterns that are commonly involved in mate recognition and therefore under sexual selection (Nadeau et al. 2012). At the same time, Heliconius species experience frequency-dependent selection on local patterns due to their role in predator avoidance (Nadeau et al. 2012). Taking the temporal course of divergence into account, a study on Dicentrarchus labrax found several regions responsible for initial reproductive isolation (RI) evolved during allopatry in the two lineages that were eroded by gene flow during secondary contact (Duranton et al. 2018).
Genomic divergence is thus a progressive process, initiated by selection targeting local genomic regions and expanding from there over time, eventually leading to speciation with gene flow (Feder et al. 2012b;Via 2012). Feder et al. (2012b) recognized four phases in this process. (i) During the first phase, direct selection targets (independent) loci inducing divergence, whereas the rest of the genome is still exchanged by gene flow. (ii) The second phase is characterized by DH elevating divergence in the proximity of loci under selection. (iii) In the third phase, DH is either already widespread throughout the genome or the selection is strong enough to decrease gene flow on a genome-wide scale. This leads to the fixation of further divergent loci resulting in "continents of divergence." (iv) This so-called genome hitchhik-ing (GH) eventually leads to a fourth phase in which gene flow increasingly ceases throughout the genome.
However, this sequence of events does not appear to be an inevitably progressive process. Several empirical studies (Schwenk et al. 2000;Faure et al. 2009;Canestrelli et al. 2017) as well as theoretical models (Yeaman and Whitlock 2011;Flaxman et al. 2014;Rafajlović et al. 2016) suggest the possibility of intermediate constant equilibrium states where certain parts of the genome remain diverged ("islands" or "continents of divergence"), whereas others are freely exchanged among closely related species without ever reaching complete genomic isolation (GI). Theoretical prerequisites for such a migration-selectiondrift equilibrium, effectively freezing the divergence process in some phase of the Feder et al. model (2012b), are beneficial mutations locally connected through linkage and either reduced recombination or increased selection in genomic regions differing in divergent taxa (Yeaman and Whitlock 2011;Rafajlović et al. 2016). However, it remains unclear how frequent such equilibria among hybridizing taxa occur.
We investigated the state of divergence and the possibility of a stable migration-selection-drift equilibrium in two cryptic, nonbiting midges sister taxa Chironomus riparius (Meigen, 1803) and Chironomus piger (Strenzke, 1956). Both of these multivoltine taxa are found in freshwater habitats such as small streams, ditches, and puddles (Pfenninger and Nowak 2008) throughout the temperate Holarctic (Strenzke 1957;Oppold 2017) and spend most of their life in one of four aquatic larval stages before they pupate, emerge, and mate in swarms (Armitage et al. 1995). The two taxa often co-occur throughout their largely overlapping distribution ranges with no known differential centers of distribution (Kiknadze et al. 1991;Gunderina and Salina 2003;Petrova et al. 2014;Pedrosa et al. 2017). Usually, one taxon prevails due to differential adaptation to ecological niches (Pfenninger and Nowak 2008). Chironomus piger is known, for example, to tolerate higher nitrite pollution and salt concentrations (Pfenninger and Nowak 2008), which have negative impact on larval development in both taxa (Haas and Strenzke 1957;Neumann et al. 2001). Thus, C. piger usually predominates in lentic waterbodies near agricultural landscapes (Foucault et al. 2019), whereas C. riparius prevails in slow-flowing streams.
Chironomus riparius and C. piger differ remarkably in genome size: the C. riparius genome is around 30% larger compared to the sister taxon. This increased genome size of C. riparius can be explained in large part by the expansion of a certain transposable element, the so-called Cla-element. In addition, this element has been suggested to be involved in the initiation of the speciation process (Schmidt 1984), which was estimated to have occurred about 1.3-1.8 million years ago (Schmidt et al. 2013). Karyotype comparisons with more distantly related Chironomus species further support the assumption that the smaller genome of C. piger is the ancestral state (Keyl 1965). The suspicion that Cla-elements are drivers of divergence is strengthened by strong selection against heterozygous population-specific insertions in interpopulation crosses of C. riparius . Despite this obvious genomic difference, hybridization between the two taxa has been observed both in the lab (Foucault et al. 2019) and in nature (Petrova et al. 2014;Pedrosa et al. 2017;Foucault et al. 2019). However, substantial pre-and postzygotic RI barriers, namely, optical swarming markers (Miehlbradt and Neumann 1976) and hybrid dysgenesis syndromes (Hägele 1984(Hägele , 1999Armitage et al. 1995), lead to generally lower fitness in all (back)crossing directions (Foucault et al. 2019). This results in strong selection against introgression reflected in the low abundance of recently admixed individuals in the field (Pedrosa et al. 2017;Foucault et al. 2019). Given the estimated divergence time of about 1.5 million years ago (Schmidt et al. 2013), which corresponds to at least 9 million generations (Oppold et al. 2016), it is remarkable that hybridization is still ongoing. Strong divergent selection, long-standing divergence, and persistent hybridization make the sister taxa C. riparius and C. piger a promising system to investigate the extent of GI as well as the state of the divergence process, in particular the possibility of a stable selectionmigration-drift equilibrium.

SAMPLING, SEQUENCING, AND DATA PROCESSING
We whole genome sequenced 36 individual Chironomus specimen from five different sites across Europe ( Fig. S1; Table S1). Sequencing was performed on an Illumina HiSeq4000 platform, using the KAPA library preparation kit, resulting in around 25× coverage 150 bp paired-end reads per sample.
After quality checking, trimming, and mapping to the latest C. riparius reference genome (Schmidt et al. 2020), single nucleotide polymorphisms (SNPs) were called with GATK (De-Pristo et al. 2011). The filtered VCF file contained 5,496,457 variable positions. For further details on tools, options, parameters, and commands, see Supporting Information 1.1.

ADMIXTURE ANALYSES
Using bcftools version 1.9 (Li 2011), we filtered out strongly linked loci with r 2 > 0.8 in 100-kb windows to prune the dataset for SNPs in linkage disequilibrium (LD; Alexander et al. 2009;Baran et al. 2013). This retained 900,345 loci. ADMIXTURE version 1.3.0 (Alexander et al. 2009) was used to infer ancestry proportions based on maximum likelihood estimations using K = 2 with default settings. In addition, we used the R (R Core Team 2017) package FACTOEXTRA (Kassambara and Mundt 2017) to calculate a principal component analysis (PCA) on unlinked SNPs successfully called in all individuals.

REGIONS
We developed a PCA-based approach for identifying genomic regions that are mutually isolated between two hybridizing taxa, which requires no a priori known admixture-free reference individuals for the respective species. PCA is not based on any particular population genetic model (Zheng and Weir 2016) and is therefore widely used to infer major trends in individual-based SNP data.
PCA was performed on all SNPs successfully genotyped in all individuals using the R (R Core Team 2017) package FAC-TOEXTRA (Kassambara and Mundt 2017). The first PC describes the largest proportion of variance in the data, meaning it splits the most divergent units, that is, the described taxa as identified by DNA barcoding (Table S1). On this axis, we identified SNPs whose cumulative sum of factor scores increased linearly (species separating SNPs [sepSNPs]; Fig. S2). The remaining SNPs contributed only negligibly to the variance among species (Fig. S2) and were flagged as residual SNPs.
Isolated genomic areas were defined as nonoverlapping windows containing more than statistically expected sepSNPs (χ 2test with Benjamini-Hochberg FDR correction for multiple testing). To retain sufficient statistical power, we chose a window size containing at least 60 SNPs in more than 95% of the windows (i.e., 10 kb). Windows smaller than the chosen size arose from windows at scaffold ends.
Additionally, we estimated the fixation index (F ST ; Weir and Cockerham 1984) in the same 10-kb windows using VCFtools version 0.1.15 (Danecek et al. 2011). We calculated the point biserial correlation between F ST values per window and the classification as isolated or nonisolated based on our PCA approach.

ISOLATION
We estimated mean LD for each 10-kb window using the maximum likelihood method on Pool-Seq data of Feder et al. (2012a) on a published pool from one of the populations sampled here (ENA project ERP115516, sample ERS4040036). This method exploits the haplotype information contained in (paired) reads and yields accurate estimates if LD decay to background levels is shorter than the read length (Feder et al. 2012a), which was 250 bp here (see Fig. S3). The extension to the entire species assumed that the recombination landscape is highly correlated within species (Samuk et al. 2020).
We conducted gene ontology (GO) term enrichment analyses on the categories Biological Process (BP) and Cellular Component (CC) on genes in isolated and nonisolated regions as well as on positively selected genes (see below). The analyses on 7951 GOannotated genes were carried out using the R package TOPGO (Alexa and Rahnenführer 2018). Only terms with more than five annotated genes were considered. See Supporting Information 1.4 for details.

HISTORIC DIVERGENCE DYNAMICS
According to the DH concept, the size of isolated genomic regions should increase with time until finally the entire genome is isolated (Feder et al. 2012b). This yields three testable predictions: i) As divergence islands tend to increase in size over time, their mean length should be larger than islands resulting from randomly distributing the same number of identified isolated windows over the genome (see Supporting Information 1.5). ii) In most models of speciation with gene flow, GI inevitably increases over time until it is complete or reaches an equilibrium state. Isolated regions should therefore increase both in number and size over time and divergence time estimates of isolated region should extend to the present (or when full GI was reached) or until a stable steady-state equilibrium is reached. To infer the isolation dynamics, we used the temporal information contained in the divergence time estimates for each isolated 10-kb window (Supporting Information 1.5). We then simulated five plausible divergence scenarios with SIMCOAL (Excoffier 2000; see supplementary code in Supporting Information 1) and inferred the best-fit scenario to the empirical divergence time distribution of isolated windows using an approximate Bayesian computation (ABC) inference framework (see Supporting Information 1.5 for details): a) Continuous divergence: Divergence is continuously acting over the entire time span with the same number of windows diverging per temporal bin. b) Early eruption, fading out: After the initial, early divergence of a few windows, the number of diverging windows increased very fast to eventually slow down but continue until today. c) Escalating divergence: Divergence started early with few windows and continuously increases over time, both in number and size of divergent windows. d) Episodic divergence: Divergence started with a few windows, increased fast, peaks at the mean observed divergence time, and completely ceased again. e) Instantaneous divergence: All divergence happened at once, at the mean observed divergence time.
The respective simulated frequency distributions for each model are depicted in Figure 1A. Please note that the simulations focused on currently diverged regions. All scenarios therefore include the possibility of secondary contact after an allopatric phase of divergence.
iii) If divergence islands grow with time, their size should be roughly proportional to their estimated divergence time. Additionally, the variance in size not explained by divergence time is expected to result from a reduction of gene flow by selection and physical linkage between pairs of SNPs (Wu 2001). To test this hypothesis, we joined adjacent isolated windows to isolated regions, calculated their length, and correlated it to the respective mean divergence time, strength of selection, and physical linkage (see Supporting Information 1.5 for details).

PROCESSES DRIVING DIVERGENCE
We calculated Tajima's D (TD) for each isolation window of each taxon separately with MEGA X (Kumar et al. 2018) and compared the resulting distribution with a random sample of 1000 nonisolated windows in each taxon. Isolated windows with TD values falling beyond the upper or lower 95% confidence limits of the nonisolated windows distribution were considered as driven by nonneutral processes (Feulner et al. 2015). We considered outlier windows with significantly negative TD in both taxa as under divergent selection, while occurring in only a single taxon as positively selected in the respective taxon (Pfenninger et al. 2015). We also evaluated the influence of the Cla elements on divergence. It is known that these repeat elements are enriched in heterochromatic regions (Schmidt et al 2020). As these are the regions that could not be resolved in the current genome assembly, scaffold ends tend to neighbor the clusters of Cla elements. To test for potential influence of Cla elements on divergence, we analyzed whether isolated regions were more often closer to scaffold ends than expected by chance.
To relate adaptive protein evolution to the divergence process, we extracted the coding sequence of all genes in the isolated regions. Using PopGenome (Pfeifer et al. 2014), we calculated the approximate version of the McDonald-Kreitman test (MKT). Significance was assessed with a Fisher's exact test. For all significant genes, we calculated the neutrality index and the proportion of divergent SNPs fixed by positive selection (α; Smith and Eyre-Walker 2002). For comparative reasons only, we did the same for genes in nonisolated regions, even though the MKT test assumption of divergence was by definition not met (Supporting Information 1.6).

GENOME-WIDE (RECENT) ANCESTRY INFERENCE
Both ADMIXTURE (Alexander et al. 2009) and a PCA approach confirmed the DNA barcoding-based taxa assignment to C. riparius and C. piger and did not find indications of recent admixture (Table S1; Fig. S5). ADMIXTURE estimated admixture proportions below 1% across all samples. The PCA clustered the individual samples according to the a priori assignment and did not find any intermediate individuals (Fig. S6).

CURRENT STATE OF MUTUAL GI
The PCA for identification of highly differentiated genomic regions was calculated based on the entire genotype dataset containing 3,270,391 SNPs successfully genotyped for all samples ( Fig. 2A). PC1 explained 53% of variance and split C. riparius from C. piger. PC2, explaining 4.4% of variance, reflected intraspecific divergence of C. riparius ( Fig. 2A). Cumulative locus factor scores of the 700,000 most contributing positions explained 90.88% of the variance of PC1 (sepSNPs). Using FDR-corrected χ 2 -tests to identify 10-kb windows with significantly accumulated sepSNPs, 4917 out of 18,203 windows were identified as isolated (Supporting Information 2). These isolated windows were spread over the entire genome (Fig. 2B) and characterized by increased sequence divergence (D xy ) compared to nonisolated ones (Figs. 2B and 2C). Joining immediately adjacent isolated windows reduced this number to 1626 isolated regions ranging from 1166 bp to 350 kb length (Fig. 2D). Isolated regions were found on 488 of 752 scaffolds of the reference genome.
Our PCA approach estimated thus 28.44% of the genome to be mutually isolated. F ST values were positively correlated with isolated regions identified by the PCA approach (Point-biserial correlation, r pb = 0.74, P < 0.001).
Isolated regions were also significantly further away from scaffold ends and thus from large Cla clusters than is expected by chance (empirical distance: 140 ± 138 kb [mean ± SD], simulated distance: 59 ± 1.3 kb [mean ± SD], P < 0.001).

CHARACTERIZATION OF ISOLATED WINDOWS
The inferred isolated windows differed significantly from the nonisolated ones in all parameters estimated across the genome. The largest differences were found in F ST , D xy, and genetic diversity (θ), all having large effect sizes (Figs. 3B, S9B, and S9C). LD as a proxy for recombination rate only had a small effect (Supporting Information 1.5; Fig. S9A).
Observed isolated regions were with on average 30.2 ± 32.8 kb significantly larger than expected by chance (simulated mean 13.3 ± 0.1 kb, P < 0.0001, Supporting Information 1.5). The isolated regions contained 6723 out of 13,449 annotated genes, which is almost twice the expected number, assuming an even distribution of genes across the genome (χ 2 = 2957, P < 0.0001). The genes in isolated regions were significantly enriched for GO terms associated with protein synthesis, ranging from transcription to protein modification in the category Biological Function (Table 1). In the category Cellular Component, mainly GO terms relating to complex cell structures were found to be enriched, whose function requires close molecular interaction of many constituting protein or RNA components, such as ribosome, mitochondrion and nucleus (Table 1).

TEMPORAL DIVERGENCE DYNAMICS
The mean divergence time for isolated 10-kb windows was estimated to 4,598,663 ± 1,588,340 generations ago, corresponding to 3.06 ± 1.05 N e generations, using a molecular clock approach. The frequency distribution of estimated divergence times was comparatively centered and deviated significantly from normality (Shapiro-Wilk's W = 0.99, P < 0.0001; Fig. 1B). The ABC inference framework identified the episodic model (Fig. 1B) as producing most similar results compared to the empirical data, indicated by a Bayes factor of 1.

NEGATIVE TD THROUGHOUT THE GENOME
TD was skewed toward negative values in both isolated and nonisolated windows. In nonisolated windows, values ranged between 0.064 and −2.611 in C. riparius and between 0.47 and −3.48 in C. piger (Fig. 3B). The 5th and 95th percentiles were −1.39, −2.36, and −0.79, −2.32, respectively. A total of 365 Gene ontology (GO) term enrichment analysis in isolated parts of the genome. Significantly (P < 0.05) overrepresented GO terms in the categories biological process (BP) and cellular component (CC) for isolated regions. The column "Annotated" shows the number of genes annotated with the respective GO term, the column "Observed" the number of these genes in the isolated region, and the column "Expected" the number expected if these genes were equally distributed over the entire genome. Probability values for the number of observed genes being within random expectations are based on the weight Fisher algorithm.  isolated windows in C. riparius and 464 in C. piger were below the lower threshold. Of these isolated windows, 68 occurred in both taxa, which was significantly more than expected by chance (randomization test [Pfenninger et al. 2015]: mean expected = 31.5, P < 0.001). The upper expectation threshold for TD was exceeded by 297 isolated windows in C. riparius and 396 in C. piger. The overlap of 134 windows was much higher than expected by chance (randomization test [Pfenninger et al. 2015]: mean expected = 22, P < 0.001). A list of genes in the windows considered to be driven by nonneutral processes (see methods above) can be found in Supporting Information 3.

ISOLATED REGIONS
For the 6723 genes in the isolated regions, McDonald-Kreitman tests showed a significant accumulation of either synonymous or nonsynonymous divergent SNPs for 451 genes. Neutrality index indicated that 421 were under negative selection, whereas 30 showed signs of positive selection. No significantly enriched GO terms were found among the putatively positive selected genes. From the 392 fixed divergent nonsynonymous SNPs in positively selected genes, a proportion of 0.89 was estimated to have become fixed by this process. The list of positively selected genes in isolated regions can be found in Supporting Information 3.

Discussion
In this study, we focused on evolution and dynamics of the divergence landscape in the face of gene flow between two sister species of nonbiting midges, C. riparius and C. piger. We used a novel PCA approach not requiring admixture-free reference individuals to unravel the heterogeneous genomic differentiation across the genome. We were able to distinguish regions of significantly reduced gene flow ("isolated") from apparently freely exchangeable ones ("nonisolated"). We analyzed the driving forces behind formation of islands of divergence, inferred the temporal divergence dynamics, and investigated whether the divergence processes are still ongoing or if a stable steady-state selectionmigration-drift equilibrium was reached.

METHODOLOGICAL CONSIDERATIONS
Sampling five different C. riparius populations and comparing them to a single C. piger population to derive taxa-wide genomic differentiation patterns might have been problematic; however, this approach proved to be sufficient. The equidistance between the respective C. piger and all C. riparius populations in the PCA ( Fig. 2A) showed that the specific differences are rather categorical and less dependent on geographical proximity, local adaptation, or differential introgression. The homogeneous, specific differences may result from homogenizing intraspecific gene flow (Kumar et al. 2017), leading to similar isolation patterns across the taxa's range, even in populations never in direct contact with the sister taxon. High intraspecific gene flow among C. riparius populations (Waldvogel et al. 2018) coupled with low hybridization rates between the taxa (Pedrosa et al. 2017;Foucault et al. 2019) supports this explanation. We could reliably identify the heterogeneous divergence landscape among the recurrently hybridizing taxa. To identify mutually isolated and therefore highly differentiated genomic regions between hybridizing taxa, most studies use summary statistics such as F ST or similar (Nosil et al. 2009;Wolf and Ellegren 2017). However, these statistics are criticized for losing explanatory power due to multiple factors such as selection, recombination, and mutation rate, which all vary across genomes (Hodgkinson and Eyre-Walker 2011;Ravinet et al. 2017). This is particularly true when using whole genome data at late stages of divergence (Meirmans and Hedrick 2011) as these effects accumulate with time and number of loci considered.
Alternatively, ancestry of genome blocks and thus regions of introgression and isolation can be inferred on the basis of (putatively) nonadmixed reference individuals (Price et al. 2009;Schaefer et al. 2017). However, these individuals probably do not exist in many cases for a variety of reasons, for example, recurrent hybridization since initial divergence as in the current case. Moreover, even in populations never in direct contact with the sister lineage, genomic integrity may be compromised by intraspecific gene flow (Kumar et al. 2017). Therefore, the basic assumption of reference-based admixture analysis is likely unwittingly violated in many instances.
By maximizing the variance among correlated multivariate variables, a SNP-based PCA naturally finds on its first axis the SNPs contributing most to the differences between the highest order taxa. If these SNPs are spatially clustered and not homogeneously distributed across the genome, it is reasonable to assume that gene flow rates differ locally along the genome among the inferred taxa. Due to the time necessary for a given region to accumulate a certain number of divergent mutations, however, the method has natural limits. Simulations showed that for the most recent 2.5 million (1.7 N e ) generations, likely not all isolated windows were identified as such, thus underestimating the number and size of isolated regions for this period (false negative rate, Fig. S7A). On the other hand, windows flagged as isolated had a very low false positive rate even for relatively recently isolated windows (>200,000 generations; Fig. S7B).
As expected, identified isolated and nonisolated windows differed significantly in net divergence (Figs. 3A and S9C). Additionally, we found a strong positive correlation between the identified isolated windows and F ST values. This suggests that an F ST outlier approach (Nosil et al. 2009;Wolf and Ellegren 2017) would have resulted in a comparable inference of the heterogeneous divergence landscape. The remaining differences could be explained by F ST evaluating relative allele frequency differences, whereas the PCA approach assesses reduction in gene flow as well as divergent mutation accumulation in certain genomic regions.
As in all window approaches, the spatial resolution for delimiting isolated genomic regions is influenced by the chosen window size. It was necessary to find a balance between statistical power, which increases with increasing window size, and a biologically relevant resolution, which depends on the local recombination rate (Teo et al. 2009). Using LD as a proxy for the recombination rate showed that average haplotype length in C. ri-parius was significantly smaller than the chosen 10-kb windows (Schmidt et al. 2020). Consequently, some windows flagged as either isolated or nonisolated may actually represent a patchwork of both. This renders the estimation of all parameters calculated based on windows less accurate. However, choosing a much smaller window size would have made several parameter estimations impossible, due to the lack of informative sites, and increased variance of estimates. Choosing a suitable window size therefore strongly depends on the system under investigation and cannot be generalized here.

EXTENT OF GI
Mutually isolated regions between C. riparius and C. piger were dispersed across almost two thirds of all scaffolds, but accounted for only 28% of the total assembled genome. Interestingly, however, the isolated regions contained half of all annotated genes. Conversely, that means that 72% of the genome, containing the remaining 50% of genes, appears to be freely exchangeable between the taxa. Therefore, the nonisolated parts, comprising the majority of the genome, resemble a shared gene pool, which is reflected in higher N e (represented by θ) compared to isolated regions (Fig. S9B). Species identity hence seems to be based on only half of the genes.
Acknowledging the restrictions and limitations of GO annotations and their interpretations (Sangar et al. 2007), we found conspicuously many genes in isolated regions associated with processes or structures that require close molecular interactions with protein or RNA coding genes. These included genes associated with mitochondrial functions, protein synthesis at ribosomes or chromosome organization, and replication in the nucleus. Selection for efficiency in processes requiring close molecular or functional interaction of several genomic regions (e.g., oxidative phosphorylation in eukaryotes; Zhang and Broughton 2013) could be the reason. Disruption of these interactions has been shown to reduce hybrid fitness (Ellison et al. 2008;Zhang and Broughton 2013) and thereby foster divergence. It is therefore not surprising that strong postzygotic RI barriers exist (Hägele 1984(Hägele , 1999Armitage et al. 1995), leading to reduced fitness in all (back)crossing directions (Foucault et al. 2019). In particular, the divergent genes associated with chromosome organization and replication could account for observed chromosome aberrations in hybrids (Hägele 1984).
Other prime candidates often involved in divergence are sexlinked genes (Qvarnström and Bailey 2009), which fits the observation that the sex determining region (SDR; Supporting Information 1.3) was isolated and even showed above average D xy for isolated regions Also, the enriched GO term detection of stimulus involved in sensory perception in isolated regions is in accordance to the observation that genes in possible relation to intraspecific communication play a disproportionately large role in divergence (Harr 2006;Lawniczak et al. 2010;Nadeau et al. 2012).

DYNAMICS OF THE DIVERGENCE PROCESS
Temporal divergence dynamics has been mostly inferred from multiple species or population pairs in different stages of divergence (Martin et al. 2013;Burri et al. 2015;Riesch et al. 2017). We analyzed the divergence dynamic from a single pair of taxa by using a molecular clock approach. Once gene flow and thus subsequent recombination is suppressed by selection against introgression, the affected genome regions should start to accumulate mutations in each of the divergent lineages independently. Knowing the mutation rate (Oppold and Pfenninger 2017), it was therefore possible to estimate the approximate time of divergence for each isolated region. Even though the variance of molecular clock estimates is large (Kumar and Blair Hedges 2016) and mutation rates (Hodgkinson and Eyre-Walker 2011), recombination rates (Nachman 2002) as well as selection are probably not constant over the entire genome, the distribution of estimates should nevertheless entail information on the temporal dynamics of divergence.
We used a coalescent approach to simulate expected divergence time distributions for five possible divergence scenarios and compared them with the empirically derived divergence time estimates. Applying an ABC inference framework, our approach shares the characteristic of all model selection approaches that it can only distinguish between the models tested, which may or may not include the true one (Johnson and Omland 2004).
The empirical data showed highest similarity with the episodic model, which assumed a short divergence period in the past (Fig. 1B). Calculating divergence time per window suggested that most of the isolated windows emerged during a relatively short time span about 4.6 million (3 N e ) generations ago. Assuming eight to 10 generations per year on average, this would place the divergence into the Mid-Pleistocene as a rough estimate. As the expected mean time (as well as the standard deviation) to coalescence of two alleles in a diploid population is 2 N e (Nei and Takahata 1993), we expect that lineage sorting among isolated loci is complete for the majority, but certainly not for all loci. This showed that the divergence process among the sister species is quite advanced and not recent. Such a scenario, reminding of a "punctuated equilibrium" as described by Eldredge and Gould (1972), in the case of divergence with gene flow is supported by theory: Simulations by Flaxman et al. (2014) have shown that initially gradual adaptive change can generate nonlinear transitions causing rapid emergence of complete RI, resulting in distinct bursts of speciation. Although our findings suggested that the divergence process conformed this model, it did not result in complete RI. It remained, however, unclear what triggered the divergence. An allopatric phase with differential ecological adaptation followed by secondary contact cannot be excluded. Another, nonexclusive explanation is that the beginning expansion of the Cla-element caused reduced gene flow by intragenomic incompatibilities ). However, our data did not suggest a close spatial association between large Cla clusters (or other repetitive elements) and divergence islands.
Factors influencing the length of isolated genome regions were further examined using a LM with selection, recombination, and divergence time as independent variables, as suggested by the Feder et al. model (2012b). Contrary to expectations (Martin et al. 2019), the recombination rate proxy explained only about 1% of length distribution variance. This seemed contradictory to the isolation of scaffolds known to contain the SDR and thus a region of known reduced recombination. However, reduced recombination may have simply contributed to maintain isolation of divergently selected sex-linked genes. Selection had the biggest impact, explaining about 9% of variance, followed by divergence time (6%). Albeit significant, all interactions between the variables had minor influence (<1%). The model as a whole explained only about 18%. This may have several, nonexclusive reasons. First, θ is known to be unbiased under neutrality and constant population growth (population growth rate constantly >1 through time) only (Subramanian 2016). These assumptions are probably not met in multivoltine taxa under a multivariate, fluctuating selection regime that undergo seasonal boom-bust cycles (Pfenninger and Foucault 2020). Also, LD and divergence time estimates lose accuracy under nonneutral conditions (Kumar 2005;Slatkin 2008). In addition, DH and GH may inherently lead to a heterogeneous differentiation pattern in divergence with gene flow scenarios Wolf and Ellegren 2017).
Likewise, pervasive selection by adaptive tracking in conjunction with the seasonally fluctuating population size in both species with an exponential increase over several generations in summer with a setback during a single generation in winter (Pfenninger and Foucault 2020) may explain the strongly skewed TD distribution toward negative values both in isolated and nonisolated regions. Simulations have shown that demographic fluctuations alone can produce such a pattern (Fig. S8), which would be amplified by constant selection (Excoffier et al. 2009). This bias made traces of divergently selected sites, typical for DH (Via 2012), harder to detect in the site frequency spectrum because they require more extreme negative TD values. Hence, our analysis might underestimate the portion of isolated windows showing signs of recent divergent selection (7% and 9% in C. riparius and C. piger, respectively). Other mutually nonexclusive explanations are (i) the majority of divergent selective sweeps occurred so long ago that most of the isolated regions are back in mutation-drift equilibrium or (ii) the sweeps were soft and thus did not produce a strong signal. Against this overall trend, a substantial fraction of isolated windows in both species, with a highly significant overlap, showed a lack of rare alleles (TD > 0 and exceeding the upper expectation threshold), indicating balancing selection. Both may be linked to multivoltinity of the species in which successive generations experience very different selection regimes throughout the seasons and over time (Foucault et al. 2019;Pfenninger and Foucault 2020).
MKT results indicated that most of the genes within isolated regions were not affected by (divergent) protein evolution. Thus, most of the genes within divergent regions either hitchhiked with the positively selected genes, which is most likely given the high gene density of the Chironomus genome, and/or hitchhiked with selection acting on expression control in noncoding parts.
Overall, parts of the isolated regions (still) showed signs of divergent and/or positive selection that fostered divergence. However, the isolation likely is kept up for a larger portion by background selection. Only a very small fraction of isolated windows (1.3%) showed signs of recent divergent selection in both species. It appears that the taxa are in a stable steady-state selectionmigration-drift equilibrium corresponding to phase 2 (DH) of the four-phase model of speciation with gene flow (Feder et al. 2012b. Similarly, Burri et al. (2015) found linked selection to lead to the formation of heterogeneous divergence landscapes in diverging Ficedula flycatcher species. But at the same time, they identified local recombination rate to have a strong effect on the build-up of isolation throughout the genome and not only in regions with known reduced recombination like, for example, the SDR in Chironomus. Avian species are known for their particularly variable recombination rate suggesting that its influence on isolation may be considered taxon specific rather than general (Kawakami et al. 2014). Similar to our results, studies on stick insects (Riesch et al. 2017) and flycatchers (Burri et al. 2015) question the concept of constantly growing isolation regions during the divergence process inevitably reaching complete RI, contrasting to theoretical expectations (Wu 2001;Feder et al. 2012b). It seems that, despite incomplete RI, species identity in these cases is maintained by other pre-and postzygotic reproductive barriers such as ecological selection (Schwenk et al. 2000), behavioral differences (Canestrelli et al. 2017), or genetic incompatibilities (Foucault et al. 2019).

CONCLUSION
The steadily accumulating evidence for speciation with gene flow reaching a stable intermediate steady state suggests the widespread occurrence of selection-migration-drift equilibria. Closely related taxa may thus share significant parts of their genomes over extended evolutionary times while retaining their specific identity.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1: Map of Europe displaying sample sites: MG: Hesse, Germany (50.188297, 9.214170); NMF: Lorraine, France (49.008298, 6.223934); MF: Lyon, France (46.017340, 4.911689); SI: Piemont, Italy (45.248064, 7.610818); SS: Andalusia, Spain (37.807905, -4.765055). Figure S2: Cumulative factor score sum of SNPs associated with (species splitting) PC1 in the PCA used to identify the isolated windows (see Fig. 2a). Figure S3: LD decay curve of maximum likelihood estimates of r2. Figure S4: Experimental design for coalescence simulations. Figure S5: Maximum likelihood estimation of individual ancestry based on 900,345 SNPs (quality filtered, LD-pruned) using ADMIXTURE version 1.3.0 with K = 2. Figure S6: Principal component analysis of the LD-pruned SNP dataset. Figure S7: Inference of lower temporal resolution limit with coalescence simulations. Figure S8: Influence of fluctuating population size on distribution of Tajima's D. Figure S9: Comparisons of nonisolated versus isolated parts of the genome in view of factors influencing the length of isolated regions. Table S1: Barcoding results for all 36 samples based on universal mitochondrial COI (Folmer et al. 1994), taxa-specific nuclear L44 (Oppold et al. 2016) as well as microsatellites markers. Table S2: Gene ontology (GO) term enrichment analysis of positively selected genes found in isolated parts of the genome. Table S3: Mean values of all four summary statistics per divergence scenario describing the divergence time distributions used in the ABC model inference approach as well as the empirical data. Table S4: Confusion matrix based on leave-on-out cross validation for 100 samples for each of the five models, applying a tolerance rate of 0.01 and using the "rejection" method.