Evolution of stickleback spines through independent cis-regulatory changes at HOXDB

Abstract

Understanding the mechanisms leading to new traits or additional features in organisms is a fundamental goal of evolutionary biology. We show that HOXDB regulatory changes have been used repeatedly in different fish genera to alter the length and number of the prominent dorsal spines used to classify stickleback species. In Gasterosteus aculeatus (typically 'three-spine sticklebacks'), a variant HOXDB allele is genetically linked to shortening an existing spine and adding an additional spine. In Apeltes quadracus (typically 'four-spine sticklebacks'), a variant HOXDB allele is associated with lengthening a spine and adding an additional spine in natural populations. The variant alleles alter the same non-coding enhancer region in the HOXDB locus but do so by diverse mechanisms, including single-nucleotide polymorphisms, deletions and transposable element insertions. The independent regulatory changes are linked to anterior expansion or contraction of HOXDB expression. We propose that associated changes in spine lengths and numbers are partial identity transformations in a repeating skeletal series that forms major defensive structures in fish. Our findings support the long-standing hypothesis that natural Hox gene variation underlies key patterning changes in wild populations and illustrate how different mutational mechanisms affecting the same region may produce opposite gene expression changes with similar phenotypic outcomes.

Fig. 1. Genetic mapping, expression and role of HOXD11B in stickleback dorsal spine development.

a, Gasterosteus mapping cross. b, QTL scan results for spine number and spine length. x axis: Gasterosteus chromosomes; y axis: LOD score for three- versus four-spine trait (top), length of DS2 (bottom). The QTL peak on chromosome 6 includes the HOXDB cluster (gene diagram at the bottom, scale bar, 1 kb). The peak on chromosome 4 includes the EDA-MSX2A-STC2A cluster described elsewhere^,. Dashed lines: genome-wide significance thresholds from permutation testing. c, Integration of GFP reporter using CRISPR–Cas9 upstream of the endogenous HOXD11B locus of low-spine Gasterosteus. Plasmid: grey; eGFP: green; basal hsp70 promoter: blue; chromosomal locus: black. Scale bar, 100 bp. TSS, transcription start site. d, eGFP expression in posterior half of fish at the stage when the dorsal spines are forming (Swarup stage 31). Scale bar, 1 mm. e, Note expression in fin fold between DS2 and DSL, DSL and dorsal fin (DF). Scale bar, 1 mm. f, X-ray of uninjected Gasterosteus (top) and Gasterosteus injected at the single-cell stage with Cas9 and sgRNA targeting the coding region of HOXD11B (bottom). Arrows: two blank pterygiophores are often located between DS2 and DSL but only in uninjected fish (insets: two blank pterygiophores in n = 5 out of 18 control and n = 0 out of 23 injected F0 mutants, two-tailed Fisher’s exact test P = 0.01). Scale bar, 5 mm. g, Length comparisons of dorsal and anal spines. Box and whisker plot: centre line, median; box limits, interquartile range (IQR); whiskers, 1.5× IQR; individual measurements shown as single points (circles: WT; triangles: mutant). y axis: residuals after accounting for standard length of fish (Extended Data Fig. 2a). DSL and AS were significantly longer in injected than uninjected fish (two-tailed t-test Bonferroni-corrected at α = 0.05, n = 18 control and n = 23 injected, DSL P_adj = 3 × 10⁻⁵, AS P_adj = 0.02). DS1 and DS2 lengths were not significantly different. Source data

Fig. 2. HOXDB genes show cis-acting expression differences in Gasterosteus spines.

