An intronic transposon insertion associates with a trans-species color polymorphism in Midas cichlid fishes

Conflicting results suggest a missing variant

In order to narrow down candidates for the causal genetic variant, we performed genome-wide association mapping separately in individual lake populations (previously, association mapping was only performed across the whole species flock5). Interestingly, despite clear association peaks in the crater lakes (Fig. 1a, b), the exact position of the top-associated variants differed by ca. 95 kb between CL Masaya [11:7,064,181] and CL Xiloá [11:7,159,113]. Even more surprisingly, both great lake populations completely lacked strongly supported marker-trait associations (Fig. 1c, d). A possible explanation for both observations is that we solely identified linked variants and did not genotype the causal variant itself. This could be, for example, due to a structural variant (SV) that is either missing completely from the reference genome or is hard to identify due to limitations of the used methods (i.e., reads cannot be mapped or not reliably genotyped using short-read sequencing). In fact, SVs often remain uncovered, although they are a substantial source of population genetic variation and can have direct functional implications8. Owing to their relatively old population ages and large population sizes, linkage blocks are rather small in the Midas cichlid great lake populations. In these situations, it is expected that almost no association to neighboring markers will be detected, even to those that are in close proximity. Moreover, this could also explain the shifted association peaks in the crater lakes, as the most-highly associated linked variants can differ depending on the demographic history of the populations.

Strong association of a newly identified structural variant

To uncover structural variants that were possibly absent from the current reference genome (derived from a dark morph individual), we generated a de novo assembly of a lab-raised heterozygous golden individual of an F4 inbred line using accurate long reads (PacBio HiFi). The assembly generated 817 primary scaffolds (N50 of 2.47 Mb), with a total size corresponding to 95% of the Midas cichlid reference genome4. Both haplotypes at the gold locus were recovered, determined based on known variants at flanking sites that differed between the parents of the cross. The assemblies of the G and d haplotypes included the candidate interval, as well as hundreds of kilobases of flanking sequences (total size of the alignment of both haplotypes was 1.3 Mb). This new long-read alignment led to the discovery of an 8.2 kb insertion in the G haplotype at position 11:7,069,471 of the reference genome (Fig. 1f)—located in close proximity to the most-highly associated SNPs from our previous genome-wide association analysis (5.7 kb downstream; 11:7,063,765) and positioned between the peaks found in CLs Masaya (5.3 kb downstream; 11:7,064,181) and Xiloá (89.6 kb upstream; 11:7,159,113).

To test whether this newly identified structural variant is more strongly associated across all populations than the previously identified SNP/short indel markers, we added its genotype information to our previous genotype call set (see methods). When we repeated the genome-wide association mapping (Fig. 1e) we now obtained substantially stronger support for this variant (LRT test, P = 2.23 × 10−76; in comparison: the previously most-highly associated SNP using the same data except the one variant had statistical support of P = 7.96 × 10−27). Moreover, we also found strong support for this variant in the source lake populations from GL Nicaragua and GL Managua (Fig. 1h, i), where no strongly associated markers were found in the initial analysis (Fig. 1c, d). The overall association is near perfect (Fig. 1g, Supplementary Fig. 1), with only very few conflicting results explained by the characteristics of the focal trait. There are eight dark individuals with one copy of the G allele that are likely not yet transformed (4.3%; 8 of 185 individuals; Fig. 1b, Supplementary Fig. 1). Notably, all of these are heterozygotes, which is consistent with previous data showing that the onset of color change from dark to gold during ontogeny is dosage-dependent, with heterozygous fish transitioning much later in life (hence untransformed adult individuals would more likely carry the Gd genotype)5. Only a single genetically homozygous dark (dd) phenotypically golden individual was detected (0.6%; only 1 in 88 individuals, Fig. 1g, Supplementary Fig. 1). This represents very likely a genetic or non-genetic phenocopy.

goldentouch is differentially expressed between morphs

