Genetic and Developmental Divergence in the Neural Crest Program between Cichlid Fish Species

Abstract

Neural crest (NC) is a vertebrate-specific embryonic progenitor cell population at the basis of important vertebrate features such as the craniofacial skeleton and pigmentation patterns. Despite the wide-ranging variation of NC-derived traits across vertebrates, the contribution of NC to species diversification remains underexplored. Here, leveraging the adaptive diversity of African Great Lakes' cichlid species, we combined comparative transcriptomics and population genomics to investigate the evolution of the NC genetic program in the context of their morphological divergence. Our analysis revealed substantial differences in transcriptional landscapes across somitogenesis, an embryonic period coinciding with NC development and migration. This included dozens of genes with described functions in the vertebrate NC gene regulatory network, several of which showed signatures of positive selection. Among candidates showing between-species expression divergence, we focused on teleost-specific paralogs of the NC-specifier sox10 (sox10a and sox10b) as prime candidates to influence NC development. These genes, expressed in NC cells, displayed remarkable spatio-temporal variation in cichlids, suggesting their contribution to interspecific morphological differences, such as craniofacial structures and pigmentation. Finally, through CRISPR/Cas9 mutagenesis, we demonstrated the functional divergence between cichlid sox10 paralogs, with the acquisition of a novel skeletogenic function by sox10a. When compared with teleost models zebrafish and medaka, our findings reveal that sox10 duplication, although retained in most teleost lineages, had variable functional fates across their phylogeny. Altogether, our study suggests that NC-related processes—particularly those controlled by sox10s—are involved in generating morphological diversification between species and lays the groundwork for further investigations into the mechanisms underpinning vertebrate NC diversification.

evolutionary developmental biology, gene duplication, neural crest cell biology, cichlid fishes

Introduction

The remarkable diversity and complexity of craniofacial structures, pigmentation patterns, and social behaviors within vertebrates is a testament to their outstanding capacity to adapt and exploit a wide range of ecological niches. Much of this phenotypic diversity is intimately connected with the emergence of the neural crest (NC; Gans and Northcutt 1983; Donoghue et al. 2008). This embryonic multipotent cell population arises from the dorsal portions of the neural tube and then migrates extensively to finally differentiate into a broad range of cell types and tissues, including neurons and glia, pigment cells, craniofacial cartilage and bone, among others (Bronner and Simões-Costa 2016; Brandon et al. 2023). These diverse cell lineages later assemble to form complex pigmentation patterns in fish, amphibians, and birds, as well as divergent head structures, such as fish jaws, bird beaks, or mammalian horns (Eames and Schneider 2005; Jheon and Schneider 2009; Nasoori 2020; Elkin et al. 2023).

NC has been primarily studied in the context of its origin, development and function, including developmental disorders involving its derivatives (neurocristopathies; Bolande 1997). Studies in model organisms have revealed that the gene regulatory networks (GRNs) and developmental processes governing NC specification, migration, and differentiation are highly conserved across distantly related species (Simões-Costa and Bronner 2015; Bronner and Simões-Costa 2016). This remarkable macroevolutionary conservation of the NC program raises key questions about its evolvability and its potential contribution to the origins of vertebrate diversity. Surprisingly, the role that NC cells may play in the evolution of species-specific traits (i.e. at the population or species level) remains largely unexplored (Abzhanov et al. 2004; Donoghue et al. 2008; Powder et al. 2014; Kratochwil et al. 2015; Brandon et al. 2023). This is despite the rapid and extensive diversification of NC-derived structures—a hallmark of adaptive radiations of multiple vertebrate clades. Striking examples include the diversification of cranial shapes of Anolis lizards, beak morphologies in Darwin's finches and, perhaps most spectacularly, the craniofacial skeletons and color patterns of cichlid fish radiations in the Great African Rift Lakes (Abzhanov et al. 2006; Albertson and Kocher 2006; Sanger et al. 2012; Santos et al. 2023).

Here, we investigated the molecular evolution of NC-related phenotypic diversity in the spectacular and recent (700 to 800 thousand years ago) Lake Malawi cichlid fishes radiation, which includes over 600 to 800 species displaying some of the lowest sequence divergence observed in vertebrates (Turner 2007; Brawand et al. 2014; Meyer et al. 2017; Malinsky et al. 2018; Vernaz et al. 2021) (also see Materials and Methods). To this end, we first examined variation in NC genetic and developmental programs between two genetically closely related, yet eco-morphologically divergent cichlid species, namely the generalist Astatotilapia calliptera “Mbaka” and the pelagic piscivore Rhamphochromis sp. “chilingali.” Both species belong to the Lake Malawi cichlid radiation and are characterized by distinct craniofacial morphologies, body pigmentation, ecologies and diets, representing eco-morphological extremes within this radiation (Turner 2007; Salzburger 2018; Edgley and Genner 2019; Santos et al. 2023) (Fig. 1a and b and supplementary fig. S1, Supplementary Material online). Our previous work identified variation in pigmentation and craniofacial shapes at the earliest stages of their overt appearance at posthatching stages (Marconi et al. 2023). Considering the direct mode of development in cichlids (i.e. without larval stage and metamorphosis, unlike zebrafish; Woltering et al. 2018; Marconi et al. 2023), the variation in these NC-derived traits likely originates from differences in early embryogenesis (Marconi et al. 2023; supplementary fig. S1a to d, Supplementary Material online).

Fig. 1.

Divergence in whole-embryo transcriptomic trajectories during early embryogenesis and NC development between two eco-morphologically distinct Lake Malawi cichlids. a) Geographical map of Lake Malawi/Nyasa with the estimated number of cichlid species and the estimated age of the radiation. b) Cichlid species part of this study exhibit distinct NC-derived craniofacial morphologies and pigmentation patterns (to scale). AC has intermediate phenotype of a generalist feeder with variable melanic patches comprising features of both bars and stripes on a uniform background, whereas RC has a flattened head and elongated jaws, typical of a piscivore, with light coloration and dark horizontal stripes. c) AC and RC exhibit different total somite numbers upon completion of somitogenesis (30 in AC, 38 in RC). d) The stages of cichlid somitogenesis (expressed as ss) examined in this study and collected for RNA sequencing range from the early stages of NC specification (4ss) through migratory NC (10 to 12ss onwards) and its differentiation during late somitogenesis stages, concluding at 30 and 38ss in AC and RC, respectively (n = 3 biological replicates per ss). Overview of NC development and migration (in blue) is based on the in situ data presented in this study. e) PCA of whole transcriptome samples reveals significant ontogenic (PC1) and species-specific (PC2) clustering. Each data point corresponds to a single replicate embryo. A, anterior; AC, Astatotilapia calliptera “Mbaka”; dpf, days postfertilization; le, lens; oc, optic cup; op, optic primordium; ov, otic vesicle; P, posterior; RC, Rhamphochromis sp. “chilingali”; s, somites; ss, somite stage; st, stage; b, tri-partite brain; V, V-shaped somites. Map a) modified from d-maps.com.

To test this hypothesis, we focused on the interspecific comparison at earlier embryonic stages concomitant with NC cell specification, migration, and onset of their differentiation (Rocha et al. 2020; Fig. 1c and d). Using whole-transcriptome time-series sequencing data, we uncovered substantial variation in coding and noncoding gene expression throughout NC development between the two species (Fig. 1d and e and supplementary fig. S1e to i, Supplementary Material online) and included divergence in expression levels and temporal trajectories of dozens of genes with key functions within the teleost NC-GRN (Rocha et al. 2020). Moreover, we show that several of these crucial genes are also associated with signatures of divergent positive selection between species, potentially contributing to species-specific phenotypes. We then focused on two SRY-box transcription factor 10 (sox10) paralogs of a key NC specifier, namely sox10a and sox10b, and showed that they both arose during teleost-specific whole genome duplication (WGD) and exhibit prevalent interspecific expression variation throughout NC development. These results suggest potential contribution of sox10 paralogs, and NC development more broadly, to species differences. Finally, we provide experimental evidence that sox10a function is essential for craniofacial skeletal development, indicating a novel role in cichlids that has not been described in any teleost to date. Taken together, our study reveals that sox10 paralogs followed divergent functional evolution across the teleost phylogeny, including gene loss (zebrafish), subfunctionalization (medaka) and neofunctionalization (cichlids). We propose that the expansion of genetic toolkit associated with NC development during genome duplication, subsequent lineage-specific divergence of paralogous genes, including the acquisition of novel functions, and regulatory and transcriptomic evolution in cichlids, may have collectively contributed to the extensive morphological diversification in this clade. Our results highlight cichlids as a unique teleost system to investigate the developmental and genetic underpinnings of adaptive phenotypic evolution.

Materials and Methods

Animal Husbandry and Embryo Culture

