Sp185/333-positive BAC clones were identified in the sea urchin large insert BAC library

The Sp185/333 gene family in the S. purpuratus genome sequence (ver. 3.1; June 15^th, 2011; (http://www.echinobase.org)) shows unexpectedly few members [21, 29]. This was likely due to the repetitive nature of the Sp185/333 genes and flanking sequences [20, 21] that led to gene collapse by the assembly algorithm and the formation of artificially hybrid Sp185/333 gene sequences. To address this problem, to identify as many of the members of the Sp185/333 gene family as possible, and to characterize the structure of the family, a screen of the large insert (~140 kb) arrayed BAC library (25X genome coverage [33]) was conducted as described [20] and 75 Sp185/333-positive BAC clones were identified. To identify a subset of BAC clones for sequencing and insert assembly, clones were initially evaluated by PCR amplification using four primer sets (F2/R6, F6/R9, R2/R9, F5/R1; Fig. 1b; Additional file 1: Table S1) and the resulting amplicon sizes were used to verify the existence and sizes of Sp185/333 genes, and to identify BAC clones that may have spanned the same genomic sequence (as expected based on library coverage). Of the 75 BAC clones that were evaluated, 27 supported PCR amplification and resulted in six different Sp185/333 gene amplicon patterns (Additional file 2: Figure S1A). Patterns 1–3 were found in six, eight and nine BAC clones, respectively, whereas patterns 4–6 were only in one or two BACs each. In accordance with [21], PCR analysis of the intergenic distances in the BACs indicated that the majority of the genes were about 3.0–3.5 kb apart (Additional file 2: Figure S1B). The 48 BAC clones that did not support PCR for Sp185/333 sequences may have lacked annealing sites for the primers, and consequently they were rescreened by Southern blot for Sp185/333 sequences. Of these BAC clones, only BACs 3020I13 and 4069G2 showed visible bands for all restriction digests (Not1 plus either SalI, Xho1 or BamHI) and BAC 4069C2 showed a weak band only in Sal1/Not1 digest (Additional file 3: Figure S2). Variations in the positions of restriction sites suggested differences in the sequences among these BACs. These results indicated that the Sp185/333 gene family in the sequenced genome of S. purpuratus was likely included in 30 BAC clones. The remaining 45 BACs were not evaluated further.

Three Sp185/333 gene clusters were identified in the Sp185/333-positive BAC clones

The Sp185/333 genes have a wide range of sizes [20], however differences in length among some genes may only be a few nucleotides (nt), which cannot be detected by standard gel electrophoresis of gene amplicons (see Additional file 2: Figure S1A). Furthermore, three duplicated copies of genes with the D1 element pattern (Fig. 1a) in a single BAC insert has been reported (Additional file 4: Figure S3A) [21]. Because amplifying duplicated genes that result in a single amplicon size does not provide reliable estimates of gene copy numbers in BAC inserts, we implemented fragment length analysis that is regularly used to identify microsatellite length polymorphisms in populations [34]. This approach increased the resolution of the estimated gene sizes (a proxy for element pattern; see Fig. 1a) and copy numbers in each BAC. The most variable region in the Sp185/333 genes, which is the 3′ half of the second exon (Fig. 1a) [20], was amplified from BAC DNA with F6 and R9-FAM primers (Additional file 1: Table S1; Additional file 2: Figure S1B) to generate amplicons within the size limits of the capillary sequencer. This approach was validated with BAC clone R3-3033E12, which had been previously sequenced, assembled and verified [21] and chromatograms of four amplicon sizes identified accurately each of the Sp185/333 gene element patterns (A2, B8, D1 and E2; Fig. 1a, Additional file 4: Figure S3). This approach also identified the three copies of the D1 genes based on the amplicon peak height and area, which was approximately three times that relative to the other peaks for this BAC clone. Based on the verification of this approach, fragment length analysis was used to evaluate the 27 BACs that supported PCR amplification of Sp185/333 genes. Gel images of non-fluorescent amplicons (Additional file 2: Figure S1A) and fragment length chromatograms of fluorescent amplicons for each BAC suggested six patterns of peaks with one to five amplicon sizes that were characteristic of different sets of BACs (Fig. 2). The signatures of the amplicon patterns for each BAC were characterized by both the amplicon length plus each fragment peak height and area, and the ratio of each peak relative to the tallest peak in the chromatogram for each BAC clone. Results were similar to the gel electrophoresis analysis showing three gene amplification signatures (Fig. 2) that were identified in 23 of 27 BACs and designated as Clusters 1, 2 and 3 (equivalent to patterns 1–3; Additional file 2: Figure S1A). In addition, three less common patterns were identified in only four of the BACs that partially overlapped either Cluster 1 or 2 and were designated as Clusters 1′, 1″ and 2′ (equivalent to patterns 4–6; Additional file 2: Figure S1A). Amplicon sizes reflect directly the element patterns of the genes and were used to predict different sets of genes. Results from this initial analysis of the Sp185/333-positive BACs suggested the Sp185/333 genes were arranged in at least three clusters.

BAC inserts were assembled using two high through put sequencing approaches

PacBio sequencing resolves the Sp185/333 gene cluster sequences

Because the insert assemblies based on Illumina reads did not resolve the sequences of the Sp185/333 gene clusters in 10 of 11 of the BAC insert assemblies, five BAC clones (10B1, 10M18, 3104P4, 4074J14, 10K9) were selected for re-sequencing by PacBio. BAC clones chosen for re-sequencing were based on whether the Illumina assemblies resulted in a single contig (Table 1) and if the gene clusters agreed with results from fragment length analysis (Fig. 2). Because PacBio sequence reads are significantly longer than both the Illumina reads and the sizes of the repeats associated with the Sp185/333 genes, reads were expected to include non-repetitive sequence flanking the repetitive regions that could be used as anchors for improved assemblies. Each bar-coded PacBio library resulted in 3 to 17 × 10³ reads with an average read length of 6.4–7.5 kb. All five PacBio BAC assemblies resulted in a single contig without gaps (Table 1) and four of the five BACs (with the exception of 10M18) showed complete Sp185/333 gene sequences of expected length based on the positions of Sp185/333 primer sequences in the assemblies. The PCR amplicon sizes for genes and intergenic regions matched to the corresponding sequence sizes in these four assemblies and alignments showed expected full-length gene sequences separated by intergenic regions (represented as dot plots in Fig. 3). Results were supported by comparisons of BACs 10B1 and 4074J14 insert assemblies (not shown) to the overlapping and previously sequenced BAC insert (GenBank BK007096) [21].

The BAC insert assemblies represent two loci in the genome

To determine whether all of the BAC insert assemblies represented a single locus as inferred from the single cluster of Sp185/333 genes in the genome sequence, or whether additional loci could be identified, all BAC sequences assembled from either PacBio and/or Illumina reads were aligned to each other. The regions flanking the gene clusters that did not contain repeats and that assembled more readily were used primarily to determine the positions of matches and overlaps among the sequences. Alignments of the flanking regions indicated either ~99 % sequence identity or no match (represented by dot plots, Additional file 5: Figure S4). When all BAC insert sequences were mapped relative to each other overlaps among the sequences were identified and indicated the presence of two separate loci (Fig. 4). Locus I encompassed nine BACs that either fully or partially overlapped and was composed of Clusters 1 and 2, including the Sp185/333 gene cluster previously reported (R3-3033E12, GenBank accession number BK007096; [21]. The remaining two BACs overlapped to compose locus II, which covered Cluster 3 and was not present in the sea urchin genome sequence.

The Sp185/333 gene repertoire in the sea urchin genome sequence

A detailed analysis of the BAC insert assemblies identified three clusters of Sp185/333 genes (Fig. 5), which agreed with results from the fragment length analysis (Fig. 2). The clusters included 15 predicted genes (seven in Cluster 1, six in Cluster 2 and two in Cluster 3; Fig. 5) that ranged in size from 1170 to 1894 nt. The identification of two loci for the Sp185/333 gene clusters based on matching flanking sequence for Clusters 1 and 2 (Fig. 4) strongly suggested that they may reside at the same locus in the genome and therefore may to be allelic. However, we could not rule out the possibility that Clusters 1 and 2 might be positioned within very large and recent duplications. Cluster 1 spanned 43 kb and contained seven genes with five different element patterns; A2, B8, three copies of D1 (designated as D1y, D1g and D1b) and E2 (see Fig. 1a for element patterns), of which all had been reported previously (indicated by gray shading in Fig. 5) [21]. An additional gene of the short 01 element pattern was located on an overlapping BAC and was ~400 nt beyond the right end of the previously reported BAC insert (see right – left orientation in Fig. 5). The five internal genes in Cluster 1 were tightly linked (3–3.4 kb apart), whereas the outer genes, A2 and 01, were separated from their nearest neighbor gene by ~12.5 and ~7.3 kb, respectively (Fig. 5). The orientation of the internal genes (D1y, D1g, D1b and B8) was in the opposite direction with respect to the flanking genes (A2, E2 and 01; Fig. 5). The allelic Cluster 2, which spanned 39 kb, had a similar structure, but with six genes of element patterns A2 and B8 (designated A2a and B8a), two copies of D1 (designated D1d and D1e) and two copies of E2 (designated E2a and E2b). The internal genes in Cluster 2 (B8a, D1d, D1e and E2a) were positioned 3.1–3.4 kb apart and the flanking genes, A2a and E2b, were separated from the nearest neighbor genes by 12.8 and 7.2 kb respectively (Fig. 5). The differences in gene composition between Clusters 1 and 2 included the number of D1 and E2 genes plus the presence of 01 in Cluster 1. Similarities included the A2 and B8 genes in addition to the order, orientation and spacing among the genes. It is noteworthy that if these two clusters represent allelic regions, as suggested from matching flanking regions, they do not have a matching number of genes, which confuses the concept of alleles. It is not clear which of the duplicated D1 and E2 genes should be considered alleles and may only be determined by more detailed alignments, percent identities and gene phylogenies (Barela Hudgell and Smith, unpublished). Consequently, because of the difficulties in knowing which of the genes constitute alleles, they will therefore all be called genes.

Cluster 3 was quite different from Clusters 1 and 2; it was much shorter, spanning 6.4 kb, and contained only two genes of element patterns C4 and D1f (Fig. 1a) that were in the same orientation and separated by 3.5 kb (Fig. 5). An allelic region for Cluster 3 was not identified in this study, which may be based on the possible bias of the BAC library, or that the allelic region may be so similar that allelic BAC sequences may appear the same. Although the three gene clusters were different with regard to sizes and numbers of genes, they showed strikingly similar organization with regard to gene orientation, their order relative to element patterns, and their intergenic spacing.

STRs flank the Sp185/333 genes

Previous results showed the existence of tandem and interspersed repeats within the second exon of the Sp185/333 genes [20] plus STRs of GA and GAT sequences on both sides of the genes [21]. To determine whether this unique structure was consistent for the entire gene family, the sequences of the three Sp185/333 gene clusters on BAC clones 10B1, 10 K9 and 3104P4 were evaluated for tandem sequence repeats of increasing size from 2 to 342 nt per repeat, which were identified in all clusters. STRs flanked every Sp185/333 gene, and the GA STRs with 12–132 repeats were located 303–790 nt from the start codon and 312–792 nt from the stop codon of each gene (Fig. 5, Additional file 6: Table S2). When the GA STRs were evaluated on a single DNA strand in the BAC sequences, reading from left to right (see orientation in Fig. 5) without regard to the Sp185/333 gene orientation, GA STRs (triangles above; Fig. 5) flanked the peripheral genes and complementary CT STRs (triangles below) flanked the internal genes for both Clusters 1 and 2. The GA STRs flanked both genes in Cluster 3, and long stretches of GA STRs with 39–1048 repeats that were positioned to the left of C4 (~3.1 and ~12.8 kb) and to the right of D1f (~3.4 kb; Fig. 5). These regions correlated with the relative locations of the A2/A2a, B8/B8a and E2a/D1b genes in Clusters 1 and 2 (Fig. 5, Additional file 6: Table S2). The large GA STRs also included some TC repeats plus up to 10 % interspersed repeats of TA, GG or AA (Additional file 6: Table S2). The large GA STRs located to the left of C4 also included regions of TC repeats and multiple smaller regions of non-STR sequences (Additional file 6: Table S2). Although tandem repeats are common in the sea urchin genome (http://www.echinobase.org/Echinobase/repeats) the STR locations on both sides of each of the 15 predicted Sp185/333 genes in the three clusters and surrounding segmental duplications were unique, had not been identified in similar orientations to genes in other families, and reflected an unusual but repeatedly identified pattern in this gene family.

Tandem segmental duplications bounded by GAT STRs that mostly included D1 gene copies were another feature of the three Sp185/333 gene clusters and were of similar size (4.2–4.6 kb; Fig. 6a). A partial segmental duplication of ~1.9 kb was present in Cluster 3 that did not include a gene (indicated as 3* in Cluster 3; Fig. 6a). All segmental duplications were positioned in the same orientation in each cluster and were located up to 100 nt apart. Phylogenetic analysis of the full length GAT bounded segmental duplications resulted in a clade composed of segmental duplications from Clusters 1 and 2 in a pectinate structure (see [35]) (Fig. 6b). The segmental duplication from Cluster 3 was most closely associated with one from Cluster 2 and both were positioned at the base of the tree in a sister clade. In general, the sequence similarities among the segmental duplications among the three gene clusters suggested that they may have originated from a single sequence.

In addition to segmental duplications that included the D1 genes, GAT STRs (4–68 repeats) also surrounded the E2 and 01 genes in Cluster 1 and the E2a gene in Cluster 2 (Fig. 5; Additional file 7: Table S3). The duplications that included E2 and E2a were about the same length (~7.1 kb) and were highly similar (a pairwise distance of 0.005 based on an alignment using ClustalW in MEGA). A single flanking GAT STR (14 repeats) was located to the left of E2b in Cluster 2 (Figs. 5 and 6; Additional file 7: Table S3), however no flanking GAT STR was identified to the right of E2b based on searches of the overlapping BAC 3020I13 (Fig. 5). Consequently, the positioning of the GAT STRs surrounding E2b did not conform to the pattern identified for the D1, E2 or E2a genes. The level of identity among the segmental duplications that included the D1 genes and for the duplications harboring the E2 and E2a genes suggested that the segmental duplications from all clusters may be recently generated.

The Sp185/333 gene family locus in the genome sequence is a hybrid of Clusters 1 and 2

The repetitive nature of the Sp185/333 gene family created a significant genome assembly problem that caused the Sp185/333 genes to be assembled artificially into five likely hybrid genes (genome assembly ver. 3.1, Scaffold 125; 02/18/2016) as predicted by [21]. Furthermore, allelic genes were purposefully omitted during genome assembly because they would have appeared as tandem gene duplications. Consequently, to determine which cluster was most similar to the Sp185/333 gene cluster in the genome sequence, the five genes (A1, B3, D1, [D1 or D8] and E2) in scaffold 125 (62.5 kb; 469500–531999 nt) plus the intergenic sequences and flanking regions were compared to the three clusters of genes reported here using YASS [36] (results not shown). The cluster in scaffold 125 had two D1 genes, a single E2 gene and did not include an 01 gene, and was therefore similar to both Clusters 1 and 2. The A1 gene sequence in scaffold 125 was a hybrid of both A2 in Cluster 1 and A2a in Cluster 2. Similarly, the E2 gene sequence in scaffold 125 was a hybrid of E2 in Cluster 1 and E2a in Cluster 2. The GA and GAT STRs were present in the expected positions in the gene cluster in scaffold 125, and the intergenic spacing was similar but not identical to the spacing in the BAC clusters. There were no similarities between Cluster 3 and the gene cluster in scaffold 125, and there was no evidence of Cluster 3 in the sequenced genome. In general, this comparison suggests that the Sp185/333 gene cluster in scaffold 125 of the sequenced genome is likely a hybrid of both Clusters 1 and 2, which supports the possibility that these two clusters may be allelic.

Sp185/333 gene repertoire varies among individual sea urchins

Sp185/333 family size prediction

The number of Sp185/333 genes predicted for a given genome using three different estimates (reviewed in [29]) is 50 ± 10 genes, which is much greater than that reported here for the BAC assemblies and for scaffold 125 in the assembled sea urchin genome sequence. We screened the sea urchin large insert BAC library comprehensively for Sp185/333-positive clones employing several methods for BAC analysis. This effort resulted in three Sp185/333 gene clusters that were repeatedly identified in a total of nine BACs that covered either all or some genes for each cluster. However, we cannot rule out the possibility that additional genes may have escaped identification, which is based on two additional genes isolated from the small insert BAC library that was constructed from the same genomic DNA as that used for the genome assembly [20]. It is feasible that both libraries had their own biases with regard to numbers of BACs with inserts harboring members of the Sp185/333 gene family.

The low copy number of genes identified in the sea urchin genome sequence may also be an outcome of the assembly of the genome into scaffolds. The 42 chromosomes for S. purpuratus fall into three size ranges and are relatively similar within a given size range, but have not been characterized extensively as to the location of genes [39]. This type of gene mapping, which is typically done by in situ hybridization, is unavailable as a template onto which to map relative gene locations and linkage groups to assemble scaffolds into chromosomes. Furthermore, the assembly of a genome that is limited to the level of scaffolds can be the outcome of failures to assemble scaffolds across specific genomic regions that contain repeats [2, 3]. This may be attributed, in part, by additional Sp185/333 gene clusters such as Cluster 3 that is not present in the genome sequence. Assembling a genome for an outbred, diploid organism is challenging and the outcome is likely to be a hybrid sequence of both haplotypes, as is the case for the hybrid gene sequences in scaffold 125 that are similar to both Clusters 1 and 2. Although the number of Sp185/333 genes in the sequenced genome is lower than expected, it is possible that the individual sea urchin that was chosen for genome sequencing had a particularly small Sp185/333 gene family. Variations in gene family size among individual sea urchins have been suggested previously [20, 29] as has been demonstrated here.

STRs may drive duplications and deletions

The genes in the three Sp185/333 clusters are organized in the same general order, orientation, and spacing, which leads to the notion that these clusters may have originated from a single ancestral gene or cluster that was significantly expanded by processes that include duplications and deletions of genes, segments and clusters into tandem or ectopic locations [7]. The proposed genomic instabilities may result, in part, from the presence of both GA and GAT STRs in the clusters. Historically, STRs have been considered “junk” DNA, however evidence suggests the functional importance of some non-coding DNA including promoting genomic instability and increased rates of changes to the DNA sequence [3–5]. The STRs themselves are prone to an estimated 10–100,000 fold higher mutation rates (that may include changes to nearby genes) compared to other regions of the genome that do not harbor STRs [4]. In the Sp185/333 gene family, certain genes appear to have been duplicated more recently, including copies of both D1 and E2 genes based on both being positioned within almost identical segmental duplications. As has been suggested previously [21], the GAT STRs may be involved in the diversification process based on their unique positioning at the edges of the D1 segmental duplications and some of the E2 segmental duplications. Although the mechanisms by which this may occur is not understood, one possibility is that the duplications are the outcome of DNA breakage followed by repair that is based on misalignments to similar but non-identical sequences on the homologous chromosome [7]. During this process, chromosomal misalignment or meiotic mispairing, may lead to the deletion and duplication of genes by unequal crossing over between homologous chromosomes resulting in changes to the gene family size and composition in the subsequent generation. Segmental duplications of regions that include the D1 genes in the BAC insert sequences are consistent with profiles of Sp185/333 genes in all 10 sea urchins in which amplicon sizes that correspond to D1 genes are present in a relatively higher copy number. On the other hand, the GA STRs that surround each of the Sp185/333 genes may have important functions in genomic gene organization, including gene duplication and deletion [21]. Deletions are inferred from the long stretches of GA STRs that flank both sides of Cluster 3 and match precisely to locations of genes in Clusters 1 and 2. Furthermore, the combination of both segmental duplication and gene deletion may have been the basis for the partial segmental duplication in Cluster 3 (identified as 3*; Fig. 6a). This region matches the 5′ regions of the GAT STR segmental duplications that include D1 genes but is missing the 3′ region of the segmental duplication that would be expected to include a gene. The missing 3′ region of this segmental duplication may have undergone a GA STR mediated deletion. However the underlying mechanisms for these processes are not understood.

The GA and GAT STRs are present in the clusters in both orientations, in accordance with the orientations of the genes, and constitute inverted repeats with complementary sequences. The formation of STR based secondary structures such as stem loops within the Sp185/333 gene clusters may result in replication fork stalling during S phase of the cell cycle [5, 40], which can induce DNA double strand breaks [4]. This damage would be followed by repair with sequence from the homologous chromosome, resulting in sequence changes, mutations and gene recombination [4] among the Sp185/333 gene clusters. Inverted complementary repeats may also be a basis, in part, for the genomic instability in the Sp185/333 gene family, and examples of expected outcomes are hybrid Sp185/333 genes resulting from recombination, which has been predicted from computational analysis of sequenced genes [31], and the mismatched alleles in clusters 1 and 2.

There are no pseudogenes in the Sp185/333 gene family

In addition to the STRs, the Sp185/333 gene family has at least three other types of repeats: direct and interspersed repeats within the second exon, elements with shared sequences among genes, and segmental duplications. Similar to STRs, each one of these types of repeats may cause chromosomal misalignment during meiosis that may lead to variations in gene copy number plus sequence changes to individual genes resulting from homologous DNA repair processes. However, because these types of repeats appear within genes and include coding sequences (in contrast to the STRs), they may also lead to gene conversion, as has been suggested previously [7, 21]. The general lack of pseudogenes in the Sp185/333 gene family (only a single pseudogene has been identified, although it was one of two intronless genes that may be retroposons [20]), is very unexpected for such tightly linked genes that share sequences and are surrounded by STRs (see [19]). One speculation is that gene conversion may correct pseudogenes in this family and that the STRs may block continued strand exchange and sequence homogenization of the entire cluster [21]. Because partial gene conversion may result in what might appear as gene recombination, potentially leading to frameshifts, missense sequence and early stop codons, the lack of pseudogenes in the Sp185/333 gene family suggests that regulatory mechanisms may be involved.

Immune gene diversification

All organisms must have means to stay ahead of and protect themselves against their pathogens and to control their commensals. This is particularly important given that microbes have multiple mechanisms for diversification [41], which are then selected for virulence based on improved ability to overcome immune responses and to promote colonization and infection. Typically, host protection in vertebrates is accomplished by diversifying the immune system, which includes the somatic recombination of the Ig gene family mediated by the recombination activating genes (RAGs) in jawed vertebrates [42] or the assembly of variable lymphocyte receptor genes in jawless vertebrates by activation induced cytidine deaminases (AID) [43, 44]. However, immune gene diversification also occurs in the absence of RAG and/or AID activity in many organisms and tends to be associated with clustered genes and genomic instability. This is evident from the KIR gene family, which is tightly clustered and shows variations in gene copy numbers for different human haplotypes in addition to sequence diversity in the population based on 15 to 112 alleles for different loci [9]. Variations in the KIR region is thought to result from genomic instability driven by shared sequences within and surrounding the genes, including minisatellite sequences in the first intron, all of which can lead to crossing over among the tandemly arranged genes, expansion and contraction of gene numbers from meiotic mispairing, in addition to recombination, gene conversion, domain shuffling, duplications and deletions and single nucleotide polymorphisms ([45, 46], reviewed in [47, 48]). Similar to the assembly problems for the Sp185/333 gene family, the attributes of the KIR family have also resulted in very poor assembly in macaque genome sequences [49]. The major histocompatibility complex, with which the KIRs interact, is a cluster of duplicated genes that includes pseudogenes and is highly diversified in human populations with multitudes of polymorphic alleles, recombination hotspots and sequence variations [50]. In invertebrates, the colonial hydroid, Hydractinia symbiolongicarpus has two allorecognition (alr) loci, which control fusion and rejection among individuals, are positioned in genomic regions with pseudogenes and duplication that show sequence similarities to the expressed alr genes, in addition to putative transposable elements and regions of repetitive sequences [11, 12]. The Fusion/Histocompatibility (FuHC) locus in the colonial tunicate, Botryllus schlosseri, controls allogeneic fusion and rejection reactions among these sessile animals and shows sequence diversity, gene clustering, partial gene duplications and an extraordinary number of polymorphic alleles in the population for a subset of the genes (reviewed in [51]). Finally, excellent examples of repeat driven diversity are the anti-pathogen R gene families in higher plants, of which many are present in clusters of duplicated similar genes [15, 52]. R gene diversity results from unequal crossing over, gene conversion, duplications and deletions, all of which may lead to hybrid R genes or new genes with advantageous anti-pathogen functions [53]. The common theme among these examples is clustered similar genes (that may be considered as large repeats), plus repeats within and/or near genes, which together all lead to genomic instability. The emergence of new genes with new functions that are selected based on successful pathogen interactions and are therefore advantageous to the fitness and survival of the host reinforce and perpetuate the gene systems that function on the instability of specific genomic regions.

We have proposed here and in previous communications [7, 21, 31] a range of possible mechanisms that could result in sequence diversification of the Sp185/333 genes in addition to variations in the size of the family among individuals. To date, diversification mechanisms that may result from regional instability in the genome have not been tested experimentally for the Sp185/333 gene family (or in any other gene system that diversifies based on genomic instability) and our speculations are based on computational analyses of the sequences within and surrounding the genes. The rate of change within the Sp185/333 genes can be deduced from the observations that there are no shared gene sequences among sea urchins, that there is a wide variety of mosaic element patterns, and that adjoining elements among genes may have different sequences [31]. Gene recombination can be inferred from the lack of correlation between the patterns of tandem repeats (elements 2–5, see Fig. 1a) compared to the patterns of interspersed repeats (elements 11–26, Fig. 1a) of the second exon among different genes. Detailed computational analyses of the Sp185/333 gene sequences has inferred i) swift evolution of the tandem type 1 repeats located in the 5′ half of the second exon (Fig. 1a) based on predictions of duplications, deletions and recombinations, ii) rates of gene recombination that may be comparable to VJ recombination for the T cell receptor α locus, and iii) recombination that may occur at any point along the length of a gene [31]. In addition, there is no correlation among cDNA sequences when comparing the sequence of the first exon to the sequence of the adjacent 5′ untranslated region [24]. When the diversification rates of the Sp185/333 genes are estimated based on molecular clock analysis, results suggest that the extant genes shared a common ancestor 2.7–10 million years ago, which overlaps with the estimate of S. purpuratus divergence from its sister species [29, 32]. This suggests ongoing recombination and that the extant genes are relatively young [31] in agreement with the Sp185/333 gene family showing differences among different individuals. Although these computational predictions have not been tested experimentally, the results are consistent with observationsfor the gene family structure.

Clusters 1 and 2 appear to be allelic based on matching sequences that flank the clusters and the sequence of genomic scaffold 125 that is a hybrid of both clusters. However, the non-matching loci within the two clusters suggests larger scale genomic modifications that alter gene copy number and changes to the members of the clusters. Repetitive regions in genomes are associated with genomic instability, including meiotic mispairing that results in significant changes to gene clusters that share sequence and mispair during meiosis leading to unequal crossing over that alters the numbers of genes in one homologous chromosome relative to the other (reviewed in [7]). The outcome of such large scale changes can be observed as mismatched alleles and variations in the gene family size among individuals in the population. Many of these changes have been reported for R genes in higher plants [14, 15, 53], and the characteristics of many R gene families parallel those of the Sp185/333 gene family [7]. When genomic changes resulting from meiotic mispairing, segmental duplications, gene duplications and deletions, gene conversion and single point mutations are considered together, the level and rate of sequence diversification for members of a gene family is likely to be astounding.

Clustered immune genes with shared sequences in a wide range of plants and animals show changes that can range in impact on the host from altered coding regions to the formation of hybrid genes and gene conversion. Although these changes lead to pseudogenes in most gene families, changes can also benefit the species when new gene sequences are selected as being advantageous based on non-self interactions that result in improvements to pathogen (or allogeneic) detection (e.g., see [54]). The genomic instability in the Sp185/333 gene family and the apparent allelic mismatch between Clusters 1 and 2 may be the result of duplication/deletion events in the specific individual sea urchin that was employed for BAC library construction and genome sequencing. The outcome of such swift diversification of the Sp185/333 gene family within the population can be deduced from the variety of gene amplicons derived from 10 different sea urchins in which only some fragment sizes are shared among individuals.

BAC DNA isolation and Sp185/333 gene verification

BACs (75) were identified by screening the large insert (~140 kb) sea urchin BAC library prepared from sperm DNA from a single animal [33] as reported in [21] and obtained from Eric Davidson and Andrew Cameron at the Center for Computational Regulatory Genomics, the California Institute of Technology. Escherichia coli (DH10B) transfected with BAC clones were cultured over night at 37 °C in 2xYT culture media [55] with 12.5 μg/ml chloramphenicol. BAC DNA was isolated by an alkaline lysis protocol as described [20], and the presence of the vector (pBACe3.6) was confirmed by PCR amplification using 3.6F and 3.6R primers (Additional file 1: Table S1). BAC inserts were released from pBACe3.6 vector by NotI (New England BioLabs) digestion, and evaluated for size by PFGE (1 % agarose, pulsed field certified, 6 V/cm; switch time of 1–15 s over 16 h) by comparison to the MidRange PFG Marker I (New England BioLabs). BAC DNA was subjected to multiple PCR amplification reactions using five sets of Sp185/333 degenerate primer pairs including four sets to amplify the genes (5′UTR/3′UTR, F2/R6, F6/R9, F2/R9) and one pair to amplify the intergenic regions (F5/R1; Additional file 1: Table S1; Additional file 2: Figure S1B). Primers were designed previously based on 121 unique Sp185/333 cloned gene sequences and had been tested for their ability to amplify different regions within and between the genes [20, 21]. PCR reactions of 20 μl and 27 cycles included 50–100 ng BAC DNA template, 1X GoTaq ready mix (Promega), 0.1 μM each primer using the appropriate annealing temperature for the primer pair (see Additional file 1: Table S1) with elongation time of 30 s to 4 min depending on the expected amplicon size.

Southern blots

BAC DNA was double digested with Not1 plus either SalI, Xho1 or BamHI (New England BioLabs) and fragments were separated on a 0.7 % agarose gel. DNA was transferred to nylon membranes (GeneScreen Hybridization Transfer Membrane; PerkinElmer) by standard capillary blotting [55] and membranes were blocked with hybridization solution according to [56]. Blots were evaluated with ³²P-riboprobes generated with either T3 or T7 RNA polymerase from three, linearized gene clone templates (2-034 (GenBank accession number EF607716), 4-1521 (EF607770) and 4-1543 (EF607784)) according to [20, 21]. The specificity of the ³²P-riboprobe to hybridize with Sp185/333 sequences was evaluated with three Sp185/333 genes of different element patterns in addition to an SpC3 cDNA [57] that served as the negative control (Additional file 3: Figure S2A). Probes were hybridized to the filters in hybridization solution at 42 °C with rocking overnight according to [56]. Membranes were washed in 4X, 2X and 1X SSC (20X SSC is 0.3 M Na citrate, 3 M NaCl, pH 7) with 1 % SDS at 65 °C, dried and exposed to X-OMAT AR film (Eastman Kodak) for 1 or 3 h at room temperature or 24 h at −70 °C with an intensifying screen.

Genomic DNA isolation

Sperm was obtained from 10 sea urchins that were spawned by exposure to a non-lethal electric shock (6–10 volts). DNA was extracted from ~5 μl of dry sperm using DNeasy Blood & Tissue Kit (Qiagen) according the manufacturer’s instructions for blood extraction.

Fragment length analysis

Sp185/333 genes of different sizes and copy number in the BAC clones (50 ng) or from genomic DNA (5 ng) were estimated by fragment length analysis using F6 and R9 modified with 6-fluorescein amidite (R9-FAM; Integrated DNA Technologies; Additional file 1: Table S1; Additional file 2: Figure S1B) to amplify the 3′ half of the second exon of the Sp185/333 genes. Labeled amplicons were analyzed by electrophoresis on a 1 % agarose gel to estimate quantities, and 0.5–2 μl of the product were mixed with 10 μl Hi-Di formamide containing 2.5 % GeneScan™ 1200 LIZ® size standard (ThermoFisher) and loaded onto an ABI PRISM 3130 Genetic Analyzer for capillary electrophoresis operated by ABI3130 Data Collection software (ThermoFisher). Peaks were analysed by GeneMarker® HID ver. 1.90 (SoftGenetics) based on length in nt, relative height and relative area of each peak according to the following calculation: peak intensity (height or area) divided by the tallest peak (or largest area) in a given sample, which resulted in fractional ratios of each individual peak in a sample.

BAC sequencing and insert assembly

BAC DNA (20 μg) was isolated using the Large-Construct Kit (Qiagen). Contamination with E. coli (DH10B) genomic DNA was detected by amplification of two regions of the genome using primer pairs Ec1F/Ec1R and Ec2F/Ec2R (Additional file 1: Table S1). The level of contamination was determined by comparing E. coli amplicon intensity to amplicons of similar size from pBACe3.6 vector using primer pairs; 3.6-1F/3.6-1R and 3.6-2F/3.2-2R (Additional file 1: Table S1) that amplified two different genomic regions of 1 kb. E. coli DNA contamination that was 10 fold less than the amount of isolated BAC DNA was deemed sufficient for sequencing. Illumina sequencing and assembly for 11 BAC clones were conducted at the J. Craig Venter Institute (Rockville, MD). Two libraries were constructed for each BAC clone: a 500 nt paired-end library and a 5 kb mate pair library into which barcodes were incorporated to allow multiplexed sequencing. Subclone libraries from all BACs were sequenced on one lane of an Illumina MiSeq sequencer with an average read length of 250 nt. Reads were processed with BOWTIE2 [58] to remove matches to E. coli (DH10B) and pBACe3.6 sequences. Reads were separated by bar-code index, and adapter and primer sequences and poor quality regions (base quality < 10) were trimmed using TRIMMOMATIC [59]. Because read coverage varied greatly across different regions of individual BACs, coverage was normalized to 45X using DIGINORM (http://ged.msu.edu/papers/2012-diginorm), which had been used for the Drosophila melanogaster genome assembly [60]. Reads were assembled with ALLPATHS-LG (ver. R47449 [61]) with ploidy = 1 and MIN_CONTIG = 200. Intra-scaffold gaps were closed when possible using GapFiller [62]. The ligation sites between the pBACe3.6 vector and the insert were identified by reads that crossed these regions. The reads were assembled to reconstruct the vector ends; vecL was next to the EcoRI site at position 10 in the vector, and vecR was next to the EcoR1 site at position 2801 in the vector. Reconstruction of the vector – insert ligation sites allowed verification that the entire insert sequence had been obtained and enabled the orientation of the insert sequence relative to the vector.

Single molecule real time (SMRT; Pacific Biosciences, Menlo Park, CA) sequencing and assembly was carried out at the University of Maryland Genome Research Center (Baltimore, MD). DNA was sheared to 6–20 kb fragments using a g-Tube (Covaris) and libraries were constructed using the DNA Template Prep Kit and the DNA/Polymerase Binding Kit (Pacific Biosciences). Libraries were loaded onto SMRT Cells, and sequenced with the DNA Sequencing Kit (Pacific Biosciences). Assemblies of reads for each BAC were performed using Celera Assembler (ver. 82; (http://wgs-assembler.sourceforge.net/wiki/index.php?title=Main_Page)) and HGAP3 (SMRT Analysis ver. 2.3.0) using standard parameter sets. The read filtration procedure HGAP3 was used to assemble each bar-coded set of reads independently on each set of de-multiplexed reads. The best assembly was selected based on manual review and metrics including contig count, contig N50, assembly size compared to BAC digests evaluated by PFGE (see above), and gene profiles obtained by fragment length analysis (see above). Sequences that overlapped the insert and the vector were detected and Minimus2 [63] was used to overlap and trim the ends of the contigs. The final contig set for each BAC was improved using Quiver [64]. pBACe3.6 vector sequences at the edges of the insert were identified after assembly using VecScreen NCBI tool (http://www.ncbi.nlm.nih.gov/tools/vecscreen/) and removed manually. Correctness of the sequence outside the Sp185/333 gene cluster was verified and corrected manually based on pairwise alignments for the SMRT and Illumina assemblies that covered the same region.

Verification of BAC insert assemblies and identification of Sp185/333 genes

The identification of Sp185/333 gene sequences and their positions within the BAC insert assemblies was established using GenePalette [65] based on the presence and locations of the Sp185/333 primer sequences (R1, F2, F5, F6, R6 and R9; Additional file 1: Table S1; Additional file 2: Figure S1B) within the two exons and the verification of the gene ends based on the conserved primer sequences at the 5′ and the 3′ ends of the genes (5′UTR, 3′UTR; Additional file 1: Table S1; Additional file 2: Figure S1B). Amplicon sizes (estimated by gel electrophoresis images and fragment length analysis) of intragenic and intergenic sequences were compared to sizes and distances predicted in the insert assemblies as calculated by the locations of the primer sequences. Pairwise alignments of the assemblies to self and to each other in addition to the generation of dot plots were done using Global align (https://blast.ncbi.nlm.nih.gov/Blast.cgi). The element pattern for each gene was identified from alignments to 121 known unique Sp185/333 gene sequences [20] using BioEdit [66] and genes were named according to matches to previously characterized gene element patterns (see Fig. 1a) [26].

Identification of segmental duplications and short tandem repeats

Segmental sequence duplications were identified by dot plot analyses of each assembly to itself using Global align (https://blast.ncbi.nlm.nih.gov/Blast.cgi). STRs were identified using Tandem Repeat Finder [67] with alignment parameters of match = 2, mismatch = 7, indels = 7, and a minimum alignment score = 40. The large STR regions located to the left and right of Cluster 3 were analyzed with the parameters of match = 2, mismatch = 3, indels = 5, and a minimum alignment score = 30. Searches for GATs used the parameters of match = 2, mismatch = 3, InDels = 5, minimum alignment score = 30, maximum period size = 3. STR sequences were screened by eye to identify internal repetitive sequences that did not conform to the major repeat sequence and that occupied up to 10 % of the STR sequence. Additional searches for GAT (three repeats) and GA (four repeats) STRs were done with GenePalette [65] including allowances for mismatches.

Alignments and phylogenic trees

The segmental duplications bounded by GAT STRs and that included the D1 genes were collected from the Sp185/333 gene clusters and the regions surrounding B8 and B8a of corresponding size and position as the D1 segmental duplications were selected as outgroup sequences. Alignments were conducted using MEGA (ver. 7; (http://www.megasoftware.net)) using default parameters for ClustalW and were not improved manually. Phylogenic trees were generated by both maximum likelihood and maximum parsimony using pre-set parameters with 500 bootstrap iterations.