Evolutionary Conservation of Five Pax6 Homologs under Strong Purifying Selection in the Horseshoe Crab Species Cluster

1. Introduction

Together with Myriapoda and Pancrustacea, chelicerates constitute one of the three monophyletic subphyla of the arthropods. Best known for fear-instilling terrestrial taxa like harvestmen, scorpions, and spiders [1], chelicerate biodiversity also includes two ancient marine lineages, i.e. sea spiders and horseshoe crabs (Xiphosura). Among these two, horseshoe crabs stand out by the conservation of canonical arthropod compound eyes [2], as part of an overall body plan design that is famous for qualifying as a “living fossil” [3]. Equally important, the Atlantic horseshoe crab Limulus polyphemus has played a pivotal role in early vision research [2].

Two vision-related gene families received exceptional attention in L. polyphemus. This includes opsins, the GPCR family of light-sensitive transmembrane receptors [4], and the developmental transcription factor subfamily Pax6 [5]. The first arthropod homologs of the Pax gene family were discovered in the fruit fly Drosophila melanogaster, which possesses two homologs of Pax6 named eyeless (ey) and twin of eyeless (toy) [6]. Comparative analyses of the past 15 years revealed that the presence of two Pax6 genes in Drosophila and other arthropods, including species from Pancrustacea, Myriapoda, and Chelicerata, is due to a gene duplication that must have taken place over 500 million years ago in the arthropod stem lineage [7]. Given this deep conservation of singleton homologs of ey and toy in arthropods, it came as a surprise to find not only one Pax6 ortholog in L. polyphemus that was originally uncovered by gene cloning efforts [5] but three orthologs of ey and two orthologs of toy in the whole genome assembly of L. polyphemus [8].

While unexpected given the generally limited range of Pax transcription factor gene family size variation in arthropods, these findings align with the comparative genomic reconstruction of three whole genome duplication (WGD) events in the xiphosuran stem lineage. Initial evidence of this exceptional trajectory emerged from the analysis of the first genome assemblies for the Atlantic horseshoe crab L. polyphemus [4,9]. In addition to L. polyphemus, extant Xiphosura includes the Mangrove horseshoe crab Carcinoscorpius rotundicauda, the Chinese horseshoe crab Tachypleus tridentatus, and the Indo-Pacific horseshoe crab Tachypleus gigas [10,11]. Comparative genomic analyses timed the earliest split of L. polyphemus from the other three species to approximately 135 million years ago and the split of C. rotundicauda from the last common ancestor of T. tridentatus and T. gigas to about 50 million years ago [12,13].

The comparative analysis of the organization of the homeobox (Hox) gene clusters in L. polyphemus, C. rotundicauda, and T. tridentatus produced evidence of at least one WGD in the xiphosuran lineage to the last common ancestor of extant horseshoe crabs [4,9]. Subsequent analyses of high-quality whole genome assemblies for all four xiphosuran species, however, produced evidence of a total of three WGDs in the xiphosuran lineage [14,15,16]. The genomes of the members of the horseshoe crab species cluster thus stand out by having been molded by one more WGD compared to the two rounds of WGD in early vertebrates [17,18,19,20] and two more rounds of WGD compared to the independent WGD discovered for terrestrial Chelicerates [21], i.e. the Arachnopulmonata (Figure 1).

While these insights shed light into the cause for the large number of Pax6 homologs in horseshoe crab species, it is unknown which evolutionary forces led to their final numbers and whether these paralogous loci serve novel or ancestral functions. Novel functions can be acquired by a process that is referred to as neofunctionalization following gene duplication [22,23]. Alternatively, the duplication of highly pleiotropic genes frequently results in the differential inheritance of ancestral functions in the descendant gene duplicates [22,23]. These basic scenarios can be distinguished by comparative analysis of substitution rates both at the amino acid and nucleotide sequence levels. Neofunctionalized paralogs, for instance, are characterized by transiently increased amino acid substitution rates as a result of positive selection on function optimizing amino acid residue replacements [24]. This acceleration is effectively measured by substitution rate comparisons between of sister gene duplicates. As one paralog continues to confer ancestral function and hence evolves under purifying selection, the neofunctionalized sister paralog will be characterized by a significantly higher relative substitution rate. In the case of subfunctionalization, by contrast, sister paralogs are less likely to differ dramatically in their protein sequence substitution rates. However, the release of constraints resulting from the reduction of pleiotropy can lead to moderate but detectable transient increases in protein substitution rates [25]. Last but not least, paralogous genes can also maintain genetically redundant functions, thereby buffering the impacts of environmental or gene regulatory variation [26,27]. One example of this is the deeply conserved redundant specification of the ocular segment region by ey and toy in arthropods [7].