The 8.2 kb insertion is located within an intron of goldentouch, a previously undescribed gene (Fig. 2a). Based on these new results and identification of the likely causal variant we hypothesized that the transposon insertion might affect expression of goldentouch (and/or possibly other genes in genomic proximity). To specifically screen for these genes as well as to identify downstream effector genes, we performed RNA-seq on scale tissue from adult gold (n = 9) and dark individuals (n = 9). Fish scales are covered by epidermal and dermal tissue that contain pigment cells including the dark melanophores that undergo cell death during the transition from dark to gold coloration6. In total, we found 247 differentially expressed genes (Fig. 2b Supplementary Data 1; Padj < 0.05; total number of coding genes: 22,495). To screen for cis-regulatory variation at the gold locus, we specifically checked for differentially expressed genes within a 4 Mb window (variant ± 2 Mb; to also include long-range cis-regulated genes9) around the insertion site. Only one of the 96 genes in this region showed differential expression (Fig. 2c, d): the gene goldentouch that harbors the transposon insertion. A previous candidate gene underlying the gold/dark polymorphism, the gene stk that is directly 5′ of the locus did not show support for differential expression (P = 0.99). To provide further confirmation for the differential expression of goldentouch, we performed quantitative PCR (qPCR) in homozygous gold and dark individuals (n = 6 for each genotype). These results confirm a 2.4 times lower expression of this gene in the gold morph (Fig. 2e; Welch Two-Sample t-test [two-sided]; P = 0.018; t = −3.04, df = 7.05, 95% CIs, −0.037, −0.004).

Fig. 2: The genetic basis of the gold/dark polymorphism.

a The structural variant that likely underlies the gold/dark polymorphism is an 8.2 kb transposon insertion in the last intron of an heretoforth undescribed gene that we here coin goldentouch. b–d RNA-seq on scales of adult fish reveals 247 differentially expressed genes between gold and dark individuals (Padj < 0.05; Benjamini–Hochberg adjusted Wald test P-value), including several pigmentation genes (b). The gene goldentouch is the only differentially expressed gene found within a 4 Mb interval around the insertion (c, d). e Quantitative PCR confirms the decreased expression of goldentouch in golden individuals (center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, data points; Welch’s two-sample t-test, two-sided). f Relative change of gene expression during ontogeny in golden individuals. Expression of goldentouch (red) remains relatively stable, compared to melanophore marker genes (blue) that decrease in expression due to cell death of the melanophores in which they are expressed (data are presented as mean values ± SD; statistical test: one-way ANOVA). g Gene ontology (GO) term enrichment analysis of genes differentially expressed between gold and dark individuals.

goldentouch expression does not change during ontogeny

As the morphological color change in Midas cichlids has an ontogenetic progression it is important and relevant to ask what genes change their expression during ontogeny. This has been investigated and many of the same pigmentation genes, as found here, have been previously identified7. We reanalyzed this data set using only genotypically golden individuals (before transition n = 5, transitional n = 6, after transition n = 6) using the Midas cichlid reference genome (the initial analysis7 used two independent approaches; a de novo assembly and an assembly against the Nile tilapia transcriptome) to test two hypotheses. First, we tested if there are additional differentially expressed genes in any of the three pairwise comparisons within or close to the gold interval (2 Mb upstream and downstream). Although we found the same pigmentation genes (Supplementary Fig. 2a) that had been identified previously with this data set7, we did not find any gene in or close to the gold interval (all genes in 4 Mb interval P > 0.05). Second, we wanted to specifically test if goldentouch expression changes during the process of color change (Fig. 2g). In contrast to melanophore marker genes (ANOVA, P = 0.0003, P = 0.018, P = 0.0001 for mreg, tyrp1b, and pmela; Fig. 2f) we did not find any significant expression changes for goldentouch throughout our ontogenetic sampling (ANOVA, P = 0.37; Fig. 2f). This lack of change was robust (ANOVA, P = 0.31; Supplementary Fig. 2b) even after accounting for genotype and only considering heterozygous Gd individuals (using markers in neighboring genes; sample size n = 5, n = 4, n = 3 for before, during and after transition). Furthermore, we found no significant correlation between melanin-related genes and goldentouch expression (linear regression, P = 0.48, P = 0.19, P = 0.28 for mreg, tyrp1b, and pmela; Supplementary Fig. 2c-h). This suggests that goldentouch is unlikely to be expressed in melanophores themselves but in the surrounding skin tissue as the expression does not decrease with declining melanophore numbers. This data also provide insights into the mechanism as it might suggest that there is a constitutive, genotype-dependent expression of goldentouch that non-cell autonomously acts on melanophore survival.