Breeding stocks of A. calliptera “Mbaka” (AC) and Rhamphochromis sp. “chilingali” (RC) were maintained under standardized conditions as previously described in (Marconi et al. 2023). Eggs used for RNA extractions and HCR in situ hybridization experiments were collected from mouthbrooding females immediately after fertilization and then reared individually in 1 mg/L of methylene blue (Sigma Aldrich) in water in 6-well plates (ThermoFisher Scientific) placed on an orbital shaker moving at slow speed at 27 °C until needed. All experiments were conducted in compliance with the UK Home Office regulations.

Whole Embryo Bulk RNA Sequencing

Sample Acquisition

For each examined species, all samples were collected from the same egg clutch. Sampling covered the entire period of somitogenesis, which coincides with NC development, with samples taken at 3h intervals. At each time point, at least four embryos were dissected and placed individually into either 250 μL of prechilled Trizol (Ambion) and stored at −80 °C until RNA extraction (at least overnight) or into 1 mL of 4% PFA in 1 × PBS for overnight fixation at 4 °C. Embryos preserved in 4% PFA were later rinsed twice in 1 × PBS and stained with 10 nM DAPI in 70% glycerol in 1 × PBS overnight at 4 °C, protected from light. Following a wash in 1 × PBS (10 min/wash, once), the embryos were mounted on microscopy slides (ThermoFisher) with Fluoromount G (Southern Biotech) and imaged with an Olympus FV3000 confocal microscope to confirm the developmental age (ss) of each sampled cohort.

RNA Extraction

All procedures were conducted on ice, unless otherwise specified. Samples stored in Trizol were thawed from −80 °C. For each sample, 100 mg of 0.1 mm zirconia/silica beads (Stratech) were added before homogenization using a TissueLyser II (Qiagen) for 120 s at 30 Hz. The samples were then topped up to 1 mL with chilled Trizol and allowed to rest for 5 min. Next, 200 μL of chloroform (ThermoFisher Scientific) was added and the samples were vigorously shaken for 15 s, briefly vortexed and incubated at room temperature for 15 min. The samples were then centrifuged at 300 × g for 20 min at 4 °C. The supernatant was carefully transferred to a fresh tube and further processed using the Direct-Zol RNA Microprep Kit (Zymo) according to the manufacturer's instructions. The quality and quantity of the extracted total RNA were assessed using Qubit (RNA HS assay, Agilent) and Tapestation (Agilent). Total RNA extracted from each embryo was submitted individually for sequencing, with quantities ranging from 135 ng to 1.3 μg per sample. All sequenced samples had eRIN values above 9.3.

NGS Library Preparation

All libraries were prepared, quality-controlled and sequenced by Novogene Corporation (China) using the Illumina NovaSeq 6,000 platform to generate paired-end reads of 150 base pairs (bp). On average, 32.49 ± 2.5 million paired-end 150 bp reads were generated per sample (supplementary table S1, Supplementary Material online).

Adapter Trimming and Quality Filtering

The adapter sequences in reads were removed, and low-quality sequences (Phred < 20) were filtered out with TrimGalore (v0.6.6).

Mapping of RNAseq Reads to Reference Genome

All RNAseq reads were mapped to the A. calliptera genome assembly (fAstCal1.2 in Ensembl 105), which contains 27,018 annotated genes, similar to the other Malawi genome, Maylandia zebra (27,187 annotated genes). Mapping was performed using STAR v.2.7.1a (Dobin et al. 2013), following established protocols for cichlids (El Taher et al. 2021), and only uniquely mapped reads were used. The very low sequence divergence among Malawi cichlid species enables the use of a single reference genome. Sequence divergence across 73 Malawi cichlid species was shown to be one of the lowest seen in vertebrates and averages ∼2.0 × 10⁻³ per bp, with the AC and RC populations showing similarly low sequence divergence (2.07 × 10⁻³ per bp) (Malinsky et al. 2018, Svardal et al. 2021). The mapping rates for uniquely mapped reads were high, yielding mean mapping rates of 90.2 ± 0.8% and 89.1 ± 0.7% for AC and RC, respectively (mean ± standard deviation; t-test, P = 4.94 × 10⁻⁵; see supplementary fig. S1e and supplementary table S1, Supplementary Material online).

Gene Expression Quantification

The number of reads mapped to each gene in the reference genome was counted in STAR using the built-in HTSeq-count option (Anders et al. 2015). Gene counts were normalized using the median of ratios method in DESeq2 (Love et al. 2014) (v1.34.0). A PCA was applied to reduce the dimensionality of the dataset using the R command prcomp (R 4.2.0).

DE Analysis and Gene Annotation

DE analysis was performed on a gene count matrix using DESeq2 (Love et al. 2014) (v.1.34.0) using pairwise comparisons of embryos of the two species matched by their total somite count, excluding 37ss RC samples as no equivalent ss exists in AC. Genes with mRNA counts < 10 per sample in each species were filtered out prior to analysis and technical replicates for each stage were collapsed. Heatmaps of scaled gene expression (Z-score was calculated for each selected gene across all samples [of both species across all ss] using mean DESeq2-normalized gene count) were generated using pheatmap (v.1.0.12) and unbiased hierarchical clustering of gene expression patterns was performed using the complete linkage method. In total, seven clusters effectively categorized the variation into distinct expression pattern groups. Scaled gene expression for each gene per cluster was then plotted using ggplot2 (v.3.3.6).

Genes were annotated using the reference genome A. calliptera (fAst.cal 1.2) with biomaRt package (Durinck et al. 2009) (v.2.54.0). Significantly DEG for each pairwise comparison using ss index were identified by filtering for log2 Fold Change above the absolute value of 0.585 (i.e. ≥ 1.5-fold difference in expression) and adjusted P-value < 0.05.

GO Analysis

GO annotation and functional enrichment analysis were carried out using gProfiler2 (Kolberg et al. 2020). The full extent of the GO annotation, including the most specific GO terms for each gene or gene product, was used rather than the commonly used GOslim annotation, which comprises only a subset of the terms belonging to each parent domain, thus reflecting the broader biological categories. To focus the candidate search on genes involved in the NC development, a dataset comprising all DE genes with GO terms broadly associated with the developmental program of the NC (including GO terms related to its development, migration, and differentiation), development of pigmentation as well as craniofacial complex and cartilage was compiled (see supplementary table S2, Supplementary Material online for GO terms used for filtering and supplementary table S3, Supplementary Material online for complete candidate gene list). The identified candidates (supplementary table S3, Supplementary Material online) were further categorized into functional tiers within the NC-GRN (Fig. 2a) based on their GO annotation (supplementary table S4, Supplementary Material online) as well as the functional information from zebrafish summarized in Rocha et al. (2020).

Fig. 2.

Comparative characterization of cichlid transcriptomes during somitogenesis and NC development. a) Heatmap of scaled expression for all differentially expressed genes (DEGs, n = 12,611) across pairwise comparisons of embryos from the two species, matched by their total somite count (see Materials and Methods) reveals seven distinct gene expression clusters. b) GO terms overrepresented among DEGs identified in pairwise and pooled (all) comparisons indicates variation in multiple embryonic processes, including NC development and migration (g:SCS significance threshold, P < 0.05). For clarity, only the top-tier parent term for metabolism-related functions (i.e. “small molecular metabolic process”) is presented. Full GO annotation is presented in supplementary table S8, Supplementary Material online. c) Genome-wide scans between Astatotilapia and Rhamphochromis populations identify 551 significant SNPs (in pink) with elevated local haplotype homozygosity, measured using between-population extended haplotype homozygosity (xpEHH; see also supplementary fig. S2b, Supplementary Material online). d) Significant SNPs in close proximity (see Materials and Methods) were grouped into 154 islands, which are associated with 74 DEGs (either with gene bodies or in proximal intergenic regions). e) NC-DEGs within cluster six show significantly elevated xpEHH values, in line with signatures of positive selection (see full gene list in supplementary fig. S2d and supplementary table S5, Supplementary Material online). The number of genes in each cluster for each GO category is shown below the graphs. DEG, differentially expressed genes; nDEG, nondifferentially expressed; NC-DEG, DEGs associated with NC development, xpEHH—extended haplotype homozygosity.

TEs and Repeat Expression Quantification

