Introgression dynamics of sex-linked chromosomal inversions shape the Malawi cichlid radiation

Abstract

Chromosomal inversions can contribute to adaptive speciation by linking coadapted alleles. By querying 1375 genomes of the species-rich Malawi cichlid fish radiation, we discovered five large inversions segregating in the benthic subradiation that each suppress recombination over more than half a chromosome. Two inversions were transferred from deepwater pelagic Diplotaxodon through admixture, whereas the others established early in the deep benthic clade. Introgression of haplotypes from lineages inside and outside the Malawi radiation coincided with bursts of species diversification. Inversions show evidence for transient sex linkage, and a notable excess of protein changing substitutions points toward selection on neurosensory, physiological, and reproductive genes. These results indicate that repeated interplay between depth adaptation and sex-specific selection on large inversions has been central to the evolution of this iconic system.

Five large chromosomal inversions contribute to the diversification of Malawi cichlids.

Inversions established in the diverse benthic subradiation. Inversion region haplotypes were exchanged through hybridisation of lineages within and outside of the Malawi radiation and contribute to ecological and habitat divergence, sensory adaptation, and sex determination.

Fig. 1. Study system and prevalence of five large inversions.

Consensus phylogeny of Malawi cichlid species used in this study (data S1) with inversion frequency based on WGS and PCR-typing (see materials and methods) (38) shown in rings around the phylogeny (the same colours are used throughout the article). The benthic subradiation is expanded to show the phylogenetic position of each species and to highlight subclades that we refer to in the main text (Shallow rocky Aulonocara, Shallow Lethrinops, Eukambuzi). Note that utaka are not monophyletic in this phylogeny. Non-benthic groups of Malawi cichlids (i.e., the pelagic subradiations of Rhamphochromis and Diplotaxodon, the subradiation of predominantly rock-dwelling mbuna, and Astatotilapia calliptera – a species distributed in rivers and margins around the lake that shares its putatively ancestral characteristics and genus assignment with riverine haplochromines outside of the radiation) are each represented by a single grey triangle approximately reflecting species richness relative to each other. Dashed lines indicate branches with unstable placement. See data S1 and S2 for full phylogenies with branch lengths and support values. Annotations next to species/clade names provide the numbers of sequenced/inversion-genotyped samples (additional samples which were inversion-genotyped with PCR are indicated as + n in the annotation). Two taxa are annotated with ‘+’ to denote polyphyletic groups: Otopharynx argyrosoma contains a single Cyrtocara moorii individual and Ctenopharynx intermedius contains two Ctenopharynx pictus individuals. Full species names are given in table S3 and inversion frequencies by species in table S2. Species names for representative photographs are given in fig. S1. Tree files are given in data S1 and S2. Species subject to further experimental investigation (see text) are highlighted in bold.

Fig. 2. Characterisation of inversions.

(A) (Top panel): Identification of genomic regions from clusters of aberrant phylogenetic patterns (see materials and methods) (38). (Bottom panel): First genetic principal component in overlapping 1 Mbp windows along chromosomes, using the same colours for the benthic subclades as in Fig. 1. Outlier regions from the top panel are highlighted and colour-labelled. Centromeric satellite regions (for inference see materials and methods, text S1, fig. S17, table S6) (38) are indicated as black rectangles on top of the X axis. (B) Representative photographs of the species used in panels C-D: Astatotilapia calliptera, a lineage of the Malawi radiation distinct from benthics from which the reference genome was produced, and Aulonocara stuartgranti, a species that genetically belongs to the deep benthic group, but lives in shallow rocky habitats (clade Shallow rocky Aulonocara in Fig. 1). According to WGS-typing, the species are expected to show opposite orientations for the chromosome 9 and 11 inversions. (C) Fluorescence in situ hybridisation (FISH) of markers on chromosome 11 left and right of the putative inversion breakpoints show the expected non-inverted orientation (upper panel) in A. calliptera. In Au. stuartgranti we see a double inversion (lower panel; see fig. S8 for FISH of chromosome 9). (D) Whole genome alignment of an ONT duplex long-read assembly of Au. stuartgranti to the A. calliptera reference assembly (which was re-scaffolded with chromosome conformation capture (Hi-C) data, see materials and methods) (38) confirms the double inversion on chromosome 11 (for other chromosomes see fig. S9). (E) Top: Windowed PC1 values of whole genome sequenced founders and progeny of an interspecific cross. Among 290 F2 and F3 individuals no crossing-over events were observed in the inversion regions of chromosomes 9 (bottom left) and chromosome 11 (bottom right), while recombination was frequent in the flanking regions and on other chromosomes.

