Large inversions in Lake Malawi cichlids are associated with habitat, lineage, and sex determination
Figures
Overview of Lake Malawi cichlids.
A summary of the Lake Malawi cichlids and the results from this paper. On the left, a phylogeny showing the seven major lineages that currently live in the lake which are referred to as ecogroups. AC is Astatotilapia calliptera, a riverine species that lives in the margins of the river/lake interface. For each ecogroup, we summarize the (1) estimated number of species in the lake, (2) the number of species we have Bionano optical read data for, (3) the number of species we have PacBio based genome assemblies, (4) the number of species with Illumina short-read sequencing data, and (5) the primary habitat these ecogroups live in. Deep and shallow benthics are distinguished by the lake bottom depths where they live while utakas live in the water column but mate over the shore. We also summarize the distribution of six large inversions we have identified, colored by the chromosome they are found on (legend). We propose that the inversions on 9, 11, and 20 originated in the pelagic lineage and spread to the benthics by hybridization very early on in the evolution of the benthic lineage. Alternatively, these inversions distribution could be explained by incomplete lineage sorting.
Generation of a new high-quality reference genome for Metriaclima zebra.
Alignment of Bionano maps derived from a single Metriaclima zebra male to two different reference sequences for Chromosome 17 (white and gray, respectively). M_zebra_UMD2a reference (white) is the existing reference genome, and M_zebra_GT3a reference (gray) was created for this manuscript using a combination of PacBio HiFi reads and Bionano molecules. For each reference, the top bar shows the expected location of each CTTAAG motif (blue) created from an in-silico annotation of the reference sequence. These are aligned (gray lines) to observed locations in maps generated from Bionano Saphyr optical molecules (bottom bar). Deviations from co-linearity indicate reference errors (green, blue, and red highlights). Yellow indicates regions of the genome that lack aligned motifs. A summary of the improvements in the new genome is found in Supplementary file 1.
Identification of six large single or double inversions segregating within Lake Malawi cichlids.
Representative alignments of inversions or double inversions identified from blood samples obtained from 30 individuals from eleven species. For each inversion, the top shows predicted CTTAAG motifs (dark blue) from an in-silico annotation of the reference genome, the bottom shows motifs identified using Bionano molecules obtained from an individual of the species indicated on the left, and the gray lines indicate matching motifs based upon predicted and observed distances. Motifs in the reference or sample tracks with no matches are shown in green and light blue, respectively. Yellow indicates regions containing no matches between motifs in the reference genome and the mapped sample. Single inversions were identified on chromosomes 2, 9, 10, 13, and 20. Tandem double inversions were identified on chromosomes 11 and 20. The estimated length and percentage of the chromosome spanned by the inversion is shown on the right. The single inversion on 20 identified in Rhamphochromis longiceps has the same position as the first inversion of the double inversion on 20 identified from Mchenga conophoros. A list of all samples and their inversion genotype is in Table 1. Note that an error or polymorphism in the reference genome in chromosome 2 is indicated in the first panel.
Bionano maps from individuals that were heterozygous for the inversions on chromosomes 9, 10, or 11.
Bionano maps that support the non-inverted haplotype are shown in the middle while Bionano maps that support the inverted haplotype are shown below.
Genome alignments for chromosome 2.
Dot plots comparing the reference M_zebra_GT3a (x-axis) to six other genomes created from Lake Malawi cichlid fish (y-axis). Alignments are colored by scaffold names (see legends). MZGT3.2f and MZGT3.3f were each created from unique Metriaclima zebra female fish. YH_GT1.1f, YH_GT1.2m, and YH_GT1.3f were created from a single female, male, and female Aulonocara sp. ‘chitande type north’ Nkhata Bay fish, respectively. fRhaChi2.1 was a previously published genome generated from a Rhamphochromis sp. chilingali. Gray bar indicates location of the inversion. Short red fragments indicate regions with spurious short sequence matches which are primarily isolated to the assembly error region in chromosome 2 of the Mzebra_GT3 genome (see Figure 3).
Genome alignments for chromosome 9.
Dot plots comparing the reference M_zebra_GT3a (x-axis) to six other genomes created from Lake Malawi cichlid fish (y-axis). Alignments are colored by scaffold names (see legends). MZGT3.2f and MZGT3.3f were each created from unique Metriaclima zebra female fish. YH_GT1.1f, YH_GT1.2m, and YH_GT1.3f were created from a single female, male, and female Aulonocara sp. ‘chitande type north’ Nkhata Bay fish, respectively. fRhaChi2.1 was a previously published genome generated from a Rhamphochromis sp. chilingali. Gray bar indicates location of the inversion.
Genome alignments for chromosome 10.
Dot plots comparing the reference M_zebra_GT3a (x-axis) to six other genomes created from Lake Malawi cichlid fish (y-axis). Alignments are colored by scaffold names (see legends). MZGT3.2f and MZGT3.3f were each created from unique Metriaclima zebra female fish. YH_GT1.1f, YH_GT1.2m, and YH_GT1.3f were created from a single female, male, and female Aulonocara sp. ‘chitande type north’ Nkhata Bay fish, respectively. fRhaChi2.1 was a previously published genome generated from a Rhamphochromis sp. chilingali. Gray bar indicates location of the inversion.
Genome alignments for chromosome 11.
Dot plots comparing the reference M_zebra_GT3a (x-axis) to six other genomes created from Lake Malawi cichlid fish (y-axis). Alignments are colored by scaffold names (see legends). MZGT3.2f and MZGT3.3f were each created from unique Metriaclima zebra female fish. YH_GT1.1f, YH_GT1.2m, and YH_GT1.3f were created from a single female, male, and female Aulonocara sp. ‘chitande type north’ Nkhata Bay fish, respectively. fRhaChi2.1 was a previously published genome generated from a Rhamphochromis sp. chilingali. Gray bar indicates location of the inversion.
Genome alignments for chromosome 13.
Dot plots comparing the reference M_zebra_GT3a (x-axis) to six other genomes created from Lake Malawi cichlid fish (y-axis). Alignments are colored by scaffold names (see legends). MZGT3.2f and MZGT3.3f were each created from unique Metriaclima zebra female fish. YH_GT1.1f, YH_GT1.2m, and YH_GT1.3f were created from a single female, male, and female Aulonocara sp. ‘chitande type north’ Nkhata Bay fish, respectively. fRhaChi2.1 was a previously published genome generated from a Rhamphochromis sp. chilingali. Gray bar indicates location of the inversion.
Genome alignments for chromosome 20.
Dot plots comparing the reference M_zebra_GT3a (x-axis) to six other genomes created from Lake Malawi cichlid fish (y-axis). Alignments are colored by scaffold names (see legends). MZGT3.2f and MZGT3.3f were each created from unique Metriaclima zebra female fish. YH_GT1.1f, YH_GT1.2m, and YH_GT1.3f were created from a single female, male, and female Aulonocara sp. ‘chitande type north’ Nkhata Bay fish, respectively. fRhaChi2.1 was a previously published genome generated from a Rhamphochromis sp. chilingali. Gray bar indicates location of the inversion.
Genome alignments to Lake Malawi outgroups.
Dot plots comparing the reference M_zebra_GT3a (x-axis) to Oreochromis niloticus (red) or Pundamilia nyererei (shades of blue) (y-axis). Inversion location listed above each plot. Gray bar indicates location of the inversions and dark gray bars indicate the location of tandem inversion breakpoints.
Distribution of inversions within the Lake Malawi ecogroups.
Principal component analysis was used to analyze single nucleotide variants (SNVs) distributions in 365 samples within each inversion to genotype the inversion haplotype for all six inversions. (A) We illustrate the approach for genotyping the chromosome 10 inversions. The Principal Component Analysis (PCA) plot for the entire genome is shown on the left, with each sample colored by its ecogroup. Samples from each ecogroup cluster together, reflecting their evolutionary history. Individuals with Bionano data are shown as gray Xs. On the right, for the SNVs that fall within the chromosome 10 inversion, the deep benthic/shallow benthic and utaka individuals split into multiple clusters. Using the samples genotyped with Bionano data, we assign each cluster to an inversion genotype, which is shown by the black, green, and red boxes. To make comparisons of the whole-genome PCA and chromosome 10 PCA plots easier, we reversed the y-axis on the right panel. The assigned genotypes of each sample can be found in Supplementary file 2. Interactive PCA plots for each inversion are included in Figure 4—source data 1 to Figure 4—source data 9. (B) The derived allele frequency was calculated for each inversion within the seven ecogroups. The number of species that were genotyped and the estimated number of species that live in Lake Malawi are listed in the whole-genome phylogeny (left).
-
Figure 4—source data 1
Whole-genome Principal Component Analysis (PCA) of Lake Malawi cichlid samples across all ecogroups.
Principal component analysis of single nucleotide variants (SNVs) across the whole-genome for 365 Lake Malawi cichlid samples is shown. Each sample is colored by its ecogroup. Samples from each ecogroup cluster together, reflecting their evolutionary history. Individuals with Bionano data are shown as gray Xs. Hovering over individual points in the interactive HTML file reveals sample metadata, including PC space coordinates, Sample ID (see Supplementary file 2), and organism species name.
- https://cdn.elifesciences.org/articles/104923/elife-104923-fig4-data1-v1.zip
-
Figure 4—source data 2
Principal Component Analysis (PCA) of Lake Malawi cichlids across the chromosome 2 inversion.
Principal component analysis of single nucleotide variants (SNVs) within the chromosome 2 inversion region for 365 Lake Malawi cichlid samples is shown. Each sample is colored by its ecogroup. Individuals with Bionano data are shown as gray Xs. These Bionano samples in each cluster allow us to assign inversion genotype to other samples within that cluster. Individuals carrying the chromosome 2 inversion separate from individuals without the inversion in the same ecogroup, and the resulting clusters correspond to inversion genotype. Hovering over individual points in the interactive HTML file reveals sample metadata, including PC space coordinates, Sample ID, and organism species name. Supplementary file 2 summarizes the inversion status for each SampleID in the PCA.
- https://cdn.elifesciences.org/articles/104923/elife-104923-fig4-data2-v1.zip
-
Figure 4—source data 3
Principal Component Analysis (PCA) of Lake Malawi cichlids across the chromosome 9 inversion.
Principal component analysis of single nucleotide variants (SNVs) within the chromosome 9 inversion region for 365 Lake Malawi cichlid samples is shown. Each sample is colored by its ecogroup. Individuals with Bionano data are shown as gray Xs. These Bionano samples in each cluster allow us to assign inversion genotype to other samples within that cluster. Individuals carrying the chromosome 9 inversion separate from individuals without the inversion in the same ecogroup, and the resulting clusters correspond to inversion genotype. Hovering over individual points in the interactive HTML file reveals sample metadata, including PC space coordinates, Sample ID, and organism species name. Supplementary file 2 summarizes the inversion status for each SampleID in the PCA.
- https://cdn.elifesciences.org/articles/104923/elife-104923-fig4-data3-v1.zip
-
Figure 4—source data 4
Principal Component Analysis (PCA) of Lake Malawi cichlids across the chromosome 10 inversion.
Principal component analysis of single nucleotide variants (SNVs) within the chromosome 10 inversion region for 365 Lake Malawi cichlid samples. Each sample is colored by its ecogroup. Individuals with Bionano data are shown as gray Xs. These Bionano samples in each cluster allow us to assign inversion genotype to other samples within that cluster. Individuals carrying the chromosome 10 inversion separate from individuals without the inversion in the same ecogroup, and the resulting clusters correspond to inversion genotype. Hovering over individual points in the interactive HTML file reveals sample metadata, including PC space coordinates, Sample ID, and organism species name. Supplementary file 2 summarizes the inversion status for each SampleID in the PCA.
- https://cdn.elifesciences.org/articles/104923/elife-104923-fig4-data4-v1.zip
-
Figure 4—source data 5
Principal Component Analysis (PCA) of Lake Malawi cichlids across the first tandem inversion on chromosome 11.
Chromosome 11 contains a tandem inversion (see Figure 3); principal component analysis was performed on the first of these tandem inversions and the second separately. Here, PCA of single nucleotide variants (SNVs) within the first tandem inversion on chromosome 11 is shown for 365 Lake Malawi cichlid samples. Each sample is colored by its ecogroup. Individuals with Bionano data are shown as gray Xs. These Bionano samples in each cluster allow us to assign inversion genotype to other samples within that cluster. Individuals carrying the first tandem chromosome 11 inversion separate from individuals without this inversion in the same ecogroup, and the resulting clusters correspond to inversion genotype. Hovering over individual points in the interactive HTML file reveals sample metadata, including PC space coordinates, Sample ID, and organism species name. Supplementary file 2 summarizes the inversion status for each SampleID in the PCA.
- https://cdn.elifesciences.org/articles/104923/elife-104923-fig4-data5-v1.zip
-
Figure 4—source data 6
Principal Component Analysis (PCA) of Lake Malawi cichlids across the second tandem inversion on chromosome 11.
Chromosome 11 contains a tandem inversion (see Figure 3); principal component analysis was performed on the first of these tandem inversions and the second separately. Here, PCA of single nucleotide variants (SNVs) within the second tandem inversion on chromosome 11 is shown for 365 Lake Malawi cichlid samples. Each sample is colored by its ecogroup. Individuals with Bionano data are shown as gray Xs. These Bionano samples in each cluster allow us to assign inversion genotype to other samples within that cluster. Individuals carrying the second tandem chromosome 11 inversion separate from individuals without this inversion in the same ecogroup, and the resulting clusters correspond to inversion genotype. Hovering over individual points in the interactive HTML file reveals sample metadata, including PC space coordinates, Sample ID, and organism species name. Supplementary file 2 summarizes the inversion status for each SampleID in the PCA.
- https://cdn.elifesciences.org/articles/104923/elife-104923-fig4-data6-v1.zip
-
Figure 4—source data 7
Principal Component Analysis (PCA) of Lake Malawi cichlids across the chromosome 13 inversion.
Principal component analysis of single nucleotide variants (SNVs) within the chromosome 13 inversion region for 365 Lake Malawi cichlid samples is shown. Each sample is colored by its ecogroup. Individuals with Bionano data are shown as gray Xs. These Bionano samples in each cluster allow us to assign inversion genotype to other samples within that cluster. Individuals carrying the chromosome 13 inversion separate from individuals without the inversion in the same ecogroup, and the resulting clusters correspond to inversion genotype. Hovering over individual points in the interactive HTML file reveals sample metadata, including PC space coordinates, Sample ID, and organism species name. Supplementary file 2 summarizes the inversion status for each SampleID in the PCA.
- https://cdn.elifesciences.org/articles/104923/elife-104923-fig4-data7-v1.zip
-
Figure 4—source data 8
Principal Component Analysis (PCA) of Lake Malawi cichlids across the first tandem inversion on chromosome 20.
Chromosome 20 contains a tandem inversion (see Figure 3); principal component analysis was performed on the first of these tandem inversions and the second separately. Here, PCA of single nucleotide variants (SNVs) within the first tandem inversion on chromosome 20 is shown for 365 Lake Malawi cichlid samples. Each sample is colored by its ecogroup. Individuals with Bionano data are shown as gray Xs. These Bionano samples in each cluster allow us to assign inversion genotype to other samples within that cluster. Individuals carrying the first tandem chromosome 20 inversion separate from individuals without this inversion in the same ecogroup, and the resulting clusters correspond to inversion genotype. Hovering over individual points in the interactive HTML file reveals sample metadata, including PC space coordinates, Sample ID, and organism species name. Supplementary file 2 summarizes the inversion status for each SampleID in the PCA.
- https://cdn.elifesciences.org/articles/104923/elife-104923-fig4-data8-v1.zip
-
Figure 4—source data 9
Principal Component Analysis (PCA) of Lake Malawi cichlids across the second tandem inversion on chromosome 20.
Chromosome 20 contains a tandem inversion (see Figure 3); principal component analysis was performed on the first of these tandem inversions and the second separately. Here, PCA of single nucleotide variants (SNVs) within the second tandem inversion on chromosome 20 is shown for 365 Lake Malawi cichlid samples. Each sample is colored by its ecogroup. Individuals with Bionano data are shown as gray Xs. These Bionano samples in each cluster allow us to assign inversion genotype to other samples within that cluster. Individuals carrying the second tandem chromosome 20 inversion separate from individuals without this inversion in the same ecogroup, and the resulting clusters correspond to inversion genotype. Hovering over individual points in the interactive HTML file reveals sample metadata, including PC space coordinates, Sample ID, and organism species name. Supplementary file 2 summarizes the inversion status for each SampleID in the PCA.
- https://cdn.elifesciences.org/articles/104923/elife-104923-fig4-data9-v1.zip
Evolutionary history of large inversions violates species tree.
(A) Maximum likelihood trees built from either whole-genome single nucleotide variants (SNVs) (top) or within the 9, 11, and 20 inversions. Between 2–25 representative samples were selected for each ecogroup (see Supplementary file 2) and individuals that were heterozygous were excluded from this analysis. Clades were collapsed by ecogroup for visualization purposes. Full phylogenies are presented in Figure 5—figure supplements 1–8. The whole-genome tree indicates Diplotaxodon and Rhamphochromis species split first and formed their own clade. For each of the displayed inversions, the first split was between samples carrying the inverted vs. non-inverted genotypes. Additionally, the Diplotaxodon individuals now formed a clade with the benthics that carry the inversion. (B) Density plots of dXY values comparing benthic species and species of the indicated ecogroup. For the displayed inversions, the Diplotaxodon species were much closer in evolutionary distance when compared to the whole-genome plot.
Whole-genome phylogeny for Lake Malawi cichlids.
A maximum likelihood tree was created from single nucleotide variant (SNV) data across the seven major ecogroups of Lake Malawi cichlids. Bootstrap support values are shown below each node. Branch lengths are proportional to the number of SNVs per base pair.
Chromosome 2 inversion phylogeny for Lake Malawi cichlids.
A maximum likelihood tree was created from SNV data across the seven major ecogroups of Lake Malawi cichlids, restricted to the inversion region on chromosome 2. Bootstrap support values are shown below each node. Branch lengths are proportional to the number of single nucleotide variants (SNVs) per base pair.
Chromosome 9 inversion phylogeny for Lake Malawi cichlids.
A maximum likelihood tree was created from single nucleotide variant (SNV) data across the seven major ecogroups of Lake Malawi cichlids, restricted to the inversion region on chromosome 9. Bootstrap support values are shown below each node. Branch lengths are proportional to the number of SNVs per base pair.
Chromosome 10 inversion phylogeny for Lake Malawi cichlids.
A maximum likelihood tree was created from single nucleotide variant (SNV) data across the seven major ecogroups of Lake Malawi cichlids, restricted to the inversion region on chromosome 10. Bootstrap support values are shown below each node. Branch lengths are proportional to the number of SNVs per base pair.
Chromosome 11 inversion phylogeny for Lake Malawi cichlids.
A maximum likelihood tree was created from single nucleotide variant (SNV) data across the seven major ecogroups of Lake Malawi cichlids, restricted to both of the tandem inversions on chromosome 11. Bootstrap support values are shown below each node. Branch lengths are proportional to the number of SNVs per base pair.
Chromosome 13 inversion phylogeny for Lake Malawi cichlids.
A maximum likelihood tree was created from single nucleotide variant (SNV) data across the seven major ecogroups of Lake Malawi cichlids, restricted to the inversion region on chromosome 13. Bootstrap support values are shown below each node. Branch lengths are proportional to the number of SNVs per base pair.
Chromosome 20 first tandem inversion phylogeny for Lake Malawi cichlids.
A maximum likelihood tree was created from single nucleotide variant (SNV) data across the seven major ecogroups of Lake Malawi cichlids, restricted to the first tandem inversion (20 a) region on chromosome Bootstrap support values are shown below each node. Branch lengths are proportional to the number of SNVs per base pair.
Chromosome 20 second tandem inversion phylogeny for Lake Malawi cichlids.
A maximum likelihood tree was created from single nucleotide variant (SNV) data across the seven major ecogroups of Lake Malawi cichlids, restricted to the second tandem (20b) inversion region on chromosome Bootstrap support values are shown below each node. Branch lengths are proportional to the number of SNVs per base pair.
A role for the chromosome 10 inversion in sex determination.
(A) Schematic of the two inversion states (normal and inverted) and their assignment as X or Y chromosomes. (B) 24 offspring (12 male and 12 female) from two separate male/female pairs were genotyped for the inversion. For 47 of 48 individuals, the XY animals were male, and the XX animals were female. One male was XX. (C) FST analysis comparing 24 males and 24 females identified chromosome 10 as the primary genetic region that segregated with sex. (D) Comparison of heterozygosity levels of 24 male and 24 female animals also identified chromosome 10 as an outlier.
Tables
Inversion genotypes of 31 samples based on Bionano data.
| Inversion | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| Species | Ecogroup | # male/female | 2 | 9 | 10 | 11 | 13 | 20 | Notes |
| Cynotilapia zebroides ‘Cobue’ | Mbuna | 1 m/1 f | - | - | - | - | - | - | |
| Labeotropheus fuelleborni | Mbuna | 1 m/1 f | - | - | - | - | - | - | |
| Labeotropheus trewavasae | Mbuna | 1 m/1 f | - | - | - | - | - | - | |
| Metriaclima zebra | Mbuna | 3 m/1 f | - | - | - | - | - | - | Reference species |
| Pseudotropheus demasoni | Mbuna | 1 m/1 f | - | - | - | - | - | - | |
| Copadichromis virginalis | Utaka | 1 m/2 f | - | X | - | -* | - | X | * One female was het on 11 |
| Protomelas taeniolatus | Shallow benthic | 1m | - | X* | - | - | - | X | * The male was het on 9 |
| Mchenga conophoros | Shallow benthic | 1 m/2 f | - | X | - | X | - | X | |
| Nyassachromis prostoma ‘Orange Cap’ | Shallow benthic | 1 m/1 f | - | X | - | X* | - | X | * The male was het on 11 |
| Aulonocara sp. ‘chitande type north’ Nkhata Bay | Deep benthic | 3 m/1 f | X | X | X* | X | X | X | * All three males were het on 10 |
| Diplotaxodon limnothrissa | Diplotaxodon | 2 m/0 f | - | X | - | X | - | X | |
| Rhamphochromis longiceps | Rhamphochromis | 3 m/1 f | - | - | - | - | - | X* | * All samples had the single left inversion on 20 |
Additional files
-
Supplementary file 1
Genome statistics for the M_zebra_UMD2a and M_zebra_GT3a reference genomes.
Per Linkage Group (LG; chromosome) anchored nucleotide content is compared between the M_zebra_UMD2a and M_zebra_GT3a genome assemblies.
- https://cdn.elifesciences.org/articles/104923/elife-104923-supp1-v1.docx
-
Supplementary file 2
Sample Metadata.
Metadata for all samples used in this study. Where known, sample sex, parents, ecogroup, species name, whether the sample was lab-reared, and average whole-genome sequencing coverage are given. Columns F–I indicate whether a sample was included in the analysis that produced the referenced figure. Columns L–R give inversion haplotypes for each of the six inversions, with the tandem inversions on chromosome 20 separated into the first (20 a) and second (20b) regions. Column S indicates which samples were used to generate the initial principal components onto which the remaining samples were projected.
- https://cdn.elifesciences.org/articles/104923/elife-104923-supp2-v1.xlsx
-
Supplementary file 3
Summary statistics for additional de novo genome assemblies for Metriaclima zebra and Aulonocara sp.
‘chitande type north’ Nkhata Bay. Assembly statistics for the five scaffold-level de novo assemblies generated from PacBio HiFi reads and Bionano optical maps. Two are from Metriaclima zebra (MZ_GT3.2, MZ_GT3.3) and three from Aulonocara sp. 'chitande type north' Nkhata Bay (YH_GT1.1, YH_GT1.2, YH_GT1.3). The number of scaffolds, total assembly size (bp), and N50 (bp) are given for each. Assembly GenBank accessions are listed in the Data Availability section.
- https://cdn.elifesciences.org/articles/104923/elife-104923-supp3-v1.docx
-
Supplementary file 4
Estimates of breakpoint locations in M_zebra_GT3a.
Breakpoint intervals for the six large inversions are given in M_zebra_GT3a coordinates (bp). Interval estimates were calculated using whole genome alignments between Aulonocara sp. ‘chitande type north’ Nkhata Bay and the M_zebra_GT3a genome. The breakpoint intervals for the tandem inversions on chromosomes 11 and 20 are separated by the first (11 a/20 a) and second (11b/20b) inversions.
- https://cdn.elifesciences.org/articles/104923/elife-104923-supp4-v1.docx
-
MDAR checklist
- https://cdn.elifesciences.org/articles/104923/elife-104923-mdarchecklist1-v1.docx