TEs and repeats were predicted using RepeatModeler (v. 2.0.2 with LTRStruct parameter) and were then annotated using the RepeatModeler custom library in the A. calliptera genome using RepeatMasker (v.4.1.4; http://repeatmasker.org/). Only transposon elements (TE) were further analyzed (simple repeats, tRNA, rRNA, scRNA, and satellites were excluded from subsequent analyses). To quantify TE gene expression, TEcount (v.1.0.1) from the TEtranscript package (Jin et al. 2015) was used following STAR mapping with the following parameters–chimSegmentMin 10–winAnchorMultimapNmax 200–outFilterMultimapNmax 100. Finally, a DESeq2-normalized gene count matrix for 1,609 different transcribed TEs was generated and PCA (centered and scaled) was produced with R (prcomp).

Population Genomics

Detection of Sites under Positive Selection

To identify signatures of positive selections, we performed genome-wide scans to detect regions with unusually high local haplotype homozygosity using extended haplotype homozygosity (XP-EHH) analysis between the two populations with REHH (v.3.2.2; Gautier and Vitalis 2012; Bourgeois and Warren 2021). This approach follows recent established protocols used in cichlids (Kautt et al. 2020), other vertebrates (Cong et al. 2022; Schield et al. 2022) including humans, and invertebrates (Smith et al. 2022). We downloaded and analyzed available VCF files containing biallelic SNPs from 43 samples from the Rhamphochromis genus (23.6 ± 8.3× sequencing depth, mean ± SD) and 45 randomly chosen A. calliptera “Masoko Benthic” (16.02 ± 1.63×) (supplementary fig. S2a, Supplementary Material online; see Materials and Methods in https://github.com/tplinderoth/cichlids/tree/master/callset and genome data generated and downloaded from ref. Munby et al. 2021) to ensure comparable sample sizes per population. Access to sequencing data for R. chilingali (RC) was limited to a few individuals (included in the analysis), we included other Rhamphochromis species to maintain statistical power. Consequently, observed genetic differences between AC and RC populations might be shared across the Rhamphochromis genus rather than being unique to RC. All samples were all wild-caught and sequenced individually (not pooled). Following the approach by (Ravinet et al. 2018), we calculated between-population extended haplotype homozygosity (xpEHH) using phased variants with a minor allele frequency > 0.05. In total, 551 significant individual sites (SNPs) were identified across all chromosomes (log[P value] ≥ 3), suggesting potential positive selection. Significant sites found within 50 kbp (±50 kbp) were grouped into “islands” using bedTools (2.29.2), resulting in 154 regions, which comprised between 1 and 39 significant SNPs (3.6 SNPs on average) and varied in size between 1 bp and 69 Mb (mean = 4.2 kb). To minimize false positives, we focused our analysis on significant peaks within gene bodies or, if located outside gene bodies, within the close proximity (±200 kbp away from the nearest gene using bedTools), which may approximate cis-regulatory regions. In total, 184 genes were linked to putative islands of selection (1 to 5 genes per island, mean = 1.2), of which 73 were DEGs (supplementary table S5, Supplementary Material online; 48, 5, and 20 islands located in gene bodies, promoter and intergenic regions, respectively). To identify enrichment in sites under potential positive selection within gene expression clusters, permutation tests were performed between the observed distribution of xpEHH P-values across candidate genes (25 kbp regions upstream of TSS) for each gene expression cluster and the expected distribution (over chance). Expected values were calculated by randomly selecting xpEHH values (median values) across size-matched genic windows (1000 × iterations; supplementary fig. S2d, Supplementary Material online). GO enrichment analysis was performed on genes associated with xpEHH peaks using G:Profiler.

Phylogenetic Reconstructions

The coding sequences of soxE gene family members were retrieved from Ensembl (108) and directly from genome assemblies from Parey et al. (2023) for inclusion in the phylogenetic analyses (supplementary table S6, Supplementary Material online). Multiple sequence alignments of soxE gene family members were constructed using Clustal Omega (Sievers and Higgins 2018). Next, the alignments were trimmed with TrimAl (Capella-Gutiérrez et al. 2009) and used to infer evolutionary relationships with the maximum likelihood method in IQ-TREE v.1.6.12 (Nguyen et al. 2015). In-built ModelFinder (Kalyaanamoorthy et al. 2017) was used to infer the best-fit substitution model (TN+F+I+G4) based on the Bayesian information criterion. The branch supports for maximum likelihood analyses were obtained using the ultrafast boot-strap (Kalyaanamoorthy et al. 2017) with 1,000 replicates. Phylogenetic tree shown in Fig. 3 was visualized using iTOL v.6 (Letunic and Bork 2019) and only nonsignificant bootstrap values are shown (≤75%).

Fig. 3.

Sox10 paralogs were duplicated in teleosts and are differentially expressed during embryonic development in cichlids. a) The topology of maximum likelihood phylogeny of sox10 paralog coding sequences across vertebrates confirms sox10 duplication occurred at the root of teleosts, followed by sox10a loss in zebrafish Danio rerio (order: Cypriniformes) and cavefish Astyanax mexicanus (order: Characiformes), both belonging to superorder Ostariophysi (green arrow). Phylogeny constructed with IQ-TREE (Nguyen et al. 2015). All significant bootstrap values, apart from the ones shown (<75). b) The expression trajectories of sox10 paralogs across NC development in Astatotilapia calliptera “Mbaka” (AC) and Rhamphochromis sp. “chillingali” (RC). Shaded bands indicate 95% confidence intervals. Significant differences in normalized gene counts between sox10 paralogs were observed at 10ss in both examined species (two-way ANOVA and Tukey HSD, P < 0.01). c) Fold changes in expression levels of sox10 paralogs between species. Only significant comparisons are shown (P-adj < 0.05). cNCC, cranial neural crest cell; NCC, neural crest cell; n.s., nonsignificant.

Genome Editing

gRNA Design and Synthesis

