s2fieldsofstudy:
Biology
Computer Science

text:
DHCR24 associates strongly with the endoplasmic reticulum beyond predicted membrane domains: implications for the activities of this multi-functional enzyme

Cholesterol synthesis occurs in the ER (endoplasmic reticulum), where most of the cholesterogenic machinery resides. As membrane-bound proteins, their topology is difficult to determine, and thus their structures are largely unknown. To help resolve this, we focused on the final enzyme in cholesterol synthesis, DHCR24 (3β-hydroxysterol Δ24-reductase). Prediction programmes and previous studies have shown conflicting results regarding which regions of DHCR24 are associated with the membrane, although there was general agreement that this was limited to only the N-terminal portion. Here, we present biochemical evidence that in fact the majority of the enzyme is associated with the ER membrane. This has important consequences for the many functions attributed to DHCR24. In particular, those that suggest DHCR24 alters its localization within the cell should be reassessed in light of this new information. Moreover, we propose that the expanding database of post-translational modifications will be a valuable resource for mapping the topology of membrane-associated proteins, such as DHCR24, that is, flagging cytosolic residues accessible to modifying enzymes such as kinases and ubiquitin ligases.

INTRODUCTION
Cholesterol synthesis enzymes are typically localized to the ER (endoplasmic reticulum), and as membrane-associated proteins, little has been done to determine their exact topology and structure. We are particularly interested in DHCR24 (3βhydroxysterol 24-reductase), the final enzyme in cholesterol synthesis, as it is a multifunctional protein, which is intricately regulated (reviewed in [1]). Initial studies characterizing DHCR24 by Greeve et al. [2], and subsequently confirmed by Wu et al. [3], found that it was targeted to the ER, and to a lesser degree, to the Golgi.
terol is dependent on FAD [7], suggesting functionality of the conserved domain. In addition, DHCR24 contains conserved binding domains for both p53 and Mdm2, required for mediating cellular responses to oncogenic and oxidative stress [3]. DHCR24 is reported to be proteolytically cleaved during apoptosis, at two highly conserved caspase recognition motifs located in the FAD-and p53-binding domains, presumably destroying them, producing a soluble, 40 kDa peptide [2]. How DHCR24 interacts with the membrane, however, is less well known.
Using simulated DHCR24 membrane models with and without substrate and cofactors, Pedretti et al. [5] predicted DHCR24 as a monotopic membrane protein; with the N-terminus partially embedded in the membrane as a stem or 'peduncle' (meaning 'little foot'), rather than traversing the membrane bilayer (bitopic). This 'peduncle' firmly anchors the protein to the membrane, with the C-terminus protruding into the cytoplasm, allowing access to substrates and cofactors [5]. The DHCR24 homologue, DWF1 in Arabidopsis thaliana, has a similar predicted structure, based on hydropathy predictions: strong membrane association and a cytoplasmic C-terminus [8].
However, Lu et al. [9] published the experimental work in support of a bitopic, single TMD (transmembrane domain) model of DHCR24: a lumenal N-terminus, followed by a single TMD (also at the N-terminus), and cytoplasmic C-terminus. The N-terminus was essential for DHCR24 membrane association, as deletion of this region (1-58, 58 DHCR24) translocated DHCR24 to the cytoplasm. The orientation of the N-and C-termini were deduced from protease protection assays of fluorescent fusion constructs; an N-terminal Ds-Red fluorescent fusion protein was preserved, indicating that it was lumenal, whereas a C-terminal green fluorescent protein was degraded by trypsin, indicating cytoplasmic localization. By contrast, Pedretti et al.'s 'peduncle' model [5] predicts that the N-terminus is buried in the membrane, and the Cterminus is cytoplasmic. However, fusion of the large (∼28 kDa), soluble, Ds-Red fluorescent protein to the N-terminus could feasibly draw the membrane-embedded N-terminus into the lumen, that is, interfere with the native membrane topology of DHCR24. Therefore we aimed to determine the validity of the two different published models [5,9], investigating how DHCR24 interacts with the ER membrane, and its membrane topology.

Cell fractionation
The cell fractionation protocol was adapted from Feramisco et al. [13]. Cells were washed and scraped into PBS on ice and centrifuged at 1000 g. The cellular pellet was resuspended in 400 μl Buffer A [10 mM Hepes-KOH (pH 7.4), 10 mM KCl, 1.5 mM MgCl 2 , 5 mM sodium EDTA, 5 mM sodium EGTA and 250 mM sucrose] and passed through an 18 gauge needle 50 times. Protein concentration of cell lysate was determined using the BCA assay. Cell lysate was centrifuged at 1000 g for 5 min at 4 • C. The 1000 g pellet was resuspended in 100 μl Buffer B [20 mM Hepes-KOH (pH 7.6), 25 % (v/v) glycerol, 0.42 M NaCl, 1.5 mM MgCl 2 , 5 mM sodium EDTA, 5 mM sodium EGTA], and incubated at 4 • C for 1 h with end over end mixing, and centrifuged at 100 000 g for 30 min at 4 • C. The supernatant from this centrifugation was collected and designated the nuclear fraction.
The supernatant from the 1000 g centrifugation was centrifuged at 100 000 g for 30 min at 4 • C. The supernatant was collected and designated the cytosolic fraction. The pellet was resuspended in 100 μl Buffer C [10 mM Tris-HC1 (pH 7.4), 100 mM NaC1, 1 % (w/v) SDS], and designated the membrane fraction.

Membrane isolation
To prepare membranes for differential solubilization and protease protection assays, cells were washed and scraped in PBS and centrifuged at 1000 g for 5 min at 4 • C. The cellular pellet was resuspended in Buffer A or D (Buffer A containing 100 mM NaCl; as indicated in figure legends), passed through an 18 gauge needle 50 times, and the cell lysate was centrifuged at 1000 g for 5 min at 4 • C. The supernatant was then centrifuged at 20 000 g for 15 min at 4 • C, and the resulting membrane pellet was resuspended in 65 μl Buffer A or D. Membrane protein concentration was determined using the BCA assay.

Differential solubilization
To determine the strength of the protein-membrane interaction, equivalent amounts of membrane were treated with 200 μl 1% SDS, Buffer A, 0.1 M Na 2 CO 3 (pH 11.5), or 1 M NaCl, and incubated at 4 • C for 30 min with end over end mixing. Treated membranes were centrifuged at 100 000 g for 30 min at 4 • C; the supernatant was designated the soluble cytoplasmic fraction (C), and the membrane pellet was resuspended in 200 μl Buffer D and designated the insoluble membrane fraction (M).

Protease protection assay
To determine the membrane orientation of the N-and C-termini, equivalent amounts of membrane were treated with trypsin as described [13]. Briefly, membranes were treated with the indicated amount of trypsin in the presence or absence of Triton X-100 in Buffer A, for 30 min at 30 • C. Reactions were stopped by the addition of loading buffer and heat inactivation at 95 • C for 10 min.

RESULTS
DHCR24 has an extremely hydrophobic N-terminus, suggesting a candidate region for TMDs To investigate the topology of DHCR24, we first used an in silico approach. TM (transmembrane) proteins are hydrophobic by nature, and can therefore be predicted based on this characteristic. The human DHCR24 sequence from UniProt (Q15392) was analysed for putative TMDs based on hydrophobicity by the method of Kyte and Doolittle ( Figure 1A) [15]. Two potential TMDs were identified at peak regions of hydrophobicity at the extreme N-terminus of DHCR24 (A, B; Figure 1A), with scores above the recommended threshold by Kyte and Doolittle (dotted red line, 1.6) [15]. As DHCR24 is not an extremely hydrophobic protein overall, having a GRAVY (grand average of hydropathicity) score of − 0.061 (blue line), other hydrophobic peaks are discernible above this value, but do not reach the Kyte and Doolittle threshold of 1. 6.
A more specific prediction of regions that lie within the ER membrane is the free energy calculation for membrane insertion by means of the Sec61 translocon ( G mi ): negative G mi values are indicative of potential TMDs of strong hydrophobicity and sufficient length; positive G mi values are less likely to be TMDs [24]. Both potential TMDs identified by the method of Kyte and Doolittle have a negative G mi , further indicating them as putative TMDs (A, B; Figure 1B). Three additional regions were identified as putative TMDs based on their G mi values (c, d, e; Figures 1A and 1B); however, their G mi values were positive, and these regions did not reach the threshold for detection of TMDs in the hydropathy plot. Therefore based on their hydrophobicity and their G mi values, they are less likely to form TMDs.