a, F1 progeny were generated in a cross between a low-spine and high-spine Gasterosteus and tissues were isolated from up to seven indicated locations (DS1, DS2, DS3 (if present), blank pterygiophore (Pter), DSL, DF and AF) to measure allele-specific gene expression in the fin fold stage. Note DS3 only developed in some F1 progeny, so this location has fewer samples (n = 6 for DS3; n = 18 for all other tissues). b, The box plots show the ratios of high-spine to low-spine allele expression at each of three HOXDB genes. The y axis is the log₂ of the high-spine versus low-spine read ratio at a SNV (black line: equal expression at a log₂ ratio of 0). The x axis shows the seven tissues collected from fin fold stage fish arranged from anterior to posterior, as well as the sample collected from earlier whole embryos (embryo). Centre line, median; box limits, IQR; whiskers, 1.5× IQR; each measurement is represented by a single point. SNVs scored for each dorsal tissue compared to the anal fin: HOXD11B, chrVI:17,756,571; HOXD9B, chrVI:17,764,664; and HOXD4B, chrVI:17,783,616 (gasAcu1-4). Differences were significant in the anterior spines for all three HOXDB genes (DS1: HOXD11B (chrVI:17,756,571) P = 3 × 10⁻⁷; HOXD9B (chrVI:17,764,664) P = 6 × 10⁻⁶; HOXD4B (chrVI:17,783,616) P = 1 × 10⁻⁵; DS2: HOXD11B (chrVI:17,756,571) P = 3 × 10⁻⁷; HOXD9B (chrVI:17,764,664) P = 4 × 10⁻⁷; HOXD4B (chrVI:17,783,616) P = 4 × 10⁻⁷; DS3: HOXD11B (chrVI:17,756,571) P = 4 × 10⁻⁴; HOXD9B (chrVI:17,764,664) P = 8 × 10⁻⁴; HOXD4B (chrVI:17,783,616) P = 6 × 10⁻⁴; all P values by two-tailed Mann–Whitney U-test). **P ≤ 1 × 10⁻³, ***P ≤ 1 × 10^−6. All alleles with 0 reads have been replaced with 0.5 for graphical representation purposes and statistical analysis. NS, not significant.

Fig. 3. Spine number and DS3 length are associated with the HOXDB locus in Apeltes.

a, X-rays of Apeltes quadracus from Louisbourg Fortress with 3–6 dorsal spines. Scale bar, 5 mm. b, Association mapping of Apeltes from Louisbourg Fortress (n = 211) and Tidnish River 3 (n = 121). Both populations show a significant association between spine number and the HOXDB locus. Apeltes from Louisbourg Fortress were also phenotyped for dorsal spine length and showed a significant association between DS3 length and the HOXDB cluster. The markers used are displayed across the top (1–12) and the peak marker (6) is highlighted in red. The dotted line represents the Bonferroni-corrected significance threshold at α = 0.05. The 95% confidence intervals (CIs) (2-LOD) for spine number are denoted by the bars at the bottom (grey for Louisbourg, black for Tidnish). The overlapping 95% CI for spine number is approximately 5,750 bp from the third exon of HOXD11B to the first exon of HOXD9B. The smallest interval shared by both spine number and spine length intervals is approximately 2 kb, including HOXD11B exon 3 and part of the intergenic region between HOXD9B and HOXD11B. Additional anatomical details and association plots for other spine lengths are shown in Extended Data Fig. 7. Source data

Fig. 4. HOXDB genes show cis-acting expression differences in Apeltes spines.

a, Outline of Apeltes fin fold stage fry. Tissues were isolated from DS1, DS2, DS3, DS4, Pter, DSL, DF and AF to measure allele-specific gene expression in F1 hybrids. (DS4 only developed in some progeny.) b, Three HOXDB haplotypes segregating in cross. Black lines: association mapping markers. Pink: regions with genotypes associated with low-spine phenotypes in Fig. 3b. Yellow: regions with genotypes associated with high-spine phenotypes in Fig. 3b. Lighter shading: regions where marker association is unknown but DNA variants are shared between haplotypes of the same colour. c, Box plots showing allele-specific expression ratios in all tissues dissected from fry heterozygous for the H and L haplotypes (n = 4). All box and whisker plots: centre line, median; box limits, IQR; whiskers, 1.5× IQR; each measurement is represented by a single point. Reads from DS3, DS4, Pter, DSL and DF were compared to reads from AF to determine significance (DS3: HOXD9B (chr06:16,028,519) P = 3 × 10⁻⁷; HOXD11B (chr06:16,020,516) P = 9 × 10⁻⁴, two-tailed Fisher’s exact test). DS1 and DS2 were not assessed because read counts were too low (Extended Data Fig. 5). d, Box plots showing allele-specific expression ratios in all the tissues dissected from fry heterozygous for LHR and H haplotypes (n = 3). Allele-specific expression was seen in DS3 compared to anal fin for both HOXD11B and HOXD9B (DS3: HOXD9B (chr06:16,027,923) P = 9 × 10⁻⁸; HOXD11B (chr06:16,020,516) P = 1 × 10⁻³, two-tailed Fisher’s exact test). e, Box plot showing allele-specific expression ratios in all tissues dissected from fry heterozygous for L and LHR haplotypes (n = 4). Only HOXD9B is shown because HOXD11B lacks informative variants between L and LHR haplotypes (DS3: HOXD9B (chr06:16,027,923) P = 0.08, two-tailed Fisher’s exact test). **P ≤ 1 × 10⁻³, ***P ≤ 1 × 10⁻⁶.