Targets for CRISPR/Cas9 editing were selected with the CHOPCHOP software online (http://chopchop.cbu.uib.no/) (Labun et al. 2019) using the Astatotilapia burtoni genome (AstBur1.0) as a reference. Sequence similarity searches with BLAST against the A. calliptera genome (AstCal1.2, Ensembl 108) were performed to confirm homology and test for off-target effects. Target sequence specificity was further confirmed using Geneious Prime (Dotmatics) software. Two sgRNAs targeting exon 1 were designed for sox10b (5′-CTCGTCGTCGGATTTGACGG-3′ and 5′-CGCGGATTCCCGCGGGGAA-3′) and sox10a (5′-CGGTCAGTCAGGTGCTGGACGGG-3′ and 5′-TCGTTTCCCGATCGGCATAA-3′), respectively, and purchased from Integrated DNA Technologies (ITD) as Alt-R CRISPR-Cas9 sgRNA (2 nmol).

Microinjection

Single-cell embryos of A. calliptera “Mbaka” were injected and maintained following the protocol described in Clark et al. (2022). We were unable to perform microinjections in a similar manner in Rhamphochromis due to their prolonged and unpredictable breeding behavior, rendering it technically infeasible.

Embryo Imaging

Injected and control embryos were imaged daily until 12 days postfertilization using a Leica M205 stereoscope with a DFC7000T camera under reflected light darkfield. All specimens were positioned in 1% low melting point agarose (Promega) and anesthetized with 0.02% MS-222 (Sigma-Aldrich) if required to immobilize during imaging.

Genotyping

Tissue samples were taken from specimens sacrificed by overdose of 0.5% MS-222 (Sigma-Aldrich) to extract genomic DNA using PCRBIO Rapid Extract Lysis Kit (PCRBiosystems). Fragments of 190 to 420 bp surrounding predicted deletion sites were amplified using PCRBIO HS Taq Mix Red (PCRBiosystems) with an annealing temperature of 56 °C using following primer pairs: sox10b 5′-CTGTCACCGGGTCATTCCTC-3′ and 5′-GCGTTCATTGGCCTCTTCAC-3′; sox10a 5′-ATGGTCACTCACTGTCACCG-3′ and 5′-CCTCCTCGATGAATGGCCTC-3′. Amplicons were purified with QIAquick PCR Purification Kit (Qiagen) before Sanger sequencing. All protocols were conducted following manufacturer’s instructions. Sequence analysis to infer CRISPR edit sites was performed using the Synthego ICE CRISPR analysis tool (https://ice.synthego.com/).

Cartilage Preparations

Embryos were stained for cartilage following the protocol of Marconi et al. (2023) with following modifications: (i) all specimens were bleached to remove melanophore pigmentation using a solution of 0.05% hydrogen peroxide (Sigma) and 0.05% formamide (ThermoFisher) for 30 to 45 min under light and (ii) samples were cleared using first 50% then 70% glycerol:water solutions until complete sinking. Specimens were stored in 70% glycerol until imaging in 80% glycerol using a Leica M205 stereoscope with a DFC7000T camera under reflected light.

Whole Mount in situ Hybridization by HCR

Reagents

The HCR probes and hairpin sets were ordered from Molecular Instruments, whereas all required buffers were made following the instructions provided by the manufacturer. The HCR probe sets (14 to 20 pairs per gene) were designed using target gene template sequences retrieved from A. calliptera genome assembly (fAstCal1.2, Ensembl 108) (supplementary table S7, Supplementary Material online). Each probe set was designed by the manufacturer to target transcript regions common to all splicing isoforms while minimizing off-target effects.

Due to a very low genetic variation in the coding sequences between study species, we used the same probe sets per each target gene for both cichlid taxa examined. Probe specificity was verified by BLAST searches against the A. calliptera genome available on Ensembl (AstCal 1.2) as well as against unpublished Rhamphochromis sp. “chilingali” assembly.

Protocol Overview

Dissected embryos were fixed overnight at 4 °C in 4% PFA in 1 × PBS. The following day, they were rinsed twice in 1× PBST (1× PBS+ 0.01% Tween-20) without incubation and washed in 1 × PBST (10 min/wash, twice) before a stepwise dehydration to 100% MetOH in prechilled solutions of 25%, 50%, and 75% MetOH:PBST (10 min/wash, at 4 °C) and stored at −20 °C until further analyses.

The embryos were then rehydrated from 100% MetOH to 1 × PBST in reciprocal series (5 min/wash, at 4 °C), followed by washes in 1 × PBST at room temperature (5 min/wash, twice).

The mRNA in situ hybridization by chain reaction (HCR) was carried out according to the protocol of Andrews et al. (2021) for whole mount amphioxus embryos with the following modifications. (i) 2 pmol of each probe mixture (1 μL of 2 μM stock) per 100 μL of probe hybridization buffer were used and (ii) 60 pmol of each fluorescently labeled hairpin (i.e. 2 μL of 3 μM stock) were applied per 100 μL of amplification buffer. Finally, the embryos were stained with 10 nM DAPI in 70% glycerol in 1 × PBS overnight at 4 °C protected from light, washed in 1 × SSCT (10 min/wash, twice) before mounting with Fluoromount G (Southern Biotech) on glass bottom dishes (Cellvis) without coverslip or microscopy slides (ThermoFisher Scientific) with #1.5 coverslips (Corning), depending on the size of the specimen. To prevent embryos from getting squashed when mounting on slides, thin strips of electrical tape were used as bridges to create space between the slide and a coverslip. Clear nail varnish was used to seal the edges of the slide and all samples were cured overnight at room temperature protected from light before imaging.

All in situ hybridization experiments were performed with multiple specimens from different clutches (at least 3 individuals per clutch, repeated at least once with specimens from alternative clutches) to fully characterize the expression patterns.

Confocal Microscopy

Imaging of dissected and stained embryos was carried out with an inverted confocal microscope Olympus FV3000 at the Imaging facility of the Department of Zoology, University of Cambridge. As the fluorescence intensity levels were only compared as relative signals within each sample (i.e. embryo), imaging was performed using optimal laser power and emission wavelength for each sample. Sequential acquisition mode was used to minimize the signal crosstalk across channels and all images were acquired at 1,024 × 768 resolution and 12-bit depth.

Image Processing for Figure Presentation

Confocal micrographs were stitched using the Olympus FV3000 software and processed with Fiji (Schindelin et al. 2012) to produce optical sections, collapse z-stacks, and adjust image brightness and contrast where necessary following guidelines by Schmied and Jambor (2021). The presented transverse optical sections are maximum projection across five adjacent slices. Auto-fluorescence and background noise in each image was removed by subtracting the average pixel intensity measured for each channel in regions of the embryo where no fluorescent signal was observed. Any remaining overexposed pixels were removed using the “Remove outliers” (radius = 0.5 pixel, threshold = 50) function in Fiji. Images presented as figures were smoothened by applying a Gaussian filter with σ = 0.5. Image look-up tables (LUTs) were taken from the BIOP plugin (https://imagej.net/plugins/biop-lookup-tables) and “ChrisLUTs” (Christophe Leterrier and Scott Harden; github.com/cleterrier/ChrisLUTs) package for Fiji/ImageJ. The images were processed for any background imperfections and assembled into figures in Adobe Photoshop 2023.

All plots and statistical analyses were made using R version 4.1.1. The MetBrewer package (https://github.com/BlakeRMills/MetBrewer) was used for color palettes throughout this work.

Results

Genes Involved in NC Development Show Species-Specific Shifts in Transcriptomic Trajectories

Our previous work has shown that phenotypic divergence between cichlid species in NC-derived craniofacial skeleton and body pigmentation is first observed at the appearance of differentiated cartilage and pigment-bearing cells at early posthatching stages, respectively (Marconi et al. 2023; supplementary fig. S1c and d, Supplementary Material online). Given that cichlids do not undergo metamorphosis and develop directly from embryo to adult (Woltering et al. 2018; Marconi et al. 2023), we hypothesize that variation in processes of NC development occurring early in ontogeny might constitute an important contributor to the adult morphological divergence between these species.

To first investigate the potential role of gene expression variation in driving divergent NC-derived phenotypes, we performed comparative transcriptome profiling (RNAseq) of whole embryos of Astatotilapia and Rhamphochromis across somitogenesis. This period of embryonic development coincides temporally with NC development (Rocha et al. 2020; Fig. 1c). In total, 32.49 ± 2.5 million paired-end 150 bp-long reads were generated for each sample (three biological replicates per somite stage, ss), and then aligned against the A. calliptera reference genome to quantify gene expression, yielding high mapping rates of 89.7 ± 1.0% across all samples (supplementary fig. S1e and supplementary table S1, Supplementary Material online and see Materials and Methods). Principal component analysis (PCA) of protein-coding transcriptomes revealed that gene expression in these species is primarily dictated by ontogeny (i.e. their developmental age in ss; PC1, 24.73%), followed by species (PC2, 12.95%; Fig. 1e and supplementary fig. S1f, Supplementary Material online). Conversely, the expression of noncoding transcripts and transcribed transposable elements (TEs) is primarily clustered by species and then by ontogeny (supplementary fig. S1g to i, Supplementary Material online), in line with previous reports of higher evolutionary rates associated with expression dynamics of noncoding genes in cichlids (El Taher et al. 2021).

We then performed differential gene expression (DE) analysis across somitogenesis stages to identify gene candidates showing species- and time-specific transcriptional patterns. In total, 12,611 differentially expressed genes (DEGs; P < 0.05) with >1.5-fold expression difference in at least one pairwise comparison of ss-matched embryos were identified (Fig. 2a, Materials and Methods and Supplementary material online). Of these, 14.7% (n = 1,857) lacked assigned names in the current assembly A. calliptera genome (Ensembl 108), likely representing genes without zebrafish orthologs or novel genes, and were excluded from downstream analyses. The remaining DEGs were then classified into seven distinct clusters based on unbiased grouping of their expression patterns (Fig. 2a and supplementary fig. S3a and b, Supplementary Material online), with gene functions significantly enriched for distinct biological processes in each cluster, including NC-related processes in clusters 4 and 7. While four of these clusters displayed consistent high or low gene expression in one species across all ss (clusters 2, 3, 5, and 6; accounting for 51.4% of all DEGs), the other clusters (1, 4, and 7) showed species-specific temporal shifts in gene expression, possibly linked to the temporal differences (heterochronies) during somitogenesis between these species (Marconi et al. 2023). We next examined the enriched Gene Ontology (GO) terms associated with the identified DEGs in discrete pairwise comparisons as well as the pooled set of all DEGs (Fig. 2b, supplementary table S8, Supplementary Material online). Overall, the highest number of DEGs were linked to functions related to metabolism, followed by transcription and signaling pathways. Genes involved in many of the developmental processes concomitant with our time series (e.g. brain and heart development, somitogenesis) were also significantly overrepresented, particularly at 15ss and across all comparisons. Notably, these two comparisons were also enriched in terms involving NC development and NC migration.

Altogether, the considerable divergence in expression dynamics between species over developmental time (e.g. temporal shifts) may imply variation in multiple developmental processes during embryogenesis, including those involving NC cells. Combined with the diverse repertoire of NC derivatives, these results implicate the potential role of this cell population in the divergence of species-specific traits.

Signatures of Positive Selection are Associated with NC-Related Genes

We then sought to assess if DEGs, and in particular NC-related genes, were potentially diverging between species by conducting genome-wide scans for regions under positive selection. Using extended haplotype homozygosity (xp-EHH, Gautier and Vitalis 2012), scans between wild-caught Astatotilapia and Rhamphochromis populations (43 to 45 whole genomes per species; sequencing data from Munby et al. unpublished), we identified 154 regions showing significant signatures of positive selection (Fig. 2c, supplementary fig. S2b and supplementary table S5, Supplementary Material online, Materials and Methods). Altogether, 74 DEGs were located near or within these putative islands of selection (Fig. 2d, supplementary fig. S2b, c and supplementary table S5, Supplementary Material online; see Materials and Methods), with functions ranging from cell differentiation, signal transduction to metabolic pathways (supplementary fig. S2e, Supplementary Material online), which might potentially contribute to aspects of their eco-morphological divergence. Additionally, several DEGs were unannotated novel genes (supplementary fig. S2f, Supplementary Material online). The similar expression patterns of these unannotated genes compared with known genes (supplementary fig. S2f, Supplementary Material online) warrant further research to investigate their roles in divergent developmental programs between AC and RC, temporally overlapping with, and possibly involving, NC.

Among the annotated DEGs located in the most extreme outlier regions were fgf8a (cluster 7), kcnk18 (cluster 1), and tspan37 (cluster 5), which exhibited significant expression differences between the two species during the early and mid-phases of somitogenesis (supplementary fig. S2b and c, Supplementary Material online). tspan37 encodes an integral membrane protein involved in cellular signaling, while kcnk18 encodes a potassium channel expressed in the brain and eye of hatchling and larval zebrafish (Rahm et al. 2014). To test whether these genes were involved in NC development and characterize their expression patterns, we performed in situ Hybridization Chain Reaction (HCR) (Choi et al. 2018). Expression of kcnk18 and tspan37 was not detected in whole mount embryos at the stages of differential expression, perhaps due to overall low expression levels (supplementary fig. S2c, Supplementary Material online). Furthermore, fgf8a, albeit known for its diverse roles during embryogenesis, including in chondrogenesis of the NC-derived cranial skeleton (Gebuijs et al. 2019), showed expression only in the developing brain and posterior-most notochord during somitogenesis in the examined cichlids (supplementary fig. S4, Supplementary Material online), consistent with findings in zebrafish. This highly localized expression pattern suggests an unlikely role in craniofacial development and patterning at this stage. Altogether, the top divergent outlier genes are not likely to play a role in NC-derived divergence between these two cichlid species. However, NC-related DEGs (identified based on their GO annotation involving NC, supplementary tables S2 to S3, Supplementary Material online, Materials and Methods) belonging to cluster 6 and showing overall lower expression in Rhamphochromis throughout somitogenesis (supplementary fig. S3b, Supplementary Material online) displayed significant enrichment for sites under potential positive selection (Fig. 2e and supplementary fig. S2d, Supplementary Material online). These included the transcription factor pax3a, already implicated in cichlid interspecific pigment pattern variation (Albertson et al. 2014), and the cellular nucleic acid-binding cnbpa, involved in craniofacial development in fish (Weiner et al. 2007), among others. Collectively, these findings highlight a strong association between sequence and transcriptional differences between the two species, with several DEGs with functions related to NC processes showing enrichment for sites under potential positive selection.

Variation across the NC-GRN is Particularly Associated with NC Cell Migration

As our comparative transcriptional and selection analyses highlighted divergence in the genetic program of NC development between species, we next sought to identify genes, processes and stages associated with early NC ontogeny potentially implicated in the evolution of morphological diversity. To this end, we examined only DEGs with known functions (based on their GO annotation, supplementary tables S2, S3 and S8, Supplementary Material online and Materials and Methods) in NC development and differentiation and its two extensively diversified derivatives, namely pigmentation and craniofacial skeleton (Fig. 4a).

Fig. 4.

Prevalent transcriptomic variation between cichlids exists across the NC genetic program. a) Identified DEGs belong to different tiers of the teleost NC-GRN (Rocha et al. 2020), from specification to migration and differentiation, including development of NC-derived pigmentation and craniofacial skeleton. b) Distribution of identified NC-DEGs per function along the RNAseq time-course. c) Fold changes in expression levels of candidate genes involved in NC development between cichlids. Genes involved in NCC migration highlighted in red and white cells correspond to nonsignificant differences between species. d) In situ HCR image showing prdm1a expression at 4ss, representative of n ≥ 2 per species. Anterior to the top of the figure. AC, Astatotilapia calliptera “Mbaka”; ANS, autonomic nervous system; NCC, neural crest cells; RC, Rhamphochromis sp. “chilingali”; SNS, sympathetic nervous system; ss, somite stage. Scale bar = 100 μm.