RNA-seq reveals downstream effectors of goldentouch

To further explore the potential mechanisms behind the effect of goldentouch on melanophore (dark pigment cells) survival, we also examined potential downstream effects on genes outside the gold interval. Consistent with the observation that color transition involves the loss of cell types during the ontogeny in genetically gold individuals (e.g., melanophores), most differentially expressed genes (89%; 220 genes) were more highly expressed in the dark morph. Only 27 genes (11%) had higher expression in gold individuals (Fig. 2b). Indeed, many of the genes with lower expression in gold individuals are well-known melanophore-specific genes and key players in melanin synthesis (Fig. 2b), which is consistent with previous findings7 and provides internal validation for our results. A gene ontology (GO) term analysis of the differentially expressed genes (Fig. 2g; Supplementary Table 1 and Supplementary Data 1) reveals an overrepresentation of genes associated with metabolism and extracellular space suggesting that non-cell-autonomous effects might indirectly lead to melanophore cell death. While the involvement of metabolic and proteolytic processes is more difficult to explain (it might indicate changes in skin homeostasis that indirectly affect melanophores10), interactions between cells and between cells and extracellular matrix are important for pigment cell migration and survival11,12,13. Yet, the results generally hint at a complicated molecular mechanism underlying the morphological color change that might be rather due to a cumulative effect of non-cell-autonomous factors and their interaction with maturation and aging processes.

Phylogenetic and structural analysis of goldentouch

As the goldentouch gene had not been described before, we compared its coding sequence to the non-redundant NCBI nucleotide collection (nr/nt) using the tblastx algorithm to investigate whether (i) it is indeed a coding gene and (ii) whether it is also present in other species. We could identify orthologs of goldentouch in more than 20 teleost species, all of which belong to the Order Percomorpha (sensu ref. 14). Thus, we suggest that the gene is a Percomorpha-specific gene that originated ~120 million years ago14 (Fig. 3a). Owing to the lack of similarity to any other known genes it is possible that its origin was de novo. Next, we predicted the structure of goldentouch (Fig. 3b) and screened for more distant protein homologies and known domains. Our analysis gives weak support for a hydrolase/transport protein domain (confidence that the match between your sequence and this template is a homology: 55.0 on a scale from 0 to 100; Proportion of protein residues equivalenced to identical template residues in the generated alignment: 53%). However, the exact function of the protein within the cells of the dermis surrounding the melanophores remains unclear. Lastly, we analyzed the expression pattern of this gene across several tissues. Its expression was almost absent from most internal organs and highest in scales, supporting a specific function in the outer skin layers. (Fig. 3c).

Fig. 3: Phylogenetic, structural, and expression analyses of goldentouch.

a Bayesian phylogenetic reconstruction of goldentouch, which was only found in the order Percomorpha and could, for example, not be detected in zebrafish, Danio rerio, a non-Percomorpha. b Three-dimensional structure prediction of goldentouch. c Expression of goldentouch is mostly restricted to scales, with some expression in the underlying skin and in spleen tissue (center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers).

The insertion is an inverted repeat with a 4.1 kb stem