Fig. 5. The genomic region between HOXD11B and HOXD9B contains a conserved AxE showing sequence changes in both Gasterosteus and Apeltes.

a, The protein-coding exons of HOXD11B and HOXD9B are shown in Gasterosteus (gasAcu1) genomic coordinates. Sequence conservation: phastCons conserved sequence regions identified in exons and in an approximately 500 bp intergenic region from comparisons between fish genomes. The conserved non-coding region overlaps an assay for transposase-accessible chromatin using sequencing (ATAC-seq) peak from Medaka embryonic stage 19 (ref. ) and partially overlaps the genomic intervals defined by spine phenotype and RNA expression changes in Apeltes. In high-spine Gasterosteus, the conserved region AxE is deleted (as indicated by a black line between the two grey boxes) and ERV and LINE sequences are inserted (in red and yellow, respectively). b, Approximately 600 bp AxE regions from low-spine and high-spine Apeltes were cloned into a Tol2 GFP expression construct and injected into Gasterosteus embryos. Both versions drove expression in the tailbud of embryos (left) and the fin fold, spines, and dorsal, anal and caudal fins of stage 31 fry (right), confirming that the region acts as an enhancer. Scale bar, 1 mm. c, There are four sequence differences in the AxE region of high- and low-spine Apeltes alleles: one microsatellite variation; an 18 bp indel; a SNP; and 2 adjacent SNPs. Only the single and two adjacent SNPs are within the region implicated by Apeltes recombination and RNA-seq differences (pink bar).

Fig. 6. Model of dorsal spine pattern changes resulting from coding and expression differences in laboratory and wild stickleback populations.

Typical Gasterosteus aculeatus (three-spine sticklebacks) have a convex pattern of spine lengths where the middle/DS2 spine is the longest). All three HOXDB genes have no or low expression in the anterior, DS1 and the strongest expression in the posterior DSL regions. CRISPR-induced coding region mutations in HOXD11B lengthen DSL. This is expected since the loss of a posterior Hox gene typically results in the anteriorization of structures; in this case, DSL is longer and thus looks more like DS2 (left). In contrast, the naturally occurring Boulton high-spine allele has expanded HOXDB gene expression anteriorly, such that all the HOXDB genes are expressed in all of the spines. Corresponding morphological effects seen in the QTL cross (shortening of DS2 and addition of a new spine on a blank pterygiophore between DS2 and DSL) are both consistent with posteriorization of dorsal midline structures (middle bottom). These morphological and gene expression changes are linked to cis-regulatory changes at the HOXDB locus, including the deletion of an enhancer (AxE) and insertion of two transposable elements. Typical Apeltes quadracus fish have a concave pattern of spine lengths (middle/DS3 spine the shortest). In low-spine Apeltes, HOXDB genes are strongly expressed from DS3 to the posterior. In high-spine Apeltes, the anterior boundary of HOXDB expression is shifted to the posterior and thus DS3 has lower HOXDB expression. Loss of expression of HOXDB genes is associated with lengthening of DS3 and the formation of a new spine on a pterygiophore between DS3 and DSL, both consistent with anteriorization of dorsal midline identities.

Extended Data Fig. 1. Spine phenotypic variation in Gasterosteus populations and F2 hybrids.