The identified NC-DEGs belonged to all tiers of the NC-GRN (Sauka-Spengler and Bronner-Fraser 2008; Betancur et al. 2010; Simões-Costa and Bronner 2015; Fig. 4a, see supplementary tables S2 and S4, Supplementary Material online for NC-GRN tier classification details). These included genes involved in NC induction, specification, migration and differentiation. This last category included numerous genes associated with the development of different pigment cell lineages, such as melanophores, xanthophores, and iridophores (Howard et al. 2021) as well as factors contributing to formation of the embryonic cranial skeleton. Interestingly, expression of several ligands and receptors of four signaling pathways—Bmp, Fgf, Hedgehog, and Wnt—also diverged between species, suggesting potential unexplored differences in diverse developmental processes controlled by these pathways. Almost 40% of the identified NC-DEGs (Fig. 4b) were associated with NC cell migration, and most were differentially expressed at 4ss and 15ss (Fig. 4b). These genes (highlighted in red in Fig. 4c) exhibited variation in expression levels over time, often displaying large differences in relative expression between species at individual stages, such as 4, 15, and 18ss (Fig. 4c). In contrast, almost no NC-related genes were DE at 10ss (Fig. 4b), suggesting a degree of developmental constraint during the molecular transition from specification (4ss) to active migration (15ss) resulting in less pronounced variation between species at this stage.

The considerable number of DEGs involved in NC migration and differentiation might result from broad knock-on effects of divergence at early stages (i.e. during NC specification) within the NC program. We identified three candidate genes—sox10, prdm1a, and dicer1—known to perform multiple functions in NC development, including in specification and migration of NC cells and in differentiation of pigment cells and craniofacial cartilages. All three genes were differentially expressed at the 4ss stage, coinciding with NC specification. sox10 is a key regulator of NC specification, maintenance, migration and differentiation into multiple cell lineages, primarily neuronal and pigment cells, across vertebrates (Dutton et al. 2001; Carney et al. 2006; Kelsh 2006). prdm1a controls NC cell formation by activating foxd3 (an early NC specifier gene) and regulating sox10 in zebrafish (Hernandez-Lagunas et al. 2005; Olesnicky et al. 2010; Powell et al. 2013). dicer1 is required for craniofacial and pigment cell development together with miRNAs. It is involved in the regulation of sox10 during melanophore differentiation in zebrafish (Weiner et al. 2019). Beyond their key functions within NC-GRN, we focused on these genes to specifically investigate whether differences in expression of dicer1 and prdm1a (sox10 regulators) could potentially influence large fold change value of sox10 at 4ss and differences observed at later stages (Fig. 4c).

Both prdm1a and dicer1 were differentially expressed at 4ss, a stage concurrent with NC specification (supplementary fig. S2c, Supplementary Material online), although dicer1 was not detected in whole-mount specimens. In Astatotilapia, prdm1a was highly expressed in the prechordal plate (anterior-most tip of the neural tube) and at lower levels along the neural tube, whereas in Rhamphochromis, it was expressed at lower and uniform levels in the prechordal plate on both sides of the anterior neural tube (Fig. 4d). Considering the positive regulatory role of prdm1a on sox10 (Powell et al. 2013), interspecific expression differences in early NC ontogeny could thus influence later behavior of migratory NC cells and their differentiation in lineages regulated by sox10, such as pigment cells and cartilage. To further test this hypothesis, we next investigated in more detail the expression of sox10 in cichlid embryos to examine the behavior of migratory sox10-labeled NC cells (Dutton et al. 2001; Drerup et al. 2009).

sox10 Paralogs Originated in Teleost-Specific WGD while sox10a was Lost in Zebrafish and Cavefish

Similar to other members of soxE gene family, sox10 gene is present in two copies (sox10a and sox10b) in the genomes of cichlids and the majority of other teleosts (Lang et al. 2006; Nagao et al. 2018), with the notable exception of zebrafish—a widely utilized teleost model in biomedical research, which possesses only a single copy, sox10b (Braasch et al. 2007; Voldoire et al. 2017). Using recent genomic data of basal teleosts, eels and tarpons (Parey et al. 2023), we confirmed that the sox10 paralogs originated during the teleost-specific whole-genome duplication event and were subsequently lost in the zebrafish (order: Cypriniformes) and cavefish (order: Characiformes), suggesting a possible single loss event during evolution of ostariophysians (Fig. 3a). Such expansion of the genetic toolkit could provide opportunity for functional divergence, potentially contributing to the process of teleost diversification.

To examine the development of migratory sox10-labeled NC cells, we thus focused on both sox10a and sox10b paralogs. Our transcriptomic profiling revealed that expression of sox10 duplicates followed similar trajectories over time in both species, with a significant difference in transcript levels between paralogs observed at 10ss (two-way ANOVA and Tukey HSD, P < 0.001) (Fig. 3b). Further, the sox10 paralogs were differentially expressed between species across multiple stages of NC development, with both genes showing consistent upregulation in Rhamphochromis. Significant fold expression differences in sox10b levels were observed at earlier stages compared with sox10a and decreased for both genes with developmental time (Fig. 3c). The differences in expression levels of sox10 paralogs between and within species during embryonic NC development (Fig. 3) suggest that the NC developmental program, and more specifically NC migration, may be divergent between these two closely related species.

Differences in Expression Patterns between sox10 Paralogs Suggest Divergence in Cranial NC Development between Species