An alternative approach to probing for innovative versus conservative sequence evolution trajectories in the comparison of non-synonymous versus synonymous sequence changes at the nucleotide sequence level [28]. Grounded in the neutral theory of molecular evolution [29], strong positive selection on proteins is predicted to lead to a higher rate in the fixation of non-synonymous nucleotide substitutions (dN) compared to the base rate of silent substitution accumulation (dS). By the same token, strong purifying selection on protein sequences is indicated by a lower dN compared to dS.

In this study, we pursued two aims. First, we asked whether all five Pax6 paralogs discovered in L. polyphemus were conserved in other horseshoe crab species. Our second aim was to probe for sequence evidence of innovative versus conservative gene duplication outcomes by analyzing dN/dS ratios along the horseshoe crab species cluster tree. Our efforts revealed the conservation of Pax6 gene subfamily size in the horseshoe crabs and that all five xiphosuran Pax6 homologs have evolved under strong purifying selection identifying subfunctionalization as their likeliest post-WGD trajectories.

2. Materials and Methods

Genome assemblies

Homolog searches were conducted in the GenBank assembly GCF_000517525.1 of L. polyphemus, GenBank assembly GCA_011833715.1 for the IMCB_SINMHF_001 isolate of C. rotundicauda, GenBank assembly GCA_004102145.1 for T. tridentatus isolate BBG1, and GenBank assembly GCA_014155125.1 for T. gigas isolate IMCB_SINCHF_001.

Homolog searches

BLASTn and BLASTp searches were conducted in the NCBI genome assembly user interfaces.

Gene model annotation

Genome models were hand-annotated in downloaded chromosome and scaffold text files.

Gene tree reconstruction

We generated a multiple sequence alignment of the conceptual protein sequences with T-Coffee at default settings [30]. Given the high degree of overall sequence conservation, the alignment used as direct input for gene tree estimation with Randomized Accelerated Maximum Likelihood (RAxML) [31] as made available in the CIPRES Science Gateway online user interface [32]. We applied the Jones-Taylor-Thornton (JTT) model of amino acid sequence evolution and implemented gamma-distributed substitution rates across sites with 3 categories.

dN/dS analysis