a. Bodega Bay male parent and b. Boulton Lake female parent of the QTL cross. Note that the anterior-most dorsal spine in Boulton is located at a position intermediate between the first two spines in typical Gasterosteus and is here labeled DS because it cannot be unambiguously assigned as DS1 or DS2. Representative F2 progeny showing a range of spine numbers, including: c. and d., two dorsal spines; e. and f., three dorsal spines; g. and h., four dorsal spines. Scale bar is 5 mm. i. The distribution of DS2 residual lengths (DS2.res) in F2 progeny of the QTL cross (n = 340). In the box and whisker plot: center line, median; box limits, interquartile range; whiskers, 1.5x interquartile range; each measurement is represented by a single point. Red and blue points are included for DS2 and DS of the Bodega Bay and Boulton parents of the cross, respectively. The residual was calculated with respect to standard length and sex. j. Distribution of DS2.res as a function of genotype at the peak marker on the distal end of chr6 (Fig. 1b).

Extended Data Fig. 2. QTL mapping of other spine lengths and axial traits.

a. Schematic of Gasterosteus anatomical features. Most Gasterosteus have three dorsal spines that in this study are referred to as dorsal spine 1 (DS1), dorsal spine 2 (DS2), and dorsal spine last (DSL). The dorsal side of the fish also has median bony plates known as pterygiophores, some of which underlie dorsal spines. Typical A-P midline pattern: two non-spine bearing/blank pterygiophores (Pter1 and Pter 2), dorsal spine 1 on pterygiophore 3 (Pter3), dorsal spine 2 on pterygiophore 4 (Pter4), non-spine bearing pterygiophore 5 (Pter5), and dorsal spine last on pterygiophore 6 (Pter6). The three unpaired fins are shown in light gray: dorsal fin (DF), caudal fin (CF), and anal fin (AF). The anal spine (AS) is also indicated on the ventral side of the fish. The standard length shown with the dotted line is from the anterior tip of the jaw to the posterior of the hypural plates. b. QTL plot of DS1 length c. QTL plot of DSL length d. QTL plot of pterygiophore number e. QTL plot of total vertebral number f. QTL plot of caudal vertebral number. Dotted lines represent genome-wide significance thresholds. Abdominal vertebral number and anal spine length were also tested, but they did not result in any peaks that passed the genome wide significance threshold. The significance threshold (dashed line) is based on LOD scores obtained in 1,000 permutations of the phenotype data (α = 0.05).

Extended Data Fig. 3. Embryonic expression of Gasterosteus HOXDB genes.

In situ hybridization of Gasterosteus aculeatus embryos at Swarup stage 19/20 a. HOXD4B probe. The arrowheads point to the hindbrain, the asterisk indicates the neural tube, and the arrows point to the anterior somites; b. HOXD9B probe; c. HOXD11B probe. d. eGFP expression at embryonic Swarup stage 19–20 in the tailbud somites from the reporter gene integrated upstream of HOXD11B in low-spine Gasterosteus (Fig. 1). The eGFP pattern recapitulates the in situ hybridization results for HOXD11B in panel C. Dotted circle indicates the location of the eye. Scale bar is 1 mm.

Extended Data Fig. 4. Coding mutations in HOXD11B cause length changes in stickleback spines in a second anadromous Gasterosteus population.

a. X-rays of an uninjected sibling control RABS Gasterosteus (top) and a RABS Gasterosteus that was injected at the single cell stage with Cas9 and an sgRNA targeting the coding region of HOXD11B (bottom). Scale bar is 5 mm. b. Quantification of spine length changes. In the box and whisker plot: center line, median; box limits, interquartile range; whiskers, 1.5x interquartile range; each measurement is represented by a single point (circle for wild type and triangle for mutant). DS1 and DS2 were not significantly different between controls and HOXD11B mutant fish. DSL and AS were significantly longer in the F0 mutants compared to the controls (two-tailed t-test Bonferroni-corrected at α = 0.05, DSL padj = 4E-13, AS padj = 2E-07, n = 38 injected and n = 30 control from 3 clutches combined). The y-axis is the residual after accounting for the standard length of fish. Source data

Extended Data Fig. 5. Hox gene expression patterns in Gasterosteus and Apeltes spines and fins.