Next, we set out to characterize the precise expression patterns of sox10a and sox10b in cichlid embryos during the course of NC development. Consistent with findings in medaka (Nagao et al. 2018), the expression of both sox10 paralogs was generally observed in NC cells at all stages examined, encompassing the processes of NC specification and migration. Furthermore, this expression was detected across different axial levels, corresponding to distinct NC subpopulations (Rocha et al. 2020), as well as in the otic vesicle (Fig. 5). These findings indicate that both cichlid sox10 paralogs are likely to perform functions in NC development.

Fig. 5.

Gene- and species-specific variation in embryonic expression of sox10 paralogs during NC development in Malawi cichlids. In situ HCR images of foxd3 and sox10 paralogs in representative somite stage-matched embryos of AC and RC. a-aⁱⁱⁱ) early sox10 paralog expression partially overlaps with foxd3, supporting their expression in bona fide NCCs, whereas unique somitic domains suggest expression in non-NC cells. b) Schematic representation of sox10 paralog expression patterns at NC specification. c-cⁱⁱ) During mid-somitogenesis, sox10a is also expressed in extra-embryonic tissues in both species, here shown in RC. d) Schematic representation of four migratory streams of cranial NC in cichlids (based on expression of sox10) in lateral view. e to h) sox10 paralog expression in the cichlid embryonic head shows differences between genes and between species, implicating differences in the migratory patterns of NCCs. Schematic representations show expression overlays, highlighting overlapping and unique sox10a/sox10b-expressing cell populations. All images present maximum intensity projections of dorsal or lateral (vertical panels in eⁱ and gⁱ) views of dissected embryos, with anterior toward the top of the figure. Embryos and sections are outlined in white. Color-coded vertical bars in a) represent the AP extent of expression of each presented gene. Embryos presented in d and e representative of n ≥ 3. AC, Astatotilapia calliptera “Mbaka”; D, dorsal; NC, neural crest; nt, neural tube; RC, Rhamphochromis sp. “chilingali”; ss, somite stage; V, ventral. (Scale bar = 50 μm, 10 μm in cⁱ).

We identified variation in expression between genes within species during cranial NC specification (4ss), which manifested in temporal and spatial aspects common to both examined species. First, sox10a was expressed earlier than sox10b and concomitant with foxd3 (4ss in Fig. 5a, supplementary fig. S5a and b, Supplementary Material online). Second, sox10a was detected in cells residing bilaterally in the dorsal region of somites and co-expressed neither sox10b nor foxd3 (Fig. 5a-aⁱⁱⁱ). Notably, expression of sox10a in the somitic region was also observed at later stages in both species and localized into the extraembryonic membranes, extending laterally from the embryo proper enveloping the yolk (Fig. 5c). We propose that these sox10a+ cells may contribute to the solitary pigmented melanophores that populate yolk during somitogenesis, as they first appear on both sides of the somites in the anterior trunk region at mid-segmentation stages, prior to extensive migration. However, the embryonic origin and function of this novel population remain to be elucidated. Major differences were also observed between species, including multiple sox10a domains distributed along the anterior-posterior embryo axis at 4ss in Astatotilapia but not Rhamphochromis (Fig. 5a). These multiple distinct expression domains of sox10a, including in the somitic region, and distinct from those of its paralog sox10b, have not been reported for any of the soxE family genes in other teleost species to date (Nagao et al. 2018; Takamiya et al. 2020; Tsunogai et al. 2021). Our results thus suggest that sox10a could have been co-opted to function in both NC and non-NC cells during early somitogenesis in cichlids, whereas sox10b expression during NC specification resembles that of other vertebrates (Southard-Smith et al. 1998; Cheng et al. 2000; Aoki et al. 2003).

As the development progressed, the variation in expression levels and spatial arrangement of sox10 paralog domains visibly decreased, leading to largely overlapping patterns during cranial and trunk NC migration in both species (Fig. 5e to h, supplementary fig. S5c, Supplementary Material online). We primarily focused on cranial migration, as differences in expression patterns between paralogs and between species were most pronounced in that aspect. During cranial migration, NC cells migrate along highly conserved pathways across all vertebrates, forming four main streams: the anterior-most fronto-nasal stream, followed by the mandibular, hyoid, and postotic branchial stream (Fig. 5d; Steventon et al. 2014). Among these, subsets of cells continued to express only one of the paralogs, for example in Rhamphochromis at 30ss, the cells migrating ventrally in the mandibular stream (mdb on Fig. 5g and h) expressed only sox10a (empty blue arrowhead, Fig. 5h), whereas in Astatotilapia, cells migrating in the same stream expressed both sox10a and sox10b (blue arrowhead, Fig. 5h). Furthermore, although both paralogs were expressed in the otic vesicles (but not necessarily co-expressed by the same cells, see supplementary fig. S5d, Supplementary Material online), only sox10b was expressed in a group of cells on the dorsal surface of the forebrain, corresponding to oligodendrocytes derived from neural stem cells (Schebesta and Serluca 2009) (gray arrowhead, Fig. 5f and h). Given the conservation of this pattern across vertebrates (Suzuki et al. 2017), these findings suggest that, compared with sox10a, cichlid sox10b likely performs a broader range of the conserved vertebrate sox10 roles in fate regulation of some non-NC derived lineages, such as those derived from neural stem cells.

Several differences were also observed between species in the spatial arrangement and migratory behavior of the cranial NC cell subpopulations labeled by one or both of the paralogs until the end of somitogenesis in Astatotilapia (30ss) (Fig. 5). The most pronounced divergence involved the extents of their migration into and around head structures (e.g. otic vesicles, eyes, Fig. 5). For instance, at 18ss, sox10b+ cells migrating in the hyoid stream (hy on Fig. 5f and h), were present further ventral-laterally in Astatotilapia compared with Rhamphochromis (green arrows on Fig. 5eⁱ). This pattern persisted at 30ss (Fig. 5gⁱ). In contrast, at the same stage, cells co-expressing sox10a/sox10b and migrating in the fronto-nasal stream (fn on Fig. 5h) were found further ventrally in Rhamphochromis (purple arrowheads on Fig. 5gi and h).

In summary, cichlid sox10a and sox10b were expressed by overlapping and divergent subsets of cranial NC cells migrating along the stereotypical pathways. Notably, we observed fine-scale variation in migratory patterns of all streams, except for branchial. Such spatial differences could potentially lead to divergence in fine-scale patterning of the structures derived from cranial NC cells, i.e. craniofacial cartilages and bones, connective tissues and pigment cells (Schilling and Kimmel 1994; Wada et al. 2005; Kague et al. 2012). Combined with differences in expression levels from whole-embryo RNAseq (Fig. 3a and b), these differences could be also related to size variation of the populations expressing each paralog, especially considering the differences in embryo sizes between these species (Marconi et al. 2023). Quantitative in situ hybridization studies, such as with d- or qHCR (Choi et al. 2020), could be used in the future to compare the sizes of migratory NC populations between species. Taken together, the variation in expression patterns of sox10 paralogs in both overlapping and distinct domains indicates divergence in NC development between cichlid species and could reflect potential divergence in developmental functions between two genes.

Genome Editing Reveals Functional Divergence between sox10 Paralogs, including a Novel Craniofacial Skeletal Function of sox10a in Cichlid Fishes

Given the extent of divergence between sox10a and sox10b expression, we set out to characterize their function in cichlids. sox10 paralog function remains uncharacterized beyond the zebrafish (Danio rerio) and medaka (Oryzias latipes) model systems. In zebrafish, sox10b function is limited to pigmentation and neural derivatives (Kelsh 2006), while in medaka both paralogs have redundant functions in pigmentation development (Nagao et al. 2018). To test whether sox10 genes perform divergent roles in cichlids, we deployed CRISPR/Cas9 system in Astatotilapia to induce indel mutations in the coding sequence (exon 1) of sox10a and sox10b in turn (Fig. 6, supplementary figs. S6 to S7, Supplementary Material online).

Fig. 6.

Paralog-specific KO phenotypes indicate functional divergence of sox10a and sox10b in cichlids. a) Melanophore pigmentation defects, particularly decreased abundance, are observed in sox10b but not in sox10a A. calliptera mutants. b) sox10b-CRISPR embryos have significantly reduced melanophore counts at 7 dpf compared with WT and control injections (boxplots showing melanophore count within the dorsal head area outlined in a.; P-values are shown for Tukey HSD). c) Representative craniofacial phenotypes of sox10a-CRISPR embryos in live animals and their corresponding cartilage preparations in lateral and ventral views. White dashed lines show the reduced frontal slope of the brain case. Arrowheads and asterisks in panels showing cartilage preparations indicate reduced and missing cartilages, respectively, in the mutants. Fractions indicate the frequency of presented phenotype across all surviving embryos on that day. Data shown are representative of at least biological replicates per target gene. Ba, branchial arches; bhy, basihyal; dpf, days postfertilization; op, orbital process; pf, pectoral fin; soc, super-orbital cartilages; WT, wild-type. Scale bars = 1 mm.