A further important question is how the transposon insertion might downregulate goldentouch expression. Comparing the sequence against the non-redundant NCBI nucleotide collection (nr/nt) using blastx revealed that the insertion is transposon-derived and shows high similarity to a DNA or class II transposon, PiggyBac Transposable Element Derived 4. Interestingly, the transposon is duplicated in the sequence and thereby appears as a symmetrical structure in a plot of the blastx result: a dotplot revealed a secondary structure with two hairpins with a 4.1 kb stem (Fig. 4a). Such structures that form cruciforms (Fig. 4b) have been shown to affect DNA supercoiling and nucleosome positioning15. The stalling of Polymerase II elongation would lead to mRNA without poly-adenylation likely resulting in unstable, quickly degraded mRNA, explaining the lower expression in golden individuals. While RNA without poly-adenylation would have not at all been picked up by our RNA-seq experiment, the formation of a large hairpin on the pre-mRNA would be expected to decrease translation. This is supported by our qPCR results, suggesting that it affects mRNA stability (and or expression itself), either through lack of the polyA tail or through premature termination due to the hairpin structure of the repeat. It remains to be investigated whether this SV also affects other features of this locus including somatic as well as meiotic mutation rates as well as recombination rate as had been suggested as a consequence for long inverted repeat (LIR) sequences16. DNA properties (in this case Z-DNA formation) have been, for example, recently shown to greatly affect the dynamics and to facilitate adaptive repeated evolution in sticklebacks17.

Fig. 4: Structural properties of the structural variant in gold Midas cichlids.

a A dotplot of the golden allele at the genomic position of the goldentouch gene. The structural variant we identified based on a haplotype-resolved PacBio HiFi assembly of a heterozygous golden individual is an 8.2 kb transposon-derived inverted repeat with identical copies of a piggyBac transposable element (TE). b The two alleles likely differ in their secondary structure as inverted repeats are known to form cruciform structures. The part of the two copies that form the stem of the DNA cruciform have sequence similarity and differ only at two positions (two single-base-pair indels).

Additionally, independent of the effect on DNA conformation, the presence of the transposon within an actively transcribed region might affect goldentouch expression. Antisense transcription of transposons is also effective in blocking sense transcription18. Moreover, many mechanisms exist (including heterochromatin spreading and by stalling Polymerase II elongation) by which host cells actively silence TE expression, as well as that of the genes into which they integrated19. Furthermore, given the established and drastic implications of such an inverted-repeat TE, it is likely that a multitude of mechanisms, such as those that we proposed above, play a role. Functional validation, including knockout of the transposable element—which is challenging given its size—would help to establish a causal link that could confirm our findings.

The developmental mechanism of the gold/dark polymorphism

Although the causal link needs further confirmation and the exact molecular mechanisms remains elusive up to now, our results seem to support the notion that secondary DNA structures and/or transposable elements can be an important driver of (adaptive) evolution. Other cases in the last decades have illustrated how transposons can affect gene regulation and phenotypic evolution: Transposon insertions have been, for example, associated with insecticide resistance in Drosophila20, industrial melanism in the peppered moth21, and sex-linked color polymorphisms in birds22. What makes the Midas cichlid case so interesting and seemingly unique until now is that it seems to be such a surprisingly complex and intricate molecular process that leads to the loss of pigment cells in adult cichlids. After all, many other mutations could contribute to the process of turning Midas cichlids gold, including melanophore markers, such as tyrosinase, oca2, or members of the Agouti/Mc1r pathway that have been shown to abolish melanin production in cichlids23,24,25. Ultimately, the Midas cichlid gold locus will help address longstanding questions on the proximate and ultimate causes of stable (color) polymorphisms. Future research will show why this insertion and the mechanism by which it acts to trigger pigment loss at a specific ontogenetic phase has been maintained by selection in so many Midas cichlid populations. Our preliminary data on the ecological relevance of this, to the human eye at least, much more conspicuous Gold morph in this color polymorphism might suggest that the gold morphs are maintained at low frequencies in all populations because predators do not seem to prefer them26.

In summary, we identified what might constitute the genetic basis of the trans-species gold/dark polymorphism that is eponymous of Midas cichlids: a transposon-derived structural variant that likely causes de-regulation of a previously undescribed gene goldentouch via the formation of a huge cruciform secondary structure. We anticipate that the knowledge of the likely causal mutation will not only greatly contribute to deciphering the molecular mechanism of this fascinating case of morphological color change but also lead to exciting avenues of research on the interaction between genomic structural features, assortative mating27,28, selection regimes29, and ecological settings26 that lead to the stable persistence of intraspecific polymorphisms3,4.

Read more here: Source link