The expression patterns for each Hox gene in different stickleback Hox clusters are shown with normalized read count on the y-axis and tissue site on the x-axis. The tissues are organized by position from anterior to posterior along the dorsal side of the fish, with the anal fin at the end. The read count shown is the average across all samples for that species; the reads are normalized within each species but not between species. The genes with empty plots exist in both species but are not expressed in the tissues shown with the exception of HOXB6B and HOXB7A, which are located in gaps in the Gasterosteus assembly and thus were not scored; HOXB6B is not expressed in Apeltes, but HOXB7A is expressed. HOXA1A, HOXB1B, and HOXB1B are present in the genomes but were not expressed and are not shown. The genes differentially expressed (padj < 0.01) between the largest anterior spines (DS1 and DS2) in Gasterosteus low-spine (three-spine) fish are HOXA2A, HOXA5A, HOXA10B, HOXC3A, HOXC5A, HOXC9A, HOXC10A, HOXC11A, HOXD3A, HOXD4A, HOXD9A, HOXD10A, HOXD4B, HOXD9B, and HOXD11B; the genes differentially expressed (padj < 0.01) between DS1 and DS3 in Apeltes low-spine (four-spine) fish are HOXA10B, HOXC4A, HOXC8A, HOXC9A, HOXC10A, HOXD9A, HOXD10A, HOXD11A, HOXD4B, HOXD9B, and HOXD11B.

Extended Data Fig. 6. HOXDB Gasterosteus nucleotide variants used for allele-specific expression analysis.

With the exception of HOXD11B, variants were located in the 3’UTR of the genes. Their locations are indicated by red lines and they are numbered from 5’ to 3’. Note that two different crosses were done to generate the embryo and spine/fin samples and because the populations are not inbred the informative variants were not the same. The table below each gene diagram shows which SNVs were informative in which set of samples.

Extended Data Fig. 7. Spine anatomy and trait association mapping in Louisbourg Apeltes.

a. Schematic of anatomical structures in an Apeltes fish with five dorsal spines. Typical A-P midline pattern: two non-spine bearing pterygiophores (Pter 1 and 2), dorsal spine 1 (DS1) on pterygiophore 3 (Pter3), dorsal spine 2 (DS2) on pterygiophore 4 (Pter4), dorsal spine 3 (DS3) on pterygiophore 5 (Pter5), dorsal spine 4 (DS4) on pterygiophore 6 (Pter6), three non-spine bearing pterygiophores (Pter7–9), and dorsal spine last (DSL) on pterygiophore 10 (Pter10). The three unpaired fins are shown in light gray: dorsal fin (DF), caudal fin (CF), and anal fin (AF). The anal spine (AS) is indicated on the ventral side of the fish. The standard length shown with the dotted line is from the anterior tip of the jaw to the posterior of the hypural plates. b. The association between HOXDB genotypes and length of DS1, DS2, DS4, DSL, and AS were not statistically significant. For significant results with DS3 length, see Fig. 3.

Extended Data Fig. 8. Spine number variation in Tidnish Apeltes.

X-ray of a. two-spine, b. three-spine, c. four-spine, d. five-spine, and e. six-spine fish. Scale bar is 5 mm.

Extended Data Fig. 9. Targeting of AxE in an anadromous Gasterosteus population also causes length changes in DSL and AS.

a. Representative uninjected LITC sibling control fish (top) and injected AxE F0 mutant (bottom). Scale bar is 5 mm. b. Quantification of spine length difference. In the box and whisker plot: center line, median; box limits, interquartile range; whiskers, 1.5x interquartile range; each measurement is represented by a single point (circle for wild type and triangle for mutant). The residual after adjusting for standard length is on the y-axis and the spines ordered from anterior to posterior are on the x-axis. DS1 and DS2 do not show a significant difference in length between control and injected. DSL and AS were significantly longer in the injected compared to the control (two-tailed t-test Bonferroni-corrected at α = 0.05, DSL padj = 3E-04, AS padj = 0.02 n = 32 control and n = 24 injected). Source data

Extended Data Fig. 10. dN/dS values for HOXD11B between pairs of stickleback species.

The tree on the left shows phylogenetic relationships of extant stickleback species (branch lengths not drawn to scale,^,). The rate of non-synonymous to synonymous substitutions in HOXD11B is higher than 1 for Gasterosteus and Apeltes comparisons (yellow shading).