dN/dS ratios were computed using the online Ka/Ks Calculation tool (https://services.cbu.uib.no/tools/kaks) [33,34]. We submitted complete open reading frame nucleotide sequences and the species tree (Figure 1) and applied the discrete Grantham submatrix option without codon bias.

3. Results

Conservation of three homologs of ey in horseshoe crabs

Previous BLAST searches identified three distinct homologs of ey in L. polyphemus, which contrasted with the conservation of a singleton homolog of ey throughout all other arthropods investigated so far, including other chelicerates [8]. To probe for the conservation of all three L. polyphemus ey homologs in other horseshoe crab species, we searched the genome assemblies of C. rotundicauda, T. tridentatus, and T. gigas by tBLASTn with the L. polyphemus ey1, ey2, and ey3 homologs as queries. Candidate homologs were re-BLASTed against reference protein databases of L. polyphemus and D. melanogaster to filter false positive orthology assignments. These efforts and preliminary gene tree analyses identified 1:1 orthologs of L. polyphemus ey1, ey2, and ey3 in all of the other three horseshoe crab species (Figure 2 and Table S1). Of note, all of the putative ey homologs were characterized by the previously described lysine residue at position 64 of the Pax domain, which is diagnostic for ey throughout arthropods (Text document S1) [7,8].

Subsequently, we hand annotated all of the relevant loci to produce complete conceptual transcript, open reading frame, and protein sequences for ey1, ey2, and ey3 of C. rotundicauda, T. tridentatus, and T. gigas (Text document S2 - 4). The multiple sequence alignment of the protein sequences revealed a high level of sequence conservation among the xiphosuran ey1, ey2, and ey3 orthologs (Text document S1). Most sequence differences stemmed from single amino acid residue replacements. The ey3 homolog of L. polyphemus stood out by a two amino acid long deletion and a three amino acid residue long alanine string insertion. Both of these larger scale differences mapped to sequence variable regions of ey based on outgroup comparisons with the singleton ey orthologs of the common house spider Parasteatoda tepidariorum and the red flour beetle Tribolium castaneum (Text document S1).

The global gene tree analysis confirmed the presumed orthology relationships of the individual ey homologs but did not resolve the relationship between the ey1, ey2, and ey3 ortholog clusters (Figure 2).

Conservation of two homologs of toy in horseshoe crabs

Previous studies also identified two homologs of toy in L. polyphemus [8]. To probe for their conservation in the horseshoe crab species assembly, we proceeded in the same way as in the investigation of the conservation of the three homologs of ey. As a result of this effort, we found that also both L. polyphemus homologs of toy were conserved in C. rotundicauda, T. tridentatus, and T. gigas (Table S1; Text document S1-3). We further found that all toy homologs were consistently characterized by the previously described diagnostic arginine residue at position 64 of the Pax domain [7,8]. The global gene tree analysis confirmed the presumed orthology relationships of all individual toy homologs (Figure 2).

As in the case of the ey homologs, multiple sequence alignment of the xiphosuran toy protein sequences revealed a high level of sequence conservation (Text document S1). Most sequence differences stemmed from single amino acid residue replacements. The toy1 homolog of T. gigas was characterized by the deletion of a proline residue in a sequence variable region C-terminal of the homeodomain (Text document S1). Moreover, the toy1 homolog of L. polyphemus was characterized by the insertion of three amino acid residues in a variable sequence region even further C-terminal of the before-mentioned T. gigas variable site (Text document S1).

Strong purifying selection on all horseshoe crab Pax6 homologs

The conservation of complete open reading frames for all xiphosuran Pax6 transcription factor homologs provided strong evidence of their continued functionality in the horseshoe species cluster. To gain further insights into the functionalization trajectories of these critical regulatory genes following three WGDs and the three speciation events (Figure 1), we took advantage of the overall high conservation at the nucleotide sequence level to probe for selection pressure by dN/dS ratio analyses (Figure 3).

Consistent with strong functional conservation, all five paralogs were characterized by dN/dS values below zero. Averaged across the entire phylogeny of the xiphosuran species tree, the dN/dS ratios of toy1 and toy2 were very similar with values ranging between 0.08 and 0.12, respectively (Table S2). Compared to these points of reference, the dN/dS values of 0.02 for ey1 and 0.04 for ey3 suggested a slightly higher strength of purifying selection compared to toy1 and toy2. The ey2 homolog, finally, stood out with the highest across phylogeny dN/dS value of 0.22, indicating a possibly weaker degree of purifying selection.

Overall, these findings corroborated the functionality of all five Pax6 transcription factor paralogs in the Xiphosura and were best compatible with subfunctionalization or conserved functional redundancy following gene family expansion as a result of WGD.

4. Discussion

As part of the bilaterian gene toolkit [35], most members of the Pax transcription factor gene family perform essential regulatory gene network functions during embryonic and postembryonic tissue and cell specification processes [36]. Successful execution of these pivotal functions depends on precisely regulated expression in time and tissue context. A systematic analysis in Drosophila demonstrated the critical effects of these variables as the misexpression of Pax genes leads to dramatic morphological phenotypes in the developing head [37]. In the specific case of the arthropod Pax transcription factor paralogs ey and toy, these aspects are highlighted by the extreme nature of lack of function and gene misexpression consequences. Both, ey and toy mutant Drosophila strains are characterized by depletion of prominent components of the peripheral visual system, i.e. the compound eyes and ocelli, respectively [38]. Opposite outcomes are observed as the consequence of misexpressing ey and toy, which famously causes the formation of extra-compound eye retinal tissues in antennae and leg appendages [39].

A third critical variable underlying proper Pax transcription factor function is defined by the level of expression, i.e. dosage. This variable can be conceived to be particularly relevant in the context of gene duplication. This is because the complete duplication of the cis-regulatory content of a gene locus is predicted to result in the duplication of overall expression levels in the absence of dosage-compensating mechanisms or preceding allelic differentiation. Considering these implications, it is of little surprise to find that most chelicerates possess only a single homolog of each ey and toy despite the duplication of gene content through one round of WGD in the early Arachnopulmonata (Figure 1) [8,40,41,42]. This outcome aligns with the assumption that dosage level increases of ey and toy are likely to result in deleterious fitness consequences, which were selected against by positive selection on post-WGD individuals with deletion mutations removing one sister paralog of each ey and toy.

Given the evidence of not only one but three WGD events in the early xiphosuran lineage (Figure 1) [16], the existence of higher numbers of ey and toy homologs in extant horseshoe crabs fits expectations. Indeed, the fact that none of the ey and toy paralogs are tandem linked but are, instead, spread out over different chromosomes or scaffolds is consistent with their origins in the wake of WGDs. Applying the terminology of ohnologs for paralogous genes that originated via WGD [43], it is safe to conclude that toy1 and toy2 represent ohnologs with the respect to each other and that the same is the case for ey1, ey2, and ey3.

And yet, given the significance of dosage for proper Pax gene function, the existence of three ohnologs of ey and two ohnologs of toy in the Xiphosura still constitutes a challenging conundrum. Their endurance in the genomes of extant horseshoe crabs raises two questions. The first is how the likely fitness-reducing effects of gene duplication on expression levels were resolved. The second question is which mechanisms explain the subsequent long-term conservation of the persisting paralogs. Of course, all the while posing these questions, the five xiphosuran Pax6 genes also represent indirect evidence of the assumed fitness-reducing effects resulting from Pax homolog number increases. This is because the three WGD events in the early xiphosuran lineage must have spawned a total of eight ohnologs for each ey and toy, the majority of which were selected against to extinction.

From the copious research on ey and toy in insects and arachnids [6,40,41,42,44,45,46], we know that both ancient Pax6 paralogs are critical factors in a large number of cell- and tissue-specification processes, i.e. that they are highly pleiotropic [47]. This background leads to the straightforward prediction that the xiphosuran ey and toy homologs diversified by inheriting non-overlapping components of the highly pleiotropic functional landscape of their ancestral singleton precursor homologs. This scenario is very likely the answer to the puzzling early finding of the presumed lack of Pax6 homolog expression in the developing visual organs of L. polyphemus, which was based on the assumption of the exclusive presence of a singleton Pax6 homolog [5]. Subsequent annotation efforts identified the Pax6 homolog molecularly cloned by Blackburn et al. (2008) as the representative of the xiphosuran toy1 homolog [8] (Table S1). At this point, it is reasonable to hypothesize that one of the ey1, ey2, or ey3 paralogs is expressed in the early ocular segment and later in the peripheral visual organ precursor cell populations, while the other paralogs may be differentially involved in brain or insulin receptor specification [48]. It is further tempting to speculate that this is more likely the case for ey1 or ey3 given the lower strength of purifying selection on ey2 based on the findings in our dN/dS analysis. Tissue expression studies, ideally conducted by embryonic whole-mount hybridization as conducted by Blackburn et al. (2008), are now warranted to probe these predictions. The second possibility that awaits testing by expression studies is that some of the five xiphosuran Pax6 transcription factors engage redundantly in specific gene regulatory networks. This scenario is predicted by the deeply conserved redundant roles of ey and toy in the specification of the ocular segment during arthropod embryonic development [7].

The comparison with the functionalization fates of ey and toy in other arthropod lineages is critically informative to rank the exceptionality of the preservation and functionalization trajectories of the five horseshoe crab Pax6 genes. Similarly instructive is the comparison with the evolutionary fates of Pax6 transcription factor homologs that resulted from the two rounds of WGD in the vertebrate stem lineage (Figure 1). In this group, three paralog descendants of the Pax6 singleton ancestral precursor locus have remained conserved in select lineages, thereby defining what can be referred to as the vertebrate Pax4/6/10 subfamily [49,50]. The existence of maximally three post-WGD descendants in the well-sampled vertebrate lineages implies the early extinction of at least one Pax4/6/10 group member. In contemporary vertebrates, Pax6 is most pervasively conserved while Pax4 and Pax10 are only preserved in a small number of taxa [49,50]. Combined, these data paint a similar picture of paralog losses and subfunctionalization in comparison to the fates of ey and toy paralogs after three rounds of WGD in the horseshoe crabs.

Our study complements previous analyses of the effects of the xiphosuran WGDs on other gene families, i.e. opsins [4], the Hox gene complex members [9,14,15], the components of the JAK-STAT signaling pathway [12], and micro RNAs [16,51]. Combined, these efforts provide ample testimony of the unexpectedly eventful genomic history of horseshoe crabs, which will continue to serve as an important paradigm case for genome evolution studies.

Looking ahead, the next obvious question is how the duplication and conservation histories of other Pax transcription factor subfamilies compare to the Pax6 subgroup in the horseshoe crabs. Our ongoing analyses indicate a wide range of preservation outcomes with maximum subfamily sizes of up to seven members (Kassem et al., in preparation). We, therefore, predict that the horseshoe crab lineage will emerge as particularly useful to study the roles of constraints and opportunities in the evolutionary diversification of Pax transcription factors besides serving as an exceptionally well-defined sample case for understanding the complex consequences of WGD [52,53].