From day 6 to 7 postinjection (st. 18 to 19), coinciding with the main stages of craniofacial cartilage development and patterning in this species (Marconi et al. 2023), craniofacial malformations were observed in sox10a CRISPR mosaic embryos. Neurocranial and craniofacial deformities were prevalent among injected embryos (Fig. 6a, n = 21/76 across four clutches, supplementary table S9, Supplementary Material online) and, while ranging in severity between clutch-mates, these mutants consistently exhibited flattening of the frontal bones (brain case; indicated by white dashed lines in Fig. 6a), small and bulging, forward-facing eyes as well as protruding, unmoving jaws (white arrowheads in Fig. 6c). Alcian Blue stains for cartilage further revealed severely malformed or entirely missing super-orbital cartilages, basihyal, branchial arches, and pectoral fins (gray arrowheads and asterisks in Fig. 6c). Besides craniofacial abnormalities, sox10a mutants at this stage also displayed cardiac and circulatory system defects, reduced black melanophore pigmentation and malformed caudal fin cartilages (supplementary fig. S6, Supplementary Material online). The defects in pigmentation were not quantified due to the wide range of severity of cranial deformations, which could have had indirect effects on the pigmentation, for instance due to reduced epidermis surface area for populating chromatophores. The flattened frontal skull slope suggests that the embryonic brains were also likely adversely affected. Mosaic embryos with severe sox10a-knockout (KO) phenotypes (Fig. 6c and supplementary fig. S6, Supplementary Material online) did not survive past 9 dpf (st. 22), suggesting embryonic lethality of a complete KO.

Unlike sox10a CRISPR mosaic embryos, sox10b mutants did not show any craniofacial cartilage defects at day 7 postinjection (st. 19) nor later, but instead had significantly reduced melanophore pigmentation on the dorsal head region (Fig. 6b and c, n = 13/77 across three clutches, supplementary table S9, Supplementary Material online), the first body area consistently populated by all three differentiated pigment cell types (Marconi et al. 2023). The development of other pigment cell lineages (i.e. reflective iridophores and yellow xanthophores) appeared unaffected by the induced mutations at this point, albeit with a noticeable difference in the extent of iridophore pigmentation on the flanks and yolk cover at 12 dpf (st. 24) (supplementary fig. S7, Supplementary Material online). Despite observed sox10b expression in otic vesicles and oligodendrocytes (Fig. 5e to g) and a known role of zebrafish sox10b in glial development (Carney et al. 2006); we did not identify any discernible phenotypes in KO fish involving these tissues.

These functional analyses provide compelling evidence for divergent roles of cichlid sox10a and sox10b in the development of the NC and its derivatives, cartilage and pigment cells (melanophores), respectively. Although we observed a level of potential functional redundancy between paralogs in the differentiation of pigment lineages, which requires further investigation, we also uncovered a novel and pivotal role for sox10a in the formation of cranial skeleton (neurocranium and craniofacial cartilages)—a function so far only described in cichlids (Fig. 7). Taken together, the differences in sox10a and sox10b expression and the sox10a role in craniofacial skeletal development show that the cranial NC program is diverging between cichlid species, but whether the sox10 paralogs are causally associated with this divergence remains to be tested.

Fig. 7.

Functional repertoires of sox9 and sox10 paralogs across the bony fish phylogeny. Dashed line box indicates possible gene loss among members of the superorder Ostariophysi, including zebrafish and cavefish. References for functional analyses are presented in supplementary table S10, Supplementary Material online. TSGD, teleost-specific genome duplication; WGD, whole genome duplication. Silhouettes downloaded from http://phylopic.org.

Discussion

The genetic program orchestrating development of the NC is remarkably conserved across vertebrates, despite NC-derived structures constituting some of the most diverse phenotypic traits, especially among lineages that have undergone adaptive radiation. Our study suggests that neofunctionalization following gene duplication, together with extensive transcriptomic divergence during early NC development, may have collectively contributed to the morphological diversification of NC-derived traits, including pigmentation and craniofacial shapes, in Eastern African cichlids. On a larger evolutionary scale, we report a rare example of taxon-specific, divergent evolutionary trajectories of paralogous genes originating from a single genome duplication event in vertebrates, with sox10 paralogs showing different fates in different teleost taxa.

Interspecific Divergence in Transcriptional Landscape and Signature of Positive Selection in NC-Related Genes

Recent studies have highlighted that transcriptional evolution among closely related species and across tissues often underlies phenotypic diversity (Cardoso-Moreira et al. 2019; El Taher et al. 2021). Our findings support and expand upon this concept by uncovering large-scale, distinct transcriptomic dynamics during somitogenesis between two closely related yet eco-morphologically distinct Malawi cichlid species, affecting both coding and noncoding genes. Many protein-coding genes, including those involved in NC development and its derivatives, exhibited significant variation in expression trajectories between the two species. Notably, differences in expression of many genes could be explained by simple temporal shifts in timing of gene expression relative to ss, which in turn can be attributed to differences in developmental timing (i.e. heterochrony) between species during somitogenesis (Marconi et al. 2023). These results provide evidence that variation in developmental timing, and consequently altered gene expression dynamics, contribute to species divergence in this clade.

Using combination of genomic and transcriptomic data, we identified 74 DEGs with signatures of divergent positive selection between species, with the most extreme outliers related to cellular functions, such as fgf8a (secreted signaling molecule), tspan37 (regulator of cellular signaling), and kcnk18 (potassium channel protein), and involved in development of several organs and systems, such as brain, eye, and nervous system (Rahm et al. 2014; Gebuijs et al. 2019). The diversity of their functions during development makes it challenging to identify a specific phenotype under selection as multiple traits could be involved simultaneously. Further work could examine the functions of these genes in cichlid embryonic development and verify their expression in adult tissues to provide some insights into these results.

Furthermore, our study identified dozens of DEGs with known NC functions, ranging from specification to differentiation into pigment and cartilage cell lineages. Several of these genes have previously been associated with NC-trait variation in cichlid fishes (Albertson et al. 2014) and appear to evolve under divergent positive selection, implying that they might be playing a role in the adaptive evolution of cichlids. Although some of these DEGs perform multiple functions during embryogenesis aside from NC development, the overrepresentation of DEGs involved in NC cell development and migration, including many cell-intrinsic (e.g. transcription factors) and -extrinsic factors (e.g. signaling molecules), further underscores the causal role that variation in NC migration might play in the emergence of novel NC phenotypes and species differences in cichlids, as posited by studies across vertebrates (Tucker and Lumsden 2004; Fish et al. 2014; Powder et al. 2014).

Differences in NC Migration between Cichlid Species Potentially Associated with Trait Divergence

Our comparative analysis in two eco-morphological divergent Malawi cichlids revealed that sox10 paralogs—duplicates of a master regulator within the NC program—show species-specific temporal and spatial variation, especially in migratory cranial NC cells expressing sox10 paralogs. This variation could be in turn linked to pigmentation, craniofacial diversity and potentially differences in other NC-derived cell lineages, such as cranial sensory ganglia, Schwann cells, and cardiomyocytes (Schilling and Kimmel 1994; Sande-Melón et al. 2019).

Specifically, sox10a had consistently higher expression in Rhamphochromis in somite-stage matched embryos. sox10a+ cells in the fronto-nasal stream showed more advanced in ventral migration and earlier expression in the hyoid stream, suggesting that an earlier and, potentially, prolonged migration could lead to larger structures, consistent with the prominent jaws of the piscivore compared with moderate phenotype of an omnivore. However, it is important to note that our sox10a mutants displayed relatively normal lower jaw cartilages, while the cartilages of branchial arches (originating from NC cells migrating in the branchial stream) and the frontal slope (fronto-nasal stream) were severely affected. These findings align with patterns of sox10a expression in both of these streams throughout cranial NC migration. Further studies involving later stages (e.g. pharyngula until st. 16 when the first cartilages are present, Marconi et al. 2023) and the use of additional molecular markers (e.g. sox9 genes, dlx2a, collagen gene col2a1a) are needed to investigate the potential connection between sox10a-expressing cells in the embryo, cartilage formation and these affected structures.

Our findings posit an intriguing question regarding whether sox10 paralogs are actively driving the observed differences in NC migration, or are they merely markers of migratory populations that have been influenced by activity of other genes? Given the complexity of developmental and genetic programs governing NC, it is likely that both scenarios could be at play: sox10 paralogs may contribute to driving NC migration, particularly through their regulatory (and, in case of sox10a, at least to some extent novel) roles, while also marking NC populations that are being influenced by other genes involved in migration. Further investigation into the specific roles of sox10 paralogs, including their interactions with other genes during NC migration, would be necessary to fully understand their contribution to the observed differences.

Finally, future comparative studies on NC development across multiple divergent species, particularly those characterized by different eco-morphotypes, could provide valuable insights into the connection between NC development, migratory behaviors, and adult phenotypes, and how observed differences are associated with their adult adaptive diversity. Nonetheless, our findings highlight the variability of the NC development in closely related cichlids, suggesting unexplored variation in the NC-GRN and highlighting sox10 paralogs as potential candidates involved in trait diversification.