DHCR24 is not predicted to associate with the membrane via post-translational modifications
Hydrophobic PTMs (post-translational modifications) can facilitate membrane attachment. Common groups include the irreversible attachment of various lipid groups, such as prenyl, myristoyl or GPI (glycophosphatidylinositol) groups or the reversible attachment of a palmitoyl group. The prediction programs Myristoylator [18], NMT [19], big-PI predictor [20], PrePS [21] and CSS-Palm 2.0 [22] predict that DHCR24 is unlikely to associate with the ER membrane via PTMs; therefore, these means of membrane attachment were not considered for DHCR24.

Predicted membrane topology of DHCR24 consists of N-terminal TMDs
For advanced DHCR24 TMD prediction, TOPCONS and TM-HMM membrane topology prediction programs were employed. Both of these programs use HMM and/or neural network algorithms to integrate sequence information, and MSAs (multiple sequence alignments). TOPCONS uses multiple TMD prediction programs: SCAMPI-seq and SCAMPI-msa [25], PRODIV and PRO [26] and OCTOPUS [27]. The TOPCONS global prediction (TOPCONS consensus) [25] is a consensus of these five predictions, as well as the G mi [24] and ZPRED [28] algorithms. Both regions identified by their hydrophobicity and G mi as putative TMDs (A, B; Figure 1A) and the three regions with positive G mi values (c, d, e; Figure 1A) were predicted by at least one of the TOPCONS prediction algorithms ( Figure 1C); however, only the two N-terminal hydrophobic regions (A, B) were uniformly predicted by all five of the individual algorithms, and the global TOPCONS prediction (TOPCONS consensus; Figure 1C). Furthermore, the self-assessed reliability of the TOPCONS consensus is lowest in the region of the three other putative TMDs (c, d, e).
The TMHMM prediction [16], which also includes information on evolutionary conservation from MSAs, agreed with the TOPCONS consensus, predicting two N-terminal TMDs with high probability ( Figure 1D). However, TOPCONS predicted a cytosolic C-terminus, whereas TMHMM predicted it to be lumenal.