Fig. 3. Evolutionary history scenario of inversion haplotypes.

(A) Density plots of pairwise sequence divergence translated into divergence (coalescence) times assuming a mutation rate of 3×10^-9 bp per generation (31). The top panels show results for the genome outside the five large inversions, comparing all major clades against shallow benthics (left) and deep benthics (right). Panels below the top row show divergence in inversion regions for the non-inverted (left) and inverted (right) benthic haplotypes. (B) A simplified model for the evolutionary history of the Malawi cichlid radiation, which includes several inversion haplotype transmission events. Vertical grey connections indicate major gene flow events. Letter-labelled arrows indicate transfer of inversion-region haplotypes. For further events see fig. S34. Lineages in which re-introgressed inversion region haplotypes of ancestral orientation apparently play a Y-like role in sex determination (see Fig. 5 and main text) are indicated by 🅈. (C) Evidence for transfer of inversion haplotypes through introgressive hybridisation. Histograms of ABBA-BABA statistics D(P1, P2, P3, Outgroup) calculated outside the inversions. For the different panels, we selected those ABBA-BABA tests for which the inverted state of the respective chromosome is present in one of the two more closely related species P1 and P2 but absent in the other and ordered them such that P2 shared the inversion state (presence/absence) with P3. In such a configuration, significantly positive values are suggestive of gene flow outside of inversions between the species sharing inversion states, while significantly negative values suggest gene flow between species not sharing inversion states. Under the null hypothesis of no inversion introgression, the statistic would be symmetric around zero.

Fig. 4. Adaptive evolution of inversion haplotypes.

(A) Proportion of exonic SNPs, grouped by inversion correlation coefficient intervals, relative to all SNPs within the same interval. A positive correlation coefficient corresponds to the derived SNP allele being more common on the inverted haplotype, while a negative coefficient corresponds to the derived SNP allele being more common on the non-inverted haplotype (ancestral orientation). (B) dN/dS measured for SNPs as a function of inversion genotype correlation (see materials and methods) (38). (C) Excess of genes containing non-synonymous highly inversion correlated SNPs (nsICS) among all genes highly expressed in the main sensory (green) and nervous system (blue). Expression data is based on the single-cell expression atlas of developing zebrafish Daniocell (45). The tissues were grouped into the functional categories: vision (eye), mechanoreception (lateral line, ear), chemoreception (taste, olfaction), and nervous system (neural). (D) Non-synonymous SNPs (grey dots; if ICS: empty dots) and averaged ICS scores in 100 SNP rolling windows (markers colour-coded according to inversion) on the five inversion chromosomes. We annotated nsICS located in genes with high expression in zebrafish tissue groups related to sensory perception (same as in (C)). nsICS in candidate genes discussed in the main text are annotated with arrows (if several were located in the same gene, the highest nsICS is annotated) and color-coded by functional category (table S24).

Fig. 5. Sex association of inversions.

(A) Windowed PC analysis along chromosome 11 demonstrates perfect association of inversion genotype with sex in our sample of 28 wild-caught Copadichromis chrysonotus. (B) SNP heterozygosity among our sample of wild-caught male and female C. chrysonotus, measured as number of het SNPs per 10 kbp. (C-E) Sex-inversion associations in lab-raised populations. Per population, the number of females and males is given (separated by ‘|’) and asterisks denote significance levels of Fisher’s exact tests of inversion genotype–sex correlation (*: p < 0.05, ***: p < 0.001). Inversion genotype per sex (confirmed through gonad examination) in lab-raised broods of (C) C. chrysonotus from Lake Malombe (p < 0.001, left) and from Lake Malawi (p = 0.048, right), (D) Lethrinops chilingali from satellite lake Chilingali (p < 0.001) (E), Otopharynx tetrastigma from Lake Malombe at the outflow of Lake Malawi (p < 0.001, left) and from the northern part of Lake Malawi (p = 0.049, right). (F) Proportions of homozygous/heterozygous inversion genotypes in males and females of species with heterozygotes present, according to WGS and PCR typing of 809 samples (67 species). Asterisks denote significance levels of Fisher’s exact tests of inversion genotype–sex correlation (*: p < 0.05, ***: p < 0.001).