Evolution of the NC Program following the Teleost Genome Duplication

One consequence of WGD during teleost evolution is the duplication of entire gene networks. Retained genes post-WGD can contribute to functional innovation and evolution through coding sequence changes and rewiring of the regulatory networks controlling gene expression. The latter is particularly relevant for duplicated regulatory genes, which have a higher potential for significant impacts on gene expression and phenotypic effects (Taylor and Raes 2004). Our phylogenetic analyses confirmed that sox10 was duplicated during teleost WGD and retained as duplicates in most teleost genomes, except for two members of a major freshwater fish clade Ostariophysi, zebrafish (Cypriniformes) and cavefish (Characiformes). These results suggest that the loss of one of the sox10 duplicates likely occurred in the common ancestor of ostariophysians (i.e. single loss event) but future work involving representative species from other ostariophysian orders, such as catfish, electric eel or South American knifefish, will be required to address this hypothesis. Interestingly, our expression and functional analyses data revealed distinct fates for sox10 paralogs in cichlids compared with other teleosts. This presents a unique opportunity to study the impact of gene duplication on NC evolution in multiple evolutionary contexts, including comparisons between nonteleost fishes (i.e. pre-WGD), multiple teleost fishes and other vertebrates.

For example, the differences in expression patterns between sox10a and sox10b in cichlids are more pronounced than those observed in medaka, particularly with respect to the apparent absence of sox10a expression in the extraembryonic tissues in the latter species (Nagao et al. 2018; Tsunogai et al. 2021). The acquisition of novel and divergent expression domains in cichlids suggests the presence of new regulatory elements, warranting future studies to investigate the expression and regulation of sox10 paralogs in this clade and teleosts on a broader scale.

Moreover, we showed that sox10a KO mutant is associated with aberrant craniofacial skeletal development phenotypes in A. calliptera, consistent with the emergence of functional divergence between sox10a and sox10b in cichlid fishes and other teleosts. The modest proportion of severely affected CRISPR embryos may suggest selection for reduced mutant NC cell contribution to the mosaics, where cells with functional sox10 alleles might outcompete or compensate for the affected cells, leading to a less pronounced phenotype. Additionally, compensatory mechanisms could be at play, where other genes or paralogs mitigate the loss or reduction of sox10 function, thereby minimizing the visible impact on the phenotype in many embryos. Future work examining the impact of sox10 paralog KO mutations on NC development, including specification, migration and differentiation into pigment cells and cartilage, will provide important insights into the mechanistic basis of these phenotypes.

The observed cartilage malformations in mosaic sox10a-CRISPR cichlids contrast with the effects of complete KO mutations of both its orthologs sox10a and sox10b in medaka and sox10b in zebrafish (Kelsh and Eisen 2000; Nagao et al. 2018). Functional analyses in medaka (Nagao et al. 2018) and zebrafish (Kelsh and Eisen 2000) suggest that in these lineages sox10 paralog(s) have retained ancestral functions in pigmentation, neuron and glial cell development, consistent with observations in other vertebrate lineages (Fig. 7, supplementary table S10, Supplementary Material online). Notably, medaka sox10a and sox10b have partially redundant functions in development of pigment cells, in agreement with the most common scenario of subfunctionalization of gene duplicates, also reported for sox9 paralog in zebrafish (Yan et al. 2005; Nagao et al. 2018). In these teleosts, pigmentation was severely reduced, but cartilage development remained unaffected, similar to other vertebrates (Kapur 1999; Honoré et al. 2003). The chondrogenesis aberration seen in sox10a mutants is more reminiscent of the effects of sox9a homozygous mutation in Nile tilapia (Li et al. 2023) and zebrafish (Yan et al. 2005), as well as sox9 in mice (Wagner et al. 1994; Fig. 7, supplementary table S10, Supplementary Material online).

The expression and functional data combined thus suggest a different partitioning of functions among soxE family genes (specifically sox9 and sox10) in cichlids compared with other teleosts. We posit that sox10a acquired an essential role in chondrogenesis in the cichlid lineage, a function performed by sox9 genes in other teleosts (Cresko et al. 2003; Yan et al. 2005; Nagao et al. 2018), possibly overlapping with the role of sox9a in cichlids (Li et al. 2023; Fig. 7). In contrast, cichlid sox10b appears to have retained its function in pigmentation development, akin to its orthologs in other teleost lineages. The subtle, yet significant, pigmentation phenotypes of sox10b-CRISPR fish could be due to the mosaic nature of the induced mutations. It remains to be determined whether potential differences in head pigmentation in severe sox10a mutants reflect some functional overlap between paralogs or indirect effects of the cranial cartilage malformation on pigment cell development. Similarly, the striking eye and brain size reduction in sox10a mutants suggests that chondrogenesis of the craniofacial skeleton might also influence ocular and brain development. Together with observed expression differences of prdm1a (sox10 regulator) at the anterior neural plate border and prechordal plate, these findings may also implicate early variation in cranial placode (developmental precursors of sensory systems) development between species. More comprehensive analyses of NC development in sox10 mutants, along with examination of non-NC-derived lineages expressing sox10b, will be required to better understand the functions of these genes during cichlid embryogenesis. Future studies in a wider range of species will be necessary to assess the extent of functional overlap among other soxE family members and to delineate the roles of sox10a and sox9 paralogs in cichlid chondrogenesis, cranial placode and ocular development.

Conclusion

Our study highlights the critical role of gene expression differences during embryonic development in generating phenotypic variation, particularly in the NC, a key cell population responsible for many adaptive cichlid traits. The identification of transcriptome variation linked to NC processes, alongside signatures of positive selection, underscores the importance of this cell population in the developmental evolution of cichlid fishes.

Functional divergence of duplicated genes is widely recognized for its role in the evolution of morphological diversity, including the expansion of the pigmentation pathway in teleosts (Braasch et al. 2009). In contrast, neofunctionalization is a rare occurrence in genetic evolution, with subfunctionalization (partitioning of ancestral functions between paralogs) being the most common scenario, allowing for developmental fine-tuning (Force et al. 1999). Our findings in cichlids, however, reveal that sox10a has acquired a novel function in chondrogenesis, which has not been previously reported in any other teleost clade. This shows that these paralogs have followed divergent evolutionary fates throughout teleost evolution.

We further hypothesize that, while cichlid sox10b retained its ancestral function, the neofunctionalization of sox10a may have contributed to the rewiring of the NC genetic and developmental programs. These programs show remarkable differences between cichlid species, thus providing an opportunity and genetic raw material for the remarkable diversification of the NC-derived phenotypes in cichlids. Altogether, sox10 paralogs and their variable evolutionary fates are an ideal system for studying the evolution of NC-GRNs across multiple evolutionary scales, within cichlid radiations, and among teleosts and vertebrate species.

Supplementary Material

Supplementary material is available at Molecular Biology and Evolution online.

Acknowledgments

A.M. was supported by Wellcome Trust (WT) PhD Program in Developmental Mechanisms (215223/Z/19/Z), EMS is supported by a Natural Environment Research Council (NERC) Independent Research Fellowship (NE/R01504X/1) and G.V. is supported by an EMBO fellowship (ALTF 339-2022). The authors thank the Department of Zoology Imaging Facility for assistance with microscopy, and Dr. H. Svardal (University of Antwerp), Dr. G. Turner (Bangor University), Dr. D. Joyce and A. Smith (University of Hull) for the kind gifts of fish stocks and samples. We thank Dr. B. Steventon (University of Cambridge) and members of the Santos lab for critical comments on the manuscript. For open access, the authors have applied a CC BY public copyright license to any author-accepted manuscript version arising from this submission.

Author Contributions

A.M., G.V., and E.M.S. devised the study and analyses. M.J.N. collected wild specimens to establish laboratory stocks. R.D. contributed to whole-genome data collection. A.M. collected and performed RNA extractions on the samples, conducted CRISPR/Cas9 mutagenesis, embryo imaging, phenotyping, cartilage staining and DNA collection. A.M. and G.V. performed data analysis on RNA sequencing data. A.M. and A.K. collected embryonic material, performed HCR staining and imaging. A.M., A.M., G.V., and E.M.S. wrote and revised the manuscript. All authors read and approved the final manuscript.

Data Availability

The raw RNA sequencing reads have been deposited in the Gene Expression Omnibus under the accession GSE250387. This data is made available on an open access basis for research use only. Any person who wishes to use this data for any form of commercial purpose must first enter a commercial licensing and benefit sharing arrangement with the Government of Malawi. All scripts used to analyze the data can be accessed at: https://github.com/Santos-cichlids/neural-crest-sox10-paralogs-cichlids.

Author notes

Aleksandra Marconi and Grégoire Vernaz contributed equally to this work.

Conflict of Interest The authors declare no competing interests.