A secretory SP is predicted for DHCR24
DHCR24 contains a putative TMD at its extreme N-terminus (A; Figure 1), which is usually the location of a secretory sequence. .    Predicted TMDs based on the apparent free energy difference for ER membrane insertion by the Sec61 translocon ( G mi ) below the 1.6 threshold are marked in blue (c, d, e) (B) Putative TMDs and their probability based on G mi [24]. (C) TMD predictions by the five TOPCONS algorithms: SCAMPI-seq and SCAMPI-msa [25], PRODIV and PRO [26] and OCTOPUS [27] and the TOPCONS global prediction (TOPCONS consensus) [25]. Potential TMDs are boxed grey (cytosol → lumen) or white (lumen → cytosol). Dotted grey lines indicate DHCR24 truncations, which are used in  (D) DHCR24 membrane topology as predicted by the TMHMM server [16], with the probability of cytosol, membrane or lumenal location in red, grey or blue, respectively.
Signal sequences allow translocation of a protein to, or across, the ER membrane, which is cleaved for targeting beyond the ER membrane (SP), or retained within the membrane for ER membrane localization [SA (signal anchor)] [29]. Both SPs and SAs contain a stretch of hydrophobic residues, and therefore can be incorrectly predicted as TMDs by many prediction algorithms [30]. Most TMD prediction programs, such as TMHMM cannot discriminate between N-terminal TMDs and SPs and SAs [30].
Prediction algorithms, such as SignalP, use a neural network to distinguish N-terminal TMDs from SPs based on a combined scoring system for each position at the N-terminus [31]. The most N-terminal putative TMD for DHCR24 (A; Figure 1A) was predicted to be an SP (Supplementary Figure S1; available at http://www.bioscirep.org/bsr/034/bsr034e098add.htm), with the cleavage site at 22-23 based on all scoring systems: predicted cleavage site (C-score), SP length (S-score), and the combined cleavage site score (Y-score; a combination of the Cscore and the slope of the S-score), which is a better cleavage site prediction than the raw C-score alone.
To test whether the putative SP is required for ER membrane localization, truncated DHCR24 ( 23 DHCR24) was transfected into CHO-7 cells, and the cellular localization examined. Through cell fractionation, 23 DHCR24 was smaller than WT DHCR24 and, like WT DHCR24, 23 DHCR24 localized primarily to the membrane fraction (results not shown). This demonstrates that the putative SP is not cleaved off in WT DHCR24, and is not necessary for membrane retention.
The discounted SP and TMD lie in one of the least conserved regions of DHCR24. The C-terminal region, however, is highly conserved and contains multiple protein-binding/recognition sites for FAD [4][5][6], caspase 3 [2], p53 and Mdm2 [3]. The three predicted TMDs that did not have negative G mi values all reside close to functional domains of DHCR24 (136-156, FAD; 210-230, Mdm2; 349-369, p53), that by definition should be cytosolic. Furthermore, a single TMD model is fitting, as it would allow the C-terminus to reside in the cytoplasm, making the multiple binding sites in DHCR24 accessible to soluble-binding partners. .

Membrane orientation of the N and C termini of DHCR24
To determine the localization of the N-and C-termini of DHCR24 with respect to the ER membrane, CHO-7 cells were transfected with DHCR24 with a V5 epitope tag at either the Nterminus (V5-DHCR24) or C-terminus (DHCR24-V5), and cotransfected with Insig, an integral membrane protein with a myc epitope on the cytoplasmic C-terminus (Insig-myc) [13] (Figure 2A). Intact ER membrane vesicles were isolated and treated with trypsin in the presence or absence of Triton X-100, subjected to SDS-PAGE and immunoblotted for V5 ( Figure 2A). As controls, samples were also immunoblotted with myc for Insig, and R139 for Scap, which recognizes a lumenal loop (residues 540-707) in hamster Scap [32] (Figure 2A). Therefore when trypsin was added to microsomes, the myc epitope on Insig was digested, and the R139 epitope of Scap remained intact, as it is located within the lumen (Figure 2A, lanes 1-4). Similar to In-sig, DHCR24-V5 was digested. However, V5-DHCR24 remained stable (Figure 2A, lanes 1-4). When Triton X-100 was added to partially solubilize the membrane, allowing trypsin access to lumenal peptides, the R139 epitope of Scap was digested. However, V5-DHCR24 remained relatively stable (Figure 2A, lanes 5-8), indicating that the N-terminus was protected from trypsin digestion by the membrane.
Although there was no obvious change in the size of the V5-DHCR24 band with trypsin and Triton X-100 treatment, the band intensity slightly decreased (Figure 2A, lanes 5-8), and a small fragment was discernible with increasing amounts of Triton X-100 (∼11 kDa, Figure 2B, lanes 6-9). This smaller band was not visible in non-treated conditions (lanes 1-4) or when trypsin was inactivated with SDS (lane 10). Furthermore, the band appeared with or without Triton X-100 treatment, indicating a partially accessible cytoplasmic loop following a membrane embedded N-terminus. The calculated molecular   weight of this digested fragment indicates that the loop is located between residues 70-100, which contains six possible arginines/lysines that could be cleaved by trypsin (dotted red line; Figure 2B). Altogether, these data are indicative of a cytoplasmic C-terminus and a non-accessible N-terminus, which is either embedded within the membrane, or lumenal but inaccessible, due to strong interactions with the ER membrane ( Figure 2).

The putative TMD is not essential for DHCR24 membrane association
To determine if the hydrophobic N-terminus containing the putative TMD of DHCR24 is required for ER membrane attachment, truncated DHCR24 ( 56) was transfected into CHO-7 cells, and the cellular localization examined. DHCR24-V5, like Insig-myc, localized only to membranes ( Figure 3A). The truncation lacking the putative TMD, 56 DHCR24-V5, also localized specifically to membranes ( Figure 3A). Therefore contrary to the predictions in Figure 1, these very hydrophobic regions are not necessary for DHCR24 membrane association.

DHCR24 associates strongly with membranes
To determine if DHCR24 is an integral membrane protein, differential solubilization was performed on membranes to determine the manner in which DHCR24 associates with it. The ability of buffers (detergent, high pH, high salt) to dissociate DHCR24 from the membrane (M) into the cytoplasm (C) indicates the strength of its association with the membrane (ionic or hydrophobic).
When exposed to aqueous buffer, DHCR24 was located in the membrane fraction, but was released into the cytoplasmic fraction upon treatment with 1 % (w/v) SDS ( Figure 3B). Treatment with 0.1 M Na 2 CO 3 or 1 M NaCl, reagents known to disrupt peripheral associations by altering pH and electrostatic interactions, respectively, were not able to solubilize DHCR24, and it remained in the membrane fraction ( Figure 3B). For all treatments, Insig-1 displayed similar results, indicating DHCR24 is associated with membranes, and that it is an integral membrane protein ( Figure 3B). Similarly, 56 DHCR24-V5 behaved like an integral membrane protein ( Figure 3B). Although 56 DHCR24-V5 remained in the membrane fraction after treatment  with 0.1 M Na 2 CO 3 , the presence of a weaker signal in the cytoplasmic fraction indicated that 56 might be solubilized with this treatment. However, performing the treatment and centrifugation process a second time showed that 56 DHCR24-V5 remained membrane associated (Supplementary Figure S2; available at http://www.bioscirep.org/bsr/034/bsr034e098add.htm), suggesting that regions beyond the first 56 residues are integral for strong membrane association.
Other candidate TMDs are not essential for DHCR24 membrane association Having established that the most likely TMDs (A and B, Figure 1) was not essential for membrane localization, we next examined the other predicted TMDs. Based on the G mi values, and the TMHMM probability score (Figure 1), the two next best TMD candidates were examined (c, 136-156; d, 210-230). Truncated versions of DHCR24 ( 160, 240; Figure 4A) were transfected into CHO-7 cells, and the cellular localization examined. Similar to WT DHCR24-V5, both truncations were found primarily in the membrane fraction ( Figure 4B), demonstrating that the first 240 residues are also not necessary for membrane attachment. The appearance of secondary bands (∼22 and ∼32 kDa) in the 240 DHCR24 membrane fraction is indicative of proteolytic fragments, which are also associated with membranes, and would encompass the untested putative TMD (e). This further suggests that there are more membrane-binding site(s) closer to the Cterminus. To confirm that the truncations are not mislocalized or aggregating, immunofluorescence microscopy was performed. WT and 240 DHCR24 both localized predominantly to the ER ( Figure 5), indicating that DHCR24 remains ER membrane associated despite N-terminal deletions.

DISCUSSION
To further elucidate the membrane topology of DHCR24, we predicted how it interacts with the ER membrane through the use of numerous membrane insertion algorithms, as well as considering the current literature on DHCR24 topology. We then characterized this interaction using biochemical methods.
To create a DHCR24 topology model, we first characterized DHCR24 based on its hydropathy using the method of Kyte and Doolittle [15]. Two regions of peak hydrophobicity were identified as candidates for TMDs ( Figure 1A), which were confirmed by numerous TMD prediction algorithms ( G mi ; Figure 1B) and programs (TOPCONS and TMHMM; Figures 1C and 1D). DHCR24 was predicted to contain an SP (Supplementary Figure  S1), which generally targets proteins to or across the ER membrane, and this was located at the first putative TMD. The efficacy   The predicted and published models for DHCR24 topology are presented with features indicated, compared with our hypothesized DHCR24 membrane topology model (the present study).
of a DHCR24 SP was tested by the creation of a truncation construct, which did not contain the putative SP. This was smaller than its WT DHCR24 counterpart, and still associated with membranes. Secondly, the N-terminal V5 tagged DHCR24 construct (V5-DHCR24) was the same size as C-terminal DHCR24-V5 (V5-tagged DHCR24), indicating no cleavage of the N-terminus, which would be predicted for a true SP. These data show that the N-terminus is not cleaved, nor is it required for membrane targeting, demonstrating that DHCR24 does not contain an N-terminal SP. This is similar to cytochrome P450, which is also an integral ER membrane protein that does not require an N-terminal SP for ER targeting, but contains complex N-and C-terminal retention signals [33]. Upon deletion of the N-terminus (including the putative TMDs; 56), DHCR24 still associated with the membrane (Figure 3A), and could not be solubilized by various treatments (Figure 3B). Again, this is similar to cytochrome P450, which contains a hydrophobic N-terminus predicted to associate with the ER membrane, which retained strong membrane attachment after truncation experiments and differential solubilization [34,35]. Further N-terminal truncations in DHCR24 that deleted lower scoring TMD candidates also associated exclusively with the

Figure 7
Hypothetical membrane topology model of DHCR24 This model integrates the known PTM sites to refine our working model of DHCR24 membrane topology. As the focus of this work has been on the N-terminus, this hypothetical model does not attempt to encapsulate the membrane topology of the C-terminus, which is likely to snorkel in and out of the membrane similar to the N-terminus, and is therefore protected from trypsin digestion ( Figure  membrane ( Figure 4). These data discount the findings of Lu et al. [9] that the N-terminus is essential for membrane attachment, and that deletion of the N-terminus ( 58) results in altered cellular localization of DHCR24. In contrast, our findings demonstrate that the membrane association of DHCR24 is not reliant on TMDs at the N-terminus, but contains membrane associated region(s) beyond the hydrophobic N-terminus. We also elucidated the membrane orientation of the N-and C-termini of DHCR24, for which the TMD prediction programs gave different results. Using protease protection assays, we determined that the C-terminus was cytoplasmic, due to its accessibility to trypsin. The N-terminus, however, was strongly protected, even when trypsin was accessible to the lumen by partial solubilization of the membrane (Figure 2A). This suggests that the N-terminus is 'hidden' from trypsin; either due to being embedded within the ER membrane, which would support (in part) Pedretti et al. 's 'peduncle' model [5], or due to peripheral interactions with the ER membrane. Both of these possible models disagree with the findings of Lu et al., that the N-terminus is lumenal, which is most likely due to the difference in epitopes used in the protease protection assay [9]. As mentioned, Lu et al. [9] used a large fusion protein (∼28 kDa), whereas we used a small V5 epitope (∼1 kDa). Testing the membrane association of a C-terminal truncation using a differential solubilization assay would determine the type of interaction the N-terminus has with the membrane (integral or peripheral). Finally, the identification of a digested N-terminal ∼11 kDa fragment suggests that the membrane associated N-terminus is fol-lowed by a sterically hindered cytosolic loop at position ∼70-100 ( Figure 2B).
Overall, our findings demonstrate that DHCR24 has a strong affinity with the ER membrane, due to multiple membraneassociated regions, including beyond the N-terminus ( 240; Figure 4). The extreme hydrophobic nature of the N-terminus suggests that it is also membrane associated, but this needs to be confirmed experimentally by C-terminal truncations. We hypothesize that there are multiple membrane associated regions throughout the protein, or re-entrant loops that associate with only one leaflet of the membrane. Re-entrant loops are typically rich in small residues such as glycines and alanines [36], and one of the lower scoring putative TMDs (d, Figure 1) is rich in these amino acids, and therefore a prime candidate for a re-entrant loop. Our proposed model may also contain classical TMD(s), as although this study ruled out the most likely putative TMD, one low-scoring putative TMD was untested (e, Figure 1). As this region is beyond the first 240 residues, it is likely to explain the strong membrane localization of 240 DHCR24 (Figure 4). Our data strongly suggest that relocation of DHCR24 to other cellular compartments, proposed to occur in response to stress signals, is highly unlikely [3]. This conclusion should be re-examined in light of our results.
Altogether, our findings are at odds with prediction algorithms and published topology models ( Figure 6, [5,9]   hydrophobic domains passing through and/or anchoring DHCR24 to the hydrophobic core of the lipid bilayer of the ER. To create a hypothetical membrane topology model of DHCR24, we have also utilized in addition to our experimental findings our knowledge of known PTMs on DHCR24 [37]. Phosphorylation and ubiquitination sites must be accessible to cytosolic kinases and ubiquitin ligases respectively, and therefore we have used these data to refine our topology model of DHCR24 (Figure 7). To the best of our knowledge, we are the first to propose that the expanding databases of PTMs provide an unlikely but invaluable resource for membrane topology mapping, including for other atypical cholesterogenic membrane proteins [38]. After all, this is not the first example of a cholesterol biosynthetic enzyme that has an atypical membrane topology. Lanosterol synthase [also known as OSC (oxidosqualene cyclase)] is a monotopic ER membrane protein, containing a membrane associated region, which resides in one leaflet of the membrane, and therefore does not span the bilayer [39]. This hydrophobic catalytic domain resides in the membrane, which allows access to the hydrophobic substrate, oxidosqualene. As the substrate-binding domain of DHCR24 is as yet uncharacterized, and desmosterol is also hydrophobic, a membrane associated catalytic site is likely.
To further refine our membrane topology model of DHCR24 (Figure 7), additional experimentation is required to precisely define membrane, cytoplasmic, and lumenal regions, including whether the final putative membrane region (e) is bona fide. This commonly requires modification of the specific residues or introduction of epitopes and recognition sites. Cysteine derivitization and glycosylation site mapping have previously been successful in elucidating the membrane topology of membrane bound proteins, such as Insig [13], and these techniques could similarly be applied here.

AUTHOR CONTRIBUTION
Eser Zerenturk, Laura Sharpe and Andrew Brown designed the experiments. Eser Zerenturk and Laura Sharpe performed the experiments. Eser Zerenturk, Laura Sharpe and Andrew Brown analysed the data. Eser Zerenturk, Laura Sharpe and Andrew Brown wrote the paper.


extfieldsofstudy: Biology Medicine

text:
An update of the goat genome assembly using dense radiation hybrid maps allows detailed analysis of evolutionary rearrangements in Bovidae

Background The domestic goat (Capra hircus), an important livestock species, belongs to a clade of Ruminantia, Bovidae, together with cattle, buffalo and sheep. The history of genome evolution and chromosomal rearrangements on a small scale in ruminants remain speculative. Recently completed goat genome sequence was released but is still in a draft stage. The draft sequence used a variety of assembly packages, as well as a radiation hybrid (RH) map of chromosome 1 as part of its validation. Results Using an improved RH mapping pipeline, whole-genome dense maps of 45,953 SNP markers were constructed with statistical confidence measures and the saturated maps provided a fine map resolution of approximate 65 kb. Linking RH maps to the goat sequences showed that the assemblies of scaffolds/super-scaffolds were globally accurate. However, we observed certain flaws linked to the process of anchoring chromosome using conserved synteny with cattle. Chromosome assignments, long-range order, and orientation of the scaffolds were reassessed in an updated genome sequence version. We also present new results exploiting the updated goat genome sequence to understand genomic rearrangements and chromosome evolution between mammals during species radiations. The sequence architecture of rearrangement sites between the goat and cattle genomes presented abundant segmental duplication on regions of goat chromosome 9 and 14, as well as new insertions in homologous cattle genome regions. This complex interplay between duplicated sequences and Robertsonian translocations highlights the rearrangement mechanism of centromeric nonallelic homologous recombination (NAHR) in mammals. We observed that species-specific shifts in ANKRD26 gene duplication are coincident with breakpoint reuse in divergent lineages and this gene family may play a role in chromosome stabilization in chromosome evolution. Conclusions We generated dense maps of the complete whole goat genome. The chromosomal maps allowed us to anchor and orientate assembled genome scaffolds along the chromosomes, annotate chromosome rearrangements and thereby get a better understanding of the genome evolution of ruminants and other mammals. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-625) contains supplementary material, which is available to authorized users.

Background
The domestic goat (Capra hircus) is an important source of food and wool (cashmere) especially important in low input production systems. The ability of goat to consume a wide range of vegetation explains its significance for the agricultural economy. In 2007, the FAO estimated the world population of goats to be around 800 million animals. Despite this economic importance, goat genomic resources are not extensive compared to other livestock species, such as cattle, pig, chicken, or sheep. Goat is a ruminant and a member of the Caprinae order, a clade distinct from Bovinae in the Bovidae family. Its last common ancestor with cattle is dated between 19 and 40 million year (Myr) ago, and with sheep between 6.2 and 8.5 Myr ago [1]. Comparative mapping between cattle and goat using fluorescent in situ hybridization (FISH) showed a strong synteny conservation between the two genomes with a few large-scale rearrangements, including the well-known 9-14 translocation differentiating Bovinae and Caprinae [2] and differences in intra-chromosomal organization of the X chromosome [3,4]. However, the history of genome evolution and chromosomal rearrangements on a small scale in ruminants remain speculative, mostly due to limited number of mapped loci.
High-resolution gene maps provide basic but crucial linear information on the physical organization of the genome and typically serve as the backbone for further genomic research. For example, robust information about marker organization on chromosomes is crucial for linkage analysis and marker assisted animal breeding. Although assembling next-generation sequencing data into a draft genome comprising scaffolds is relatively straightforward, constructing a physical map of the structure of chromosomes is still difficult and costly. Recent years have witnessed a slow increase in gene mapping data in the goat since the first genetic and cytogenetic maps [2,5]. To date, approximately only 550 loci have been mapped to the goat genome using linkage maps of low resolution [6]. The initial release of the goat genome sequence is based entirely on short-read de novo sequencing of a Yunnan Black goat that yielded a~2.66 Gb assembly of 284,683 scaffolds (N50 = 3.06 Mb). Anew improved optical mapping technology was used to joined 2,090 scaffolds into 315 larger super-scaffolds (N50 = 16.3 Mb). Finally, super-scaffolds were assembled into chromosomes based on conserved synteny with cattle [7]. The goat genome sequence could yield a wealth of information about gene content and putative regulatory elements sequences; however, it is still lacking long-range continuity and its released form only gives a fragmented view of the genome.
Another genomic resource, radiation hybrid (RH) maps, played a pivotal role in the process of mapping animal genomes and validation of draft assemblies of the genome sequence, especially in mammals [8][9][10][11]. For obtaining dense and accurate next-generation RH maps, new methods were developed recently exploiting the increased availability of comparative genomic resources [12] and have proved successful for the production of RH maps based on SNP array genotyping data [13,14].
As a part of the goat genome project, we built RH maps of chromosome 1 that helped to validate the draft genome assembly. In this study, using a goat RH panel that we generated recently [15], we constructed robust physical maps of 45,953 loci for all of goat chromosomes except chromosome Y. The whole-genome goat RH map, as an independent chromosomal map, allowed us to anchor and orientate assembled genome scaffolds along the chromosomes thus permitting the annotation of chromosome rearrangements and the study of their evolutionary history using ruminants and other mammalian genomes.

Results and discussion
Genotype calling of SNPs in the RH clones In building the RH maps, we sought to provide accurate, detailed and reliable physical maps at the whole-genome scale to aid in the validation and improvement of the goat genome assembly. Genotyping SNP arrays from other species on the goat RH panel turned out to be a successful strategy for the genotype calling of RH clones. A similar strategy of marker selection from close related species had been applied for constructing buffalo RH maps [16]. Because sequence conservation between the three species considered here is high, a number of bovine and sheep SNPs presented a positive signal in the goat. Moreover, the SNPs that did not provide positive signals were particularly useful for the genotype calling procedure as they gave an internal control, within each clone, of the signal exhibited by non-retained SNPs. We believe it highly increased the robustness of the genotype calling procedure, which is ultimately proven by the quality of the RH maps produced and for future studies we would recommend this approach. The feasibility depends on the availability of a SNP array in a closely related species. To test that whether the approach can be carried out in a species, a simple experiment would consist in genotyping of whole genome sample of the species of interest with the related species array and check that (1) a number of SNPs can provide signal and (2) a number of SNPs do not provide signal. Optimally both the categories should be in similar proportion on the array. Recall that the fact that SNPs are actually polymorphic in the species of interest is not relevant for RH mapping.
The goat radiation hybrid panel [12] including 94 hybrid clones was genotyped using IlluminaBovineSNP50K BeadChip and OvineSNP50K BeadChip. RH vectors were constructed for more than 110,000 SNPs out of which 54,318 could be assigned a position on the chromosomes of the goat genome assembly CHIR_1.0. We performed genotype calls (Figure 1) independently for the sheep and the bovine SNP arrays (see Methods). Table 1 summarizes the results of our genotype calling procedure for each SNP array. The two analyzes provided consistent estimates of the panel retention fraction: 37% based on the sheep array data and 35% based on the bovine array data. False negative rates and proportion of missing data were also similar. The call rate was slightly higher for the bovine arrays (80%) than for the ovine array (73%), however, the number of SNPs that could be called was slightly lower (35,521 vs. 43,628). This is explained by the fact that there are much more SNPs (12,577) on the bovine array which do not exhibit a positive signal in the whole genome goat sample than that on the sheep array (3,969), consistent with the phylogeny of these three species.
Even low levels of error result in large losses of information about breakage probabilities, markedly increase the uncertainty in the marker ordering, and inflate estimates of inter-marker distances and total map length [17]. The high error rate of genotypes for constructing certain maps prevents from recovering sufficient signal for ordering the markers. The ability to produce useful robust maps in comparative RH mapping is impaired when the error rate in the data increases above 10% [14]. Before constructing the integrated robust map, we also checked the vector data from bovine SNPs and Ovine SNPs, respectively, and constructed two robust maps for each chromosome. These maps (data not shown) presented a very low error rate in the RH data, in consistent with the low false negative rate estimated by the Q-value R package [18]. However, we note that the slight differences in experimental conditions between genotyping RH clones on the BovineSNP50K array and on the OvineSNP50K array. It could bring certain systematic errors when merging the two datasets, which we cannot assess, since our RH vectors have not been validated by replicated genotyping. However, the availability of negative control SNPs within each clone offers a replication across clones. For increased robustness, we applied three approaches to remove possible genotyping errors before constructing the integrated robust maps.

Genome map construction and characterization
Briefly, for each chromosome, the analysis of RH mapping was done in three steps:  Imax for SNP marker genotyping. We called genotypes of clones at SNPs based on the raw signal intensities obtained from Illumina genotyping system. A number (Nneg) of SNPs showed very low intensities for both alleles in the goat (panel A). Specifically, the statistic used for genotype calling of a given SNP in a given clone was the maximum observed intensity over the two possible alleles (herein called Imax). The distribution of Imax in a given clone can be seen as a mixture of two underlying distributions, where one for the non-retained SNPs in the clone and one for the retained SNPs (e.g. panel B).
The final robust RH map contained 45,953 SNP Markers spanning totally 76,800 cR5000. Marker information of the RH maps is shown in Additional file 1. These markers are mapped at 17,628 distinct positions (38.4% of the total number of markers) hence the maps are highly saturated. A summary of the maps including chromosome size, RH map length, and the number of SNPs kept at each step for each chromosome is shown in Figure 2. The number of SNPs removed during construction of robust maps was relatively constant across chromosomes, with the notable exception of CHI23 (27.6%), CHI24 (25.6%), and CHI29 (22.1%), for which a large portion in the middle of the chromosome could not be considered in the robust maps. Marker density was also fairly uniform between chromosomes (results not shown).
Detailed characteristics of the final RH maps are provided in Additional file 2. In general, the retention frequency was stable across the genome, with the expected minor increases at the distal chromosomes. The RH maps uniformly cover >99% of CHIR_1.0 and have an average resolution of 32.6 kb/cR5000 based on a goat genome size of 2.66 GB. The map of chromosome 19 has the shortest length of 1,228 cR and does not cover the region from 40.9 Mb to 62.1 Mb in CHIR_1.0. On CHI19 we retrieved 2,377 RH vectors from raw genotyping data before any quality control. In the raw data, the retention from 40.9 Mb to 62.1 Mb is extremely high because it surrounds the selection locus TK1 (thymidine kinase 1 gene). We further observed that the regions carry TK1 gene at position of 52.7 Mb, where the anticipated highest retention near 100%. Consequently, the markers in this region could not be mapped ( Figure 3). The log of the ratio of the posterior probability of the best map to the second best map (similar to a LOD score) ranged from 0.17 to 3.54 for robust RH maps. The two chromosomes with the lowest LOD value are CHI21 and CHI4 possibly indicating a higher genotyping error rate in their RH vectors. Genotyping errors are also a vital factor to increase the distance of inter-markers unexpectedly. The total map length of CHI2, CHI7 and CHI18 are higher than expected given the sequence length of these chromosomes ( Figure 2). To further characterize the resolution of our RH maps, we studied the distribution of the length of retained fragments using a methodology suggested in a previous study [19]. We found a generally good agreement with the expected distribution of segment lengths (Figure 4), except for the class of very small segments (<2 cR) which is completely lacking in our maps. This is most likely because we are not able to order such close markers with confidence and therefore these small segments are not saw in the robust RH maps. We estimate that the RH panel can order confidently markers separated  by at least 65 Kb, which is smaller than the 90 Kb average distances between genes (considering 22,175 genes on a genome of 2.66 Gb).
The main drawback of current high throughput sequencing technologies is the short length of the reads they produce. Many of the novel assemblies including the current goat draft assembly obtained with these approaches are composed of a very large number of scaffolds. This fragmentation does not affect gene discovery, polymorphism analysis and sequence comparison between closely related species but it greatly limits the study of genome structure and chromosomal evolution. The comparative RH approach incorporated the prior knowledge of marker orders in CHIR_1.0 into the statistical model. It has the key property that a different ordering than that of the reference order will be accepted as an alternative order only  when the RH data is informative enough. As a consequence, the robust maps allow pinpointing the regions where the assembly order disagrees with an RH map order that is strongly supported by the RH data, so this approach allowed improvements over classical a framework-based mapping strategy [19]. The RH maps presented here offer a powerful complement to the available genetic/cytogenetic maps with very limited resolutions [2,6,20] and the genome sequence map [7]. Our RH panel and RH maps also revealed the potential ability to map the unplaced genes or scaffolds and to check possible misassembled regions in further research.

Validation of the draft genome assembly and anchoring scaffolds on chromosomes
The resolution of discrepancies directly addresses the question of the reliability of the order defined on one side by the genome map and on the other by the assembly. The construction of robust maps was precisely designed to address the reliability of RH maps [14]. In the process of assembling CHIR_1.0, a first draft was produced by standard in silico assembly pipelines yielding numerous scaffolds, followed by an automated optical mapping procedure to assemble scaffolds into super-scaffolds and finally conservation of synteny with cattle was assumed in unassembled regions to organize super-scaffolds into whole chromosomes. For assembling CHIR_1.0, we had constructed two RH maps of chromosome 1 using bovine SNPs and ovine SNPs separately. CHI1 is the longest chromosome and carries many long super-scaffolds. These two maps suggested that the assemblies of super-scaffolds and chromosomes in CHIR_1.0 were accurate [7].
In this study, we merged all the candidate markers from bovine SNPs and ovine SNPs to construct comprehensive RH maps for the 29 autosomes and the X chromosome. The RH marker positions were compared to their sequence positions, providing a comprehensive link between the two physical maps. Based on this link, we aimed at resolving the discrepancies between the two maps to possibly improve the goat genome sequence. Our approach was to consider contigs and scaffolds construction (a small-scale assembly problem out of range given the RH map resolution) as accurate but to allow modifying their organization along chromosomes. More specifically, to resolve discrepancies and update the goat genome sequence using RH maps, we set the priority, from highest to lowest, to the scaffolds, the RH map, the superscaffolds and the cattle genome successively. Doing so we consider (1) the scaffold data provide ultimate resolution in local region of genome in comparison of any other map, and the assembly of goat scaffold is comprehensively correct at the local level, (2) while the automated whole-genome optical mapping efficiently generated supper-scaffold data for facilitating the anchoring of scaffolds onto chromosomes of CHIR_1.0, the quality control of the maps have not fully been verified, and (3) assuming conserved synteny with cattle should only be required when no other information is available.
Markers on the 30 robust maps can be aligned to1,910 scaffolds total 2.47 Gb of the goat genome sequence (Additional file 1). In parallel, the optical mapping for assembling CHIR_1.0 joined 2,090 scaffolds into 315 larger super-scaffolds (N50 = 16.3 Mb). Therefore, these two maps provide together a mean to anchor scaffolds to the whole goat genome. Based on these comparisons we listed a set of proposed improvements to the goat genome assembly, so the modifications of the preliminary assembly (CHIR_1.0) proposed in the present study only involved reordering of scaffolds along chromosomes. In our improved version of the assembly (CHIR_1.1), the positions or/and orientations of 32 scaffolds in CHIR_1.0 were adjusted based on RH maps. 24 scaffolds within 7 super-scaffolds (super3, super240, super102, super152, super192, super93, and super210) and 8 disjunctive scaffolds were reassigned based on the RH maps. As an illustration of such an update, Figure 5 represents assembly order in a region of CHIX in the CHIR_1.0 and CHIR_1.1. For all 30 chromosomes, we provide detailed pictures illustrating comparisons of RH maps with CHIR_1.0 and with CHIR_1.1 (Additional file 3). Our RH maps recovered most of the cattle synteny information, and reassessed the order and orientation of the scaffolds. As explained above, we did not try to correct sequencing data (contigs and scaffolds); however, our approach could potentially contribute to the identification of chimeric scaffolds.
For most markers considered here we were able to assign a position on the genome sequence but still suggest that gaps remain in CHIR_1.0. In chromosome 5 (89.1 Mb between scaffold740 and scaffold666), we underlie a region spanned 72.6 cR haven't aligned to any goat scaffold data, which mean that a continuous region of over 4.2 Mb missed there in do novo assembling of short sequence. Based on predictions using the 17-mer method, the current genome is~91% of the esti-mated~2.92 Gb size. Filling gaps in the goat genome will be a challenge for both technologies of genome mapping and genome sequencing. Since, all maps have some errors in them then it will be important to integrate information from various sources to allow for these checks to be consistent, given that they are generated essentially independently of one another. An approach to resolve contradicting orders is the exploitation of additional and independent source of information. (1) Inconsistency in marker order between the RH maps and other maps is likely due to inaccuracies in independently mappings. (2) Some of the remaining inconsistencies could be biologically grounded, reflecting breed structural variations. Indeed, the reference sequence and the RH panel were generated using the DNA from different individuals and from different breeds (a Yunnan black for the reference sequence and Boer samples for the RH panels), as has already been evidenced in the goat [20].

Mapping of unplaced scaffolds on the goat genome
As an application of the RH maps, we predicted the positions of SNPs with no localization on the assembly CHIR_1.0. Herein we call these unmapped SNPs "uSNPs". On the assembly CHIR_1.0, 391 uSNPs assigned to unplaced scaffolds longer than 10Kbhad RH vectors. Firstly, the most likely chromosome was determined based on the linkage with mSNPs and then the best possible position was determined in the LKH map. We excluded uSNPs on tiny scaffolds from map construction because of the limited resolution of our RH maps (we noticed that essentially all SNPs assigned to a given scaffold had consistent predicted localization on the same RH map). We identified 32 uSNPs (2.2%) on the X chromosome, more than expected compared to autosome maps. This could imply that for CHIX more genome sequence is missing in the current assembly than for other chromosomes. We set the same principle for mapping unplaced scaffolds as we used to validate the draft genome: unplaced scaffolds could only be inserted between two contiguous placed scaffolds and we aligned the candidate scaffolds and their neighbors based on the cattle genome order.
Localizing unplaced scaffold280 (1,268 Kb with 9 uSNPs) and scaffold767 (185 Kb with 3 uSNPs) is illustrated on Figure 6. The predicted localization for this scaffold, based on SNPs predicted positions, was at the beginning of CHI4. We verified that scaffold280 and scaffold767 aligned on the putative syntenic region of the cattle genome and constructed a local RH map of the region including uSNPs of both scaffolds. We found that both scaffolds were indeed strongly linked to scaffold641 on CHI4 so that both can be placed on the goat genome confidently. Another example, illustrated in Additional file 4, revealed that RH maps can be used to improve the goat sequence map accurately in a local region by both correcting the ordering of scaffolds and adding unplaced scaffolds.
Among the uSNPs for which we could assign a position on the goat genome, we were able to align 95 scaffolds present in the goat genome sequence but currently unassigned. Of these scaffolds 52 only could be assigned a position but not orientated, either because they carried a single uSNP or because they contained multiple tightly linked SNPs whose map order was questionable due to insufficient radiation-induced breakpoints. The orientations of the 53 scaffolds were obtained using comparative information with the cattle genome (UMD3.1) assuming conserved synteny. In total, 19.45 Mb of previously unplaced genome sequences were incorporated into the new version of the genome assembly. The comparison between the RH map with super-scaffolds and the CHIR_1.1 assembly are presented in Additional file 5. Of 1,910 scaffolds that can be aligned to the 30 robust RH maps, 1,884 were used to reassembling chromosome sequences. Only a few rearrangements exist between RH map and the new version of genome assembly. The updated assembly is available as an agp file (the CHIR_1.1 assembly, at http://goat.kiz.ac.cn)

Genome rearrangements between caprine and bovine genome
We used our RH updated genome assembly to identify conserved segments between species that led to novel discoveries of ancient chromosome rearrangements and study the sequence features of evolutionary breakpoints. Mapping a large number of markers on the goat genome and cross-referencing 96.4% of these with the map locations for the markers in the bovine achieved this alignment (Additional file 1). Overall, there is a very good colinearity between the goat RH maps order and the bovine assembly order, consistent with a relatively   short divergence time between the two species. However there are few discrepancies (Table 2), in agreement with chromosome banding features [21]. X chromosomes appear to have been subjected to large number of rearrangements, compared to the autosomes, suggesting that rearrangements might better tolerated on sex chromosomes. The X-chromosome linkage group is usually conserved in placental mammals and does not interplay with autosomes in chromosome evolution. Comparison of the cytogenetic maps of caprine and bovine X chromosomes also shows some changes in loci order, defining 6 common chromosome segments [3,22]. Both the RH map and the sequence map confirmed the well-known chromosome 9-14 translocation, with CHI14 containing a small BTA9q11-q13 segment including genes COL9A1 and ETH225 [23]. Our comparative mapping did not detect inversions between CHI2 and BTA2, or CHI19 and BTA19, which were suggested by these previous cytogenetic maps. We believe the relative inaccuracy in these cytogenetic maps may account for these differences. We characterized the sequence architecture of the evolutionary breakpoints for two ruminant rearrangements, one between CHI9 and CHI14 ( Figure 7) and one in CHI13 (Additional file 6). While the CHI9-CHI14 translation is a known rearrangement between ruminant genomes [6], the CHI13 one was not described previously. Cattle shares a large region of conserved synteny with the human and dog genomes at position 13.0 Mb on cattle BTA14, whereas goat presents species-specific breakpoints and loaded to two centromeres of CHI9 and CHI14 (Figure 8), so we considered that the rearrangement occurred in the goat lineage. Breakpoint intervals from two neighboring synteny blocks between cattle and goat reveal a mosaic of new insertions at rearrangement sites, rather than a simple model of nonhomologous end-joining (NHEJ). In Figure 7, Representative repeat [δ] Genes with colors denotes transcriptional orientation (the white represents "+" and the black represent "-").
[ε] Names of ANKRD26 homologsare marked in red. Breakpoint intervals were highlighted using mixture colored regions from two neighboring synteny blocks between cattle and goat, which revealed mosaic new insertions at rearrangement sites. Self-comparisons showed that these two centromeres share larger homologous region in cattle (panel A), different from in goat (panel B) suggest two cattle centromeres underwent a much longer evolution process than two goat centromeres. The breakpoint in goat genome (panel A) presents more complex architecture than one in cattle genome (panel B). Segmental duplications are enriched in breakpoints comparing to synteny blocks suggested that in these cases duplications promoted nonallelic homologous recombination (NAHR), and thus a chromosome rearrangement. elements including LINEs, LTRs, SINEs, simple repeats and segmental duplications (SDs), as well as genes are enriched in breakpoints compared to synteny blocks. Each breakpoint region aligns with the two homologous centromeres radiantly and is sprayed with SDs. More interestingly, the breakpoint of the goat genome ( Figure 7A) presents a more complex architecture than cattle genome one ( Figure 7B). We inferred that these two cattle centromeres underwent a much longer evolution of SD than that of the two goat centromeres ( Figure 7A and B). It suggests that these two goat centromeres of CHI9 and CHI14 are formed by a Robertsonian fusion and is evolutionarily more recent. Though neocentromere formation may be independent of sequence characteristic and comply with the rules of epigenetic, DNA rearrangement is one of the primary driving forces for neocentromere formation in chromosome evolution. Abundant SDs could hypothetically provide homologous sequences that are not in allelic positions, where nonallelic homologous recombination (NAHR) or crossing-over in meiosis result in regionspecific genomic rearrangement [24]. These duplications can promote NAHR, and thus a chromosome rearrangement. Among distinctive sequence features, we noted that a gene cluster containing ANKRD26 (ankyrin repeat domain 26) was highly duplicated in ruminant breakpoints (Figure 7 and Additional file 6).
The ancestral Bovidae diploid number of 30 pairs of chromosomes has been retained in all but one (gaur) of the previous studied species of Bovina. Robertsonian fusions or centric translocations (ROBs) involving the fusions of two nonhomologous acrocentric or telocentric chromosomes at their centromeres (producing a single biarmed product with distinguishable internal centromere) is frequent in ruminant chromosome evolution and particularly in Bovidae [25][26][27]. Several hypotheses have been advanced to explain the formation of ROBs in mammals usually invoking the organization of pericentromeric satellite DNA and nonallelic homologous recombination at meiosis [28]. Breakpoints are considered to occur in regions of SDs [29] and high concentrations of repetitive elements [30]. The sequence architectures of genomic breakpoints what we identified for ROBs are found to be complex with enriched genes and other active elements, which refine the evidences of these previous studies. The fact that SDs are enriched in breakpoints compared to synteny blocks and that breakpoint specific SD have a different evolutionary history in goat and cattle suggests that in these cases duplications promoted nonallelic homologous recombination (NAHR), and thus a chromosome rearrangement.

ANKRD26 gene cluster associated with genome rearrangements
Large-scale genome rearrangements caused by chromosome breaks underlie numerous inherited diseases and are associated with species evolution. Evolutionary breakpoints are nonrandomly distributed (reused) throughout mammalian genomes [31] and many breakpoints are coincident with ancient centromere activity as well as new centromere formation [32][33][34][35][36]. It has been hypothesized that centromeric illegitimate recombination between nonhomologous chromosomes led to chromosome fusions and synteny breakpoints during rapid speciation. This hypothesis initiated from the analysis of karyotypic evolution in mammals [37] and was supported by similar studies in marsupials [32] and in plants [38].
From goat-cattle the comparison, we found that ANKRD26 has expanded in ruminant breakpoints. We then explored the genomic distribution of members of ANKRD26 homologs to other mammalian genomes. The ortholog of ANKRD26 is a highly conserved protein gene in reptiles, birds and mammals. We found that ANKRD26 have a single copy in each reptile genome or each bird genome as far as we can find in their available genome sequences, whereas it was expanded in mammals. Phylogenetic analysis of ANKRD26 homologs revealed four groups of homologous genes including ancient ANKRD26like (1), POTE (2), ANKRD18/ANKRD20 (3), and ANKRD36/ANKRD62 (4) and an over-expansion especially in primate and artiodactyla (Additional file 7). The analysis of evolutionary breakpoints using the Evolution highway program [33,39,40] was used to determinate whether a breakpoint is species-specific or lineagespecific (i.e. reused). We combined of the map of mammalian breakpoints and the genomic distribution ofANKRD26 homologs (Figure 9 and Additional file 8). The result showed that most of ANKRD26 homologs locate on species-specific or lineage-specific rearrangement regions or centromeres, whereas a very few (6/115) were located on homologous synteny blocks (HSBs). For instance, the POTE family, made of primate-specific paralogs of ANKRD26, was extensively reshaped involving segmental loss and internal duplication in human [41]. Genes in the POTE family are located near the centromeres of numerous chromosomes including the site of human-chimpanzee chromosome 2 fusion [42]. Of the 22 human autosomes, 11 centromeres including all five acrocentric ones (that always resulted from Robinson translocation) have ANKR26 paralogs. The ANKRD26, as the ancestor of this big family, is located on cattle (or goat) chromosome 13 and human chromosome 10p, and its homologous copies are spotted intensively in two neighboring breakpoint regions. One likely reason Figure 9 An example of the map of mammalian breakpoints and ANKRD26 homologs. Cattle chromosomes are used as reference chromosomes and are respectively compared to goat RH maps, pig (S. scrofa), horse (E. caballus), dog (C. familiaris), human (H. sapiens), macaque (M. mulatta), mouse (M. musculus), and rat (R. norvegicus) in the Evolution Highyway program. Large blocks of homologous synteny and a high frequency of breakpoint reuse are presented on mammalian chromosomes. We combined of this map of mammalian breakpoints and the genomic distribution ofANKRD26 family (Additional file 8) and take examples of cattle chromosome 13 in this figure. Most members of ANKRD26 family locate on species-specific or lineage-specific rearrangement regions or centromeres (red cross), whereas a very few (blue cross) were located on HSBs. The homologous copies of ANKRD26 are spotted intensively in two neighboring breakpoint regions, where ANKRD26 locate. One likely reason for of this excessive gene duplication is that the HSB between the two neighboring breakpoints have undergone repeatedly inversion, thus ANKR26 copies accumulated.
for of this excessive duplication is that the HSB between the two neighboring breakpoint have undergone repeatedly inversions, thereby provoking an accumulation of ANKR26 (Figure 9 and Additional file 6). This means that species-specific shifts in ANKRD26 duplication are coincident with breakpoint reuse in divergent lineages.
Evolutionary breakpoints are supposed to contribute importantly to new genetic variation and novel genes at these fragile sites [39]. However, the dynamics between evolutionary breakpoints and genomic elements locally have not been completely understood yet [43]. A current hypothesis is that genetic innovation and change at these sites may, as a driving force of local adaptation [44], have influences in subsequent chromosome stabilization. Indeed, most of the proteins belonging to this gene cluster conserved the Ankyrin domains near to N-terminal and the structural maintenance of chromosomes protein (SMC) domains near the C-terminal (Additional file 9). SMC proteins interact with DNA in chromosome condensation, sister-chromatic cohesion, recombination, DNA doublestrand break repair and epigenetic silencing of gene expression. It may indicate that ANKRD26 is responsive to DNA damage [45], pericentromeric cohesin and condensing [46], and even centromere formation. The members of ANKRD26 usually linked with LINEs near the C-terminal [47] may indicate their duplication by retrotransposition.
Although little is known about the regulation of ANKRD26, previous studies demonstrated that ANKRD26involved in chromosome rearrangements are responsible for several inherited disease: a gene conversion of ANKRD26 is associated with type-2 diabetes [48], which were supported by the studies of ANKRD26-inactivation [49] and familial thrombocytopenia, associated with a chromosome breakage [50] on human chromosome 10p [51] appears to be due to an ANKRD26 mutation [52,53].

Conclusions
In the present study, we constructed dense goat RH maps of SNP markers from closely related species via a novel method of genotype calling and recent methodologies. Overall, there is good agreement between the current genome assembly and the robust RH maps. We were able to propose correction to the goat genome assembly and add previously unplaced scaffolds. We combined this information in the form of a new goat genome sequence CHI_1.1.We think this updated sequence, with a more accurate ordering of scaffolds will can be useful for many genetic studies on goat, for instance, linkage disequilibrium analysis and QTL mapping. We used our new dense chromosomal maps to annotate genomic rearrangements and chromosome evolution of ruminants. We compared sequence architecture of the evolutionary breakpoints including representative repeat elements and genes in ruminant rearrangements. Breakpoint intervals from two neighboring synteny blocks in goat gnome reveal a mosaic of new insertions at rearrangement sites. The Breakpoints in the goat genome presents more complex architectures than breakpoints of the cattle genome suggesting that goat centromeres of CHI9 and CHI14 underwent much shorter evolutionary time than those two cattle centromeres. Our findings highlight the complex interplay between duplicated sequences and chromosomal rearrangements, which favored a rearrangement mechanism of centromeric non-allelic homologous recombination in mammals. We detected a gene cluster, ANKRD26, which was expanded specifically in species as well as lineage specific evolutionary breakpoints and in centromeres. Their genomic distribution associated with the biological function of their conserved protein domains suggested it might play a role in chromosome stabilization during chromosome evolution.
We performed the genotype calling procedure described below independently for each SNP array. We called genotypes of clones at SNPs based on the raw signal intensities obtained from the Illumina genotyping system. Specifically, the statistic used for genotype calling of a given SNP in a given clone was the maximum observed intensity over the two possible alleles (herein called Imax). The distribution of Imax in a given clone can be seen as a mixture of two underlying distributions, where one for the non-retained SNPs in the clone and one for the retained SNPs (e.g. Figure 1B). The difficulty of the genotype calling procedure is to evaluate from which of these two distributions a given data point (Imax) comes from. More precisely, we aimed at evaluating the p-value of the observed Imax under the hypothesis H0 that the SNP is not retained in the clone (i.e. P (T > =Imax |H0)).
To calculate this p-value an estimate of the distribution of Imax for a non-retained SNP is needed. For this, we exploited the fact that the SNP array was designed for a different species than the goat (either cattle or sheep). Looking at the intensity distribution in the positive whole genome goat sample, we observed that a number (Nneg) of SNPs showed very low intensities for both alleles in the goat ( Figure 1A). We speculate that the sequences corresponding to these negative SNPs are not conserved enough in the goat genome to provide a signal with the SNP array. Finally, the p-value that a (non-negative) SNP is not retained was calculated as the proportion of Imax in the negative SNPs that exceeds the observed Imax. If none was found, we set the p-value at 1/Nneg. After calculating p-values for all observations, we estimated their corresponding q-values using the Q-value R package [18]. This approach allows controlling the false positive and false negative rates and provides an estimate of the proportion of true null hypotheses in the data, which is in our case an estimate of one minus the panel retention fraction. We called "retained" data points at an FDR threshold of 1% and removed SNPs which retention fraction was either greater than 50% or lower than 10%. Furthermore, one RH clone (#73) showed a very high rate of retention (>90%) and was therefore removed from the analysis. Among the remaining data points, we called "missing" data points with q-values lower than 10% for the sheep array, and 5% for the bovine array. Finally, data points with q-values greater than 10% were called "un-retained". The overall approach with these thresholds allowed us to control both the missing data rate and the false negative rate (see Results).

Map construction
Our pipeline of RH mapping is very similar to that recently used for recently building porcine RH maps [13]. The Carthagene software [12,14,54] was used to compute efficiently marker ordering for candidate maps with reference maps. This approach is based on a probabilistic Bayesian model integrating the usual RH probabilistic model with a probabilistic model of breakpoint occurrences with the reference order. In this probabilistic model, breakpoints induced by chromosomal rearrangements are considered as rare events, following a Poisson law. Meaning, the model considers that genome assembly errors create rare spurious breakpoints between the RH map order and the current goat genome assembly order [12].
The map construction process was conducted in parallel to the goat genome assembly. As our first goal to build RH maps was to help ordering scaffolds of the goat genome, only SNPs that could be mapped unambiguously to a working draft of the goat genome were used to build RH maps. These summed up to 54,318 SNPs, with 35,203 SNPs from the Ovine array and 19,115 SNPs from the bovine array. We portioned 53,075 mapped SNPs (mSNPs) according to their assigned chromosome and 1243 unassigned SNPs (uSNPs) remained on assembly CHIR_1.0. For each autosome, we determined linkage groups among mSNPs, based on RH data alone. Of the 53,075 SNPs that generated RH vectors, 53,062were assigned to 30 single linkage groups based on two-point analysis at a LOD threshold of 10. In these linkage groups 13 markers did not show enough evidence for linkage.
Each RH linkage group was assigned to a chromosome based on homologous comparison between RH markers and the goat genome. (1) We manually increased the LOD threshold in a two-point analysis to remove minor marker groups (<5 markers) to remove possible error vector. (2) A primary map was first generated with Carthagene by converting the RH data into a TSP and solved using a LKH heuristic method (lkh1 1). The initial RH maps were compared to the assembly to observe the general discrepancies, and were used to inspect inter-markers intervals with extremely long distances (>25 cR). In this step, we manually removed possible questionable markers, most of which were located on map ends. (3) The command of find_errors in RHMAPPER [19] was used to flags marker/hybrid assay results that are likely to be the result of laboratory error, where the order for these results set to the primary map. A threshold of 3 was specified in the command i.e. reporting that this result is 1000:1 more likely to represent a laboratory error than a true result. Once again, we converted the RH data into a TSP and solved using the LKH heuristic method, and constructed the LKH map.
A novel approach [14] was applied next to build maps with a highly reliable ordering, called robust maps. Briefly, the principle of the method is to estimate a posterior distribution of marker ordering using an MCMC approach and then to extract from this distribution a subset of markers that show the same ordering across maps of high probability. We performed 5000 MCMC iterations and discarded the first 1000 as burning iterations in the command of Carthagene software, mcmc, using the LKH maps obtained above as starting orderings. We extracted robust maps from the posterior distribution using the metamap software [14]. An inclusion tree called metamap that summarizes the uncertainty lying in the map distribution was calculated. In the metamap tree, groups in which the best order have posterior probability higher than 95% in the robust map were include the final robust map. The distances of the robust maps were evaluated using the diploid equal retention model with an EM convergence set to 10 −9 (command cgtolerance). This approach is closely related to the idea of constructing framework maps. In the case of framework maps, an order is accepted based on a maximum likelihood ratio, also called LOD (the logarithm of odds between the best order and the second best order must exceed a preset ratio, such as 1,000:1 for example (LOD 3)). In contrast, the construction of robust maps falls into the Bayesian paradigm [14].
For mapping the unplaced scaffolds on the goat genome, we constructed primary LKH maps with mSNP and uSNP together. We removed 852 uSNPs, which assigned to scaffolds shorter than 10 Kb in RH mapping from 1,243 unassigned SNPs (uSNPs) on assembly CHIR_1.0. Then 391 uSNPs were assigned to chromosomes based on their linkage with each of 30 linkage groups. The LKH maps were constructed as above. We aligned both of the robust maps of only mSNP and the LKH map of all SNPs to the goat scaffold data, then computed the orientation of scaffolds positions on the chromosomes and the orientation of the scaffold. A comparison of the two results showed that the assignments of placed scaffolds are very similar, so we proposed that the assignments of placed scaffolds using the LKH maps are accurate.

Genome map comparisons and gene comparisons
Map discrepancies were determined by comparing the RH map with target genome data including goat scaffold sequence, goat assemblies (CHIR_1.0 and CHIR_1.1) and cattle genomes (BTA4.2 and UMD3.1). The probe sequences of markers in RH maps aligned with these target genome maps using BLASTN (version 2.2.21) with E-values < 1 × 10 −10 . We also conducted synteny-based comparisons between the goat genome sequence and other mammalian sequences. The comparisons were performed using NUCmer [55] with parameter "-c 500" (with cattle) or "-c 150" (with other mammals). The assignment of assembled sequences to chromosomes was done based on a best alignment and majority rule. The Interquartile Range method was used to filter out suspicious markers and compute the scaffold positions on chromosomes. The Least Square method was used to compute the orientation of the scaffold that had multiple markers aligned to them. Orientation was considered forward if the slope value was positive, reverse if the slope value was negative and unknown if the slope was close to be 0.
For the reconstruction of the evolutionary breakpoints, the Evolution Highway program [33,39,40] was used to compare more mammalian genomes including human, mouse, dog, horse, cattle, rat and goat RH maps. The data of orthologous gene pairs between cattle genome and other mammalian genome for generating homologous synteny block was downloaded from the Biomart database. We used lastz program and lastz-based Cassis program [56] to precisely detect sequence homology and genomic rearrangement breakpoints. RepeatMasker was run against both libraries separately to identify homologous repeats, which were classified into known classes of repeats [57]. SDs was identified using the identical whole-genome sequence comparison method [58]. The characterized elements of sequence were graphed by circus [59].
For phylogenetic analysis of gene family, the sequence of goat ANKRD26 was blasted (BlastP) against the NCBI-nr database. In the results, genes with less than 40% homology or 40% coverage were filtered out. The annotated genes of ANKRD26 in public protein databases also were retrieved. Multiple sequence alignment was performed and the Neighbor-joining tree were constructed based on Dayhoff PAM matrix method with a bootstrap value of 1000 [60].