text:
Mechanical properties of DNA origami nanoassemblies are determined by Holliday junction mechanophores

DNA nanoassemblies have demonstrated wide applications in various fields including nanomaterials, drug delivery and biosensing. In DNA origami, single-stranded DNA template is shaped into desired nanostructure by DNA staples that form Holliday junctions with the template. Limited by current methodologies, however, mechanical properties of DNA origami structures have not been adequately characterized, which hinders further applications of these materials. Using laser tweezers, here, we have described two mechanical properties of DNA nanoassemblies represented by DNA nanotubes, DNA nanopyramids and DNA nanotiles. First, mechanical stability of DNA origami structures is determined by the effective density of Holliday junctions along a particular stress direction. Second, mechanical isomerization observed between two conformations of DNA nanotubes at 10–35 pN has been ascribed to the collective actions of individual Holliday junctions, which are only possible in DNA origami with rotational symmetric arrangements of Holliday junctions, such as those in DNA nanotubes. Our results indicate that Holliday junctions control mechanical behaviors of DNA nanoassemblies. Therefore, they can be considered as ‘mechanophores’ that sustain mechanical properties of origami nanoassemblies. The mechanical properties observed here provide insights for designing better DNA nanostructures. In addition, the unprecedented mechanical isomerization process brings new strategies for the development of nano-sensors and actuators.

INTRODUCTION
Because of the highly stable and specific recognition between two complementary DNA strands, DNA has been used as an attractive component in nanoassembly. In the DNA origami nanoassembly, a long ssDNA serves as a template to fold into nanostructures through hundreds of Holliday junctions formed between short DNA staples complementary to the template sequences at particular locations (1). The simple, robust and highly efficient synthesis strategy of DNA origami has established this structure as a highly potent nanomaterial (2)(3)(4)(5)(6)(7)(8)(9)(10) yet to be fully characterized for its properties. Among these unknown territories, mechanical property is certainly a notable missing link with high significance. Mechanical stability of the connecting regions in DNA nanoassemblies is essential to sustain robust interactions between biomolecules (11) and nanoassemblies. Likewise, mechanical rigidity of DNA nanocavities is critical to define the morphology of inorganic nanoparticles contained within (12,13). In the sensing applications using DNA nanoassemblies, (14)(15)(16)(17)(18), the mechanical rigidity and stability of the DNA nanoassemblies can directly affect the accuracy in the signal output and the sensitivity in the analyte recognition (5).
So far, only a handful of investigations have been reported with a main focus on the mechanical functionalities of the DNA nanoassemblies (14,(19)(20)(21)(22). The insufficient information on the mechanical rigidity and stability of DNA nanoassemblies in general and DNA origami in particular hinders rational design of DNA nanomaterials for applications that exploit their mechanical properties. One reason for this lack of knowledge lies in the difficulty in the char-acterization of individual nano-objects for their mechanical properties. Special tools with high resolution for mechanical force measurement, such as AFM and optical tweezers, must be employed to carry out the characterization (22,23) after immobilization of nanoparticles with different shapes, which is another challenging practice.
We reasoned that due to the flexible nature in the design of a DNA origami with the single-nucleotide precision, it is possible to introduce two duplex DNA handles to tether a DNA origami. These handles are then linked to the two optically trapped polystyrene beads in laser tweezers, which allow the quantification of mechanical properties of DNA nanoassembly. Previously, with the aim to develop origami based nanomechanical devices, Sugiyama and coworkers have found that tubular designs of origami structures can have two stable conformations: a short and a long tubular forms that may stem from different isomers of Holliday junctions contained inside the origami devices (10). Given the potential applications of controlling conformations of DNA origami structures by mechanical factors as discussed above, here we wish to understand whether it is a unique feature for the tubular-shaped origami to demonstrate different structural isomers. To this purpose, we applied two duplex DNA handles to the tubular DNA origami structures for the comparison of their mechanical properties with those from other nanoassemblies, DNA nanopyramids and DNA nanotiles. We found mechanical stabilities of DNA origami structures are correlated with the effective density of Holliday junctions in a particular nanoassembly. The mechanical stability is anisotropic in nature with the short axis of a DNA nanotube resisting higher external stress than the long axis. Interestingly, mechanical isomerization between two conformations of a DNA nanoassemby at 10-35 pN external force was observed only in DNA nanotubes, which have unique symmetric arrangements of Holliday junctions. Similar to the mechanical stability, the mechanical isomerization also showed anisotropic behavior. Given that individual Holliday junctions have isomerization force in the sub-picoNewton range (24) and they are anisotropically arranged in DNA nanotubes, we have attributed the mechanical isomerization of the DNA nanotubes to the collective actions of many Holliday junctions that experience similar microenvironment. All these results indicate that Holliday junctions in DNA origami structures serve as mechanophores (25,26), which determine the mechanical property of DNA nanoassemblies, similar to chromophores and fluorophores that carry spectroscopic information in a molecule. We anticipate these new findings are instrumental to optimize the mechanical strength of DNA nanoassemblies.

Mechanical isomerization and mechanical disassembly of DNA tubes
Optical tweezers have been used to apply and measure forces in picoNewtons (27)(28)(29). Compared to AFM, it has better force resolution and therefore is particularly suitable to characterize mechanical stability of macromolecules such as proteins and nucleic acids. Here, we used set up for the mechanical isomerization of DNA origami tubes. DNA origami tube is sandwiched between two long dsDNA handles and the whole construct is tethered between two optically trapped beads with affinity linkages. The mechanochemical property of DNA tubes is revealed by moving one of the beads away from another at a loading rate of ∼5.5 pN/s. AFM image of the eight-tube DNA for longitudinal (B) and horizontal stretching (C). Schematic pictures are shown to the right of corresponding AFM images. a home-built laser tweezers instrument (30,31) to characterize the mechanical properties of individual DNA origami nanoassemblies. To this end, we introduced two doublestranded (ds) DNA handles to DNA origami structures (see Figure 1A, Supplementary Figure S1A and B and Materials and Methods) (14). DNA handles were inserted into the origami structure via overhang single-stranded staple sequence, which contains 40 nucleotides complementary to the single-stranded M13mp18 template (it has been reported that the shearing force of 30 bp duplex DNA is >60 pN) (32). The other end of the handle was labeled with a digoxigenin or biotin molecule. These DNA handles were mixed and incubated with rest of the DNA components for the origami construction (see Materials and Methods and Supplementary Figures S1-S7). AFM images have revealed successful incorporation of the handles (see Figure 1B  Using the digoxigenin and the biotin molecules labeled at the free end of the DNA handles, the modified DNA origami construct was tethered to the two optically trapped particles coated with digoxigenin antibody and streptavidin through respective affinity interactions (14). The mechanical stability of the origami was probed by force ramping experiments in which one of the optically trapped particles was moved away from the other at a loading rate of ∼5.5 pN/s in laser tweezers.
Due to the unique behavior of isomerization displayed in DNA origami tubes (10), we chose these origami structures as our first samples to evaluate mechanical properties of DNA origami structures. Using the reported procedures (10), we prepared eight-tube and six-tube DNA origami nanoassemblies that respectively consist of eight and six Holliday junctions in each of the circular layers that are spirally arranged into desired tubes. Previously, AFM images (10) have revealed that each type of tube has two equilibrated isomers, long and short tubes (similar to those shown in Figure 1B and C, Supplementary Figures S8 and  S9). Due to the limited time resolution in AFM, the isomerization between the two conformers has not been observed in real time. The dynamic nature of mechanical unfolding at the single nanostructure level may provide an unprecedented opportunity to investigate this isomerization process. To this end, by using the procedures reported earlier (14), either eight-tube (short and long) or six-tube (short and long) DNA origami structures were obtained and tethered to two optically trapped polystyrene beads as described above.
In the first set of experiments, we mechanically stretched DNA origami samples from longitudinal direction of the tube (see the eight-tube DNA in Figure 1 as an example). We observed two types of rupture transitions ( Figure 2A and B), which were assigned to mechanical isomerization and disassembly of origami tubes on the basis of the following aspects. First, the two types of transitions can be separated according to the low (10-35 pN) and high (35-60 pN) force of transitions ( Figure 2C, two left panels). The hierarchical cluster analysis (HCA) (33), which is a nonsupervised method (34) to classify groups using the distance between two data points, confirmed the presence of these two populations ( Figure 2C, two right panels). Here, we used rupture force as a distance variable in the HCA calculation (see Materials and Methods for details). When we scrutinized the properties of the transition in these two groups, we found that the single or multiple transitions in the lower force region (10-35 pN) were reproducibly observed after incubating the construct for 90 s before each force-extension cycle (Supplementary Figure S10). On the other hand, we observed successive transition events in the higher force region (35-60 pN) that were not reversible after the same incubation time. These high-force transitions, all centered in a similar force region (35-60 pN) even for different DNA tubes (Supplementary Figures S11 and S12), suggesting the disassembly of the structures with similar mechanical stabilities. Since the force range we observed here fell into that for the slide-opening of DNA duplex (35), we surmised these transitions may be due to the forceinduced disassembly of Holliday junctions (shearing breakage), whose re-assembly is slow. Consistent with this scenario, when the tension for the origami in the higher-force region was relaxed to 0 pN, we did not observe the same trajectory in subsequent force-ramping cycles. Instead, the contour length of the origami became longer even after incubation up to 90 s at 0 pN (see curves in the disassembly region in Figure 2B). Similar events were observed during the mechanical stretching of DNA nanotiles, which have been ascribed to the disassembly of DNA origami structures (14). Therefore, these high-force transitions likely represent the disintegration of a DNA origami structure through the irreversible breakage in the Holliday junctions during the time scale of the experiment.
To explain the nature of the low-force transitions, we analyzed the population of the origami that yielded re-  Figure S13B) and 75% (Supplementary Figure S10) of the origami showed these transitions in the six-tube and eight-tube DNA stretched from the longitudinal direction, respectively. During longitudinal stretching, only short tubes can isomerize into the long tubes ( Figure 2D), therefore, 42% and 75% populations represent short tube isomers. These values are close to that obtained from AFM images, which showed 45% and 85% (Table 1) for the short isomers in the six-tube and eight-tube DNA, respectively. The striking agreement between the mechanical unfolding and AFM analyses strongly suggests that these low-force transitions are likely force-induced isomerizations, which are expected to be reversible after incubation as the origami structures remain intact during the process ( Figure 2D and E). To confirm this scenario, we measured the changes-in-extension due to the rupture transitions in the lower force region. We found that the measured values were closely matched with the calculated difference in length between the short and long isomers in each type of the DNA tube (see next section for details). It is noteworthy that in the force range (10-35 pN) we observed for the isomerization, the changes-in-extension for different DNA nanotubes do not vary with force significantly (calculations see SI).
Based on these facts, we assigned the lower-force transition as mechanical isomerization while the higher-force transition as disintegration of origami structures (see schematic models in Figure 2D (nanostructure level) and Figure 2E (molecular level)).
Given that each Holliday junction consists of Watson-Crick base pairing that can be taken apart by different mechanical force from varying directions, (24,36) it is interesting to see whether or not the same force is required to dis- A set of F-X traces for another eight-tube DNA origami. The first transition event in the red curve depicts mechanical isomerization. The relaxing curve after mechanical isomerization is shown in black. After incubation, the next stretching curve is identical with the first red curve. After mechanical disassembly, the traces do not overlap with the first red curve after incubation at 0 pN. Instead, they depict longer contour lengths (blue). (C) Distribution plots for the transition forces (left two panels) and the hierarchical cluster analyses (right two panels) clearly indicate the presence of two populations. The red population with lower transition force is due to isomerization of the nanotube while the blue population with higher transition force is due to nanotube disassembly. Solid curves in the left panel represent Gaussian fittings. The green data points are identified as stand-alone groups in HCA. We assigned these points to each of the populations based on force. (D) Schematic of mechanical isomerization of the short to long eight-tube DNA origami structures followed by disassembly. Magnified sections show the rearrangement of Holliday junctions before and after mechanical isomerization. (E) Schematic of the molecular rearrangements of Holliday junctions in a DNA tube during mechanical isomerization and disassembly processes. Numbers in black and green respectively indicate the duplex layers and the Holliday junctions in DNA origami tubes. integrate the DNA origami tube from different directions. To test this scenario, we stretched the same type of tube from horizontal direction by using two dsDNA handles facing each other in the middle of the eight-tube DNA (Figure 1C and Supplementary Figure S3). Similar to what was observed for the longitudinal stretching, two regions were found for the horizontal stretching. However, the disassembly process was abrupt. While the lower force region (10-35 pN) showed reproducible features after incubation at 0 pN, the higher force region (35-60 pN) again revealed irreversible disassembly, likely due to disintegration of staples. Similar characteristics were observed during the longitudinal and horizontal stretchings of the six-tube DNA (see below for detailed discussion).

Mechanical disassembly (or stability) of a DNA origami structure is correlated with the effective density of Holliday junctions
When we scrutinized the higher-force mechanical disassembly region that reflects the mechanical stability of a DNA nanoassembly, we found that the mechanical stability for horizontal stretching of the eight-tube DNA was higher (49 pN, see Supplementary Figure S13A for F-X curves and Figure 3B for histogram) than longitudinal stretching (42 pN, Figure 2A and B for F-X curves and Figure 3A for histogram). Similar anisotropy in mechanical stability was observed for the six-tube DNA origami (see Supplementary  Figures S4 and S5 for its structure and preparation and Supplementary Figure S9 for AFM images). While the horizontal mechanical stability was 49 pN, that of longitudinal was lower (44 pN, summarized in Table 1 and Figure 3C, see representative F-X curves in Supplementary Figure S13 B&C). During disassembly, gradual deviation in the extension that corresponds to the sequential rupture of Holliday junctions was not observed in the F-X curves for either horizontal or longitudinal pulling direction. Instead, sudden and big (80-100 nm, see Figure 2A and B) rupture events that were consistent with the layer-by-layer disintegration during longitudinal pulling were observed ( Figure 3E). Peeling off each layer during layer-by-layer disassociation in the eight-tube DNA requires the breakage of eight Holliday junctions (see black arrows in Supplementary Figure S2), while it requires breaking six Holliday junctions for the same effect in the sixtube DNA (see black arrows in Supplementary Figure S4). On the other hand, with horizontal stretching, the layer-bylayer peeling is rather difficult (Supplementary Figures S3  and S5) and the disassembly of the structure is likely due to force induced melting (35,37) of Holliday junctions (Figure 3F). As a result, complete disassembly of many more Holliday junctions is necessary to break the tube. During stretching experiments, it is more likely that the loss of the Holliday junction connected to the DNA handles would cause the sudden breakage of the tether (Supplementary Figure S13A and C). These breakage forces are used to represent the lower end of the disassembly forces of the origami tubes at the horizontal direction ( Figure 3C). This geometry difference also predicts that during longitudinal stretching, more intermediates can be accommodated compared to those from horizontal direction. In fact, this has been confirmed during experiments, in which on average 6 inter- Histograms of the mechanical disassembly forces measured along longitudinal and horizontal stretching directions of DNA eight-tubes, respectively. The disassembly force was determined at the first disassembly event (the red arrow in Figure 2B). The insets depict the stretching directions.  Figure 3D).
To further characterize the mechanical stability of DNA origami structures, we performed stretching experiments of origami nanotiles (14) (see Supplementary Figure S6 for design and AFM image, S13D for F-X curves) and pyramidshaped origami (see Materials and Methods and Supplementary Figure S7 for the design and AFM image, Supplementary Table S1 for the staple sequences and Supplementary Figure S13 E for F-X curves). As shown in Figure 3C, DNA nanotiles have shown the weakest mechanical stability (disassembly force 28 pN), which was followed by a 3D structure, DNA nanopyramid (35 pN). The most compact 3D structures, horizontal eight-tube and six-tube DNA, demonstrated the strongest stabilities (49.4 and 48.5 pN respectively). Due to the open geometry in the 2D origami nanotiles (14), disassembly of the structure via stretching only involves the breakage of Holliday junctions located close to the axis of the stress. In contrast, in the 3D DNA assemblies, the stability of a structure can be contributed from distally located Holliday junctions that come close to Nucleic Acids Research, 2016, Vol. 44, No. 14 6579 the stress axis through the long-range, 3D arrangement. In fact, the density of Holliday junctions along the longitudinal pulling for the 3D tubular structure increases significantly as the Holliday junctions in the same layer becomes more compact after isomerization (see Supplementary Figure S15). Therefore, the effective density of the Holliday junctions along the stretching direction becomes a decisive factor for the mechanical stability of an origami structure along that direction. Here, we defined the density of Holliday junctions as the number of contributing Holliday junctions per nanometer of the force axis in the nanostructure. For nanotiles and nanopyramids, the effective Holliday junction densities are 0.17 and 0.18 HJ/nm, respectively. For the eight-tube and the six-tube DNA along the longitudinal pulling direction, due to the symmetric geometries, all the Holliday junctions contribute equally to the mechanical property of the origami. Therefore, they are all counted in the density calculation, which leads to much higher densities (1.43 HJ/nm for the eight-tube and 1.09 HJ/nm for the six-tube). Corresponding to this trend, the disassembly forces of the isomerizable structures along the longitudinal direction for both six-and eight-tube DNA origami are higher (35-60 pN) than those of the non-isomerizable structures such as nanotiles and nanopyramids (28-40 pN). This clearly shows that the mechanical stability of DNA origami structures may depend on the effective density of the Holliday junctions. Therefore, the observed data of mechanical stability suggest that Holliday junctions serve as mechanophores (26,38,39) to impart the mechanical stability to the entire DNA origami nanoassemblies.

Holliday junctions in DNA origami tubes are responsible for mechanical isomerization between nanotube isomers
Next, we investigated the mechanical property of the DNA origami at the lower force region (Figures 2A-C and 4), in which transition events were reproducibly observed during consecutive stretching of the DNA nanotubes (see Supplementary Figure S10). F-X curves obtained during the stretching of the 8-tube DNA showed characteristic changes-in-extension along the longitudinal ( Figure 4C) or horizontal direction ( Figure 4D). This confirms that the low-force transition events represent isomerization of the DNA tubes. To determine the nature of these rupture transitions, we measured the change-in-extension during each event of isomerization transition. The change-in-extension is described as the isomerization length (L iso ) in Figures 4  and 5. For the eight-tube DNA, we found two major transitions of ∼62 and ∼26 nm during the longitudinal stretching ( Figure 4E) and one transition of ∼11 nm during the horizontal stretching ( Figure 4F). These change-in-extension histograms show broad distributions, which is expected as the transition events are likely mediated by one or more intermediates presented in DNA origami structures with large sizes (40,41). It is also possible that some reversible disassembly or melting processes may contribute to this broad feature. For example, it has been observed that under mechanical stress, force induced melting can weaken Watson-Crick base pairing (32), which may lead to transitions with longer-than-expected sizes. To clearly determine the most likely change-in-extension values, we used Population De- Dotted yellow lines depict Holliday junction layers in different stretching orientations. Representative F-X curves were obtained during longitudinal (C) or horizontal (D) stretching of the eight-tube DNA. Arrow heads in both F-X curves indicate the isomerization events (see insets for respective schematic diagrams). The blue/red and green traces depict stretching and relaxing curves, respectively. Due to the nature of the mechanical isomerization, only the transitions from the short tube to the long tube in the longitudinal pulling, and the long tube to the short tube in the horizontal pulling can be observed. Histograms of the isomerization length (L iso ) of the eight-tube DNA are plotted for longitudinal (E) or horizontal (F) stretching (inset in E depicts the most probable populations using PoD-Nano approach, see text). Histograms of the isomerization force (F iso ) of the same DNA tube measured during longitudinal (G) or horizontal (H) stretching. Solid curves represent Gaussian fittings. N and n respectively depict the number of F-X traces and number of molecules in experiments. convolution at Nanometer resolution (PoDNano), a statistical tool based on bootstrapping analysis as described in literature (42,43). With a total of 5000 random resampling processes, we confirmed the presence of the ∼62 and ∼26 nm transitions in the change-in-extension histogram with the broadest distribution ( Figure 4E inset). In AFM imaging, we have observed two different DNA tube isomers, long and short tubular structures (Supplementary Figure  S8) (10). Compared to the expected change-in-extension between these two isomers (63 and 14 nm for longitudinal and horizontal stretching directions, respectively, see Supplementary Figure S16 for detailed calculation), the most probable values obtained from the histograms of the change-in-extension ( Table 1) strongly suggest that the transitions are due to the force induced isomerization between the long and short DNA eight-tubes ( Figure 4A and B). Stretching of the six-tube DNA origami confirmed that the transitions in the low force region are indeed due to the mechanical isomerization between the long and short DNA tubes (80 and 9.2 nm in sum of all transition events vs expected values of 83 and 12 nm for the longitudinal and horizontal stretching directions, respectively, see Figure 5, Supplementary Figures S16 and S17). It is noteworthy that multiple transitions were observed in every longitudinal stretching experiment. This suggests the presence of intermediates during isomerization. For example, a stable intermediate at ∼46 nm ( Figure 5E) can be attributed to a conformation trapped during incomplete isomerization (Supplementary Figure S18).
Interestingly, during horizontal stretching, the mechanical isomerization involves fewer intermediate states compared to the longitudinal stretching (compare between Figure 4E and F; and between Figure 5E and F). This reflects different origami arrangements along various stretching orientations. During longitudinal stretching, it travels longer (62 nm for the eight-tube and 80 nm for the six-tube) from the shorter DNA tube to the longer tube with respect to the horizontal stretching, in which the distance is much shorter between the two isomers (11 nm for the eight-tube and 9.2 nm for the six-tube). The extra transition distance in the former case (longitudinal) renders a more elastic energetic profile that allows more intermediates to populate. In the latter case (horizontal), however, the entire process becomes more cooperative due to the shorter and therefore more rigid isomerization pathways that cannot host as many intermediate populations. Such an observation is consistent with the disassembly of the origami structures at the high force region discussed above. While the DNA tube experiences more intermediates during the longitudinal, layer-by-layer disintegration, it has fewer intermediates during the horizontal disassembly with increased cooperativity ( Figure 3D). Similar observations have been previously obtained in the mechanical unfolding of short and long DNA duplexes. While the short DNA duplexes such as hairpins almost do not present intermediates, (44) the sawteeth features observed in force-extension traces demonstrate the existence of many intermediates in long DNA duplexes (41).
During the longitudinal isomerization of the eight-tube DNA origami, the isomerization force (F iso = 26.4 pN, Figure 4G and Table 1) is comparable to that of the horizontal F iso (27.8 pN, Figure 4H and Table 1). Similar forces between the two stretching orientations have been observed in the isomerization of the six-tube DNA origami (F iso = 27.3 pN for longitudinal and 29.7 pN for horizontal, see Figure  5G and H, respectively). We reasoned that the isomerization of the DNA nanotubes can be a result of collective isomerization of individual Holliday junctions at the microscopic level. Based on this assumption, we calculated the isomerization force for individual Holliday junctions as 0.12-0.14 pN (see SI). This value is in agreement with the isomerization force of a single Holliday junction (0.1-0.3 pN) predicted by Hohng et al. (24).
It is noteworthy that mechanical isomerization was not observed in DNA nanotiles or nanopyramid structures. Close inspection on all three types of origami nanoassemblies shows that the latter two structures are rather different from the tubular structures. While the DNA nanotiles have a planar geometry, nanopyramids consisting of four triangular tiles fold into an object in which the perimeter of each layer gradually reduces from the base to the apex of the pyramid. In DNA nanotubes, the circular geometry with the same tube diameter gives rise to less hierarchical contributions between the distal and proximal regions with respect to the force axis, which bring a similar microenvironment to all Holliday junctions. However, due to the loss of symmetry in either DNA nanotiles or nanopyramids, individual Holliday junctions experience different environment and therefore, behave differently. As a result, the mechanical isomerization between two origami isomers can only be carried out by many identical Holliday junctions as a collective action in DNA nanotubes, but not in other structures.

CONCLUSIONS
Serving as essential components inside DNA origami nanoassemblies, DNA staples offer stability to DNA nanoassemblies by forming Holliday junctions with the single-stranded DNA template. Using optical tweezers, we have quantified and explained two emergent properties of DNA nanoassemblies not seen in their individual components from mechanical perspective. First, we have found the mechanical stability of DNA origami structures can be geometry dependent, which is determined by the effective density of Holliday junctions along a particular stretching direction. Second, we have quantitatively ascribed the cooperative transition between short and long DNA nanotubes to the collective mechanical isomerizations of many individual Holliday junctions. These observations indicate that the Holliday junctions serve as mechanophores in DNA origami nanoassemblies investigated here. Further testing is required to validate this point for more complex DNA nanoassemblies such as DNA helix bundles (2,21,22,45). We anticipate our new findings can provide unprecedented guidelines to design DNA nanoassemblies with better mechanical properties in both thermodynamic and kinetic aspects.

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.


extfieldsofstudy: Biology Medicine

text:
The pig X and Y Chromosomes: structure, sequence, and evolution

We have generated an improved assembly and gene annotation of the pig X Chromosome, and a first draft assembly of the pig Y Chromosome, by sequencing BAC and fosmid clones from Duroc animals and incorporating information from optical mapping and fiber-FISH. The X Chromosome carries 1033 annotated genes, 690 of which are protein coding. Gene order closely matches that found in primates (including humans) and carnivores (including cats and dogs), which is inferred to be ancestral. Nevertheless, several protein-coding genes present on the human X Chromosome were absent from the pig, and 38 pig-specific X-chromosomal genes were annotated, 22 of which were olfactory receptors. The pig Y-specific Chromosome sequence generated here comprises 30 megabases (Mb). A 15-Mb subset of this sequence was assembled, revealing two clusters of male-specific low copy number genes, separated by an ampliconic region including the HSFY gene family, which together make up most of the short arm. Both clusters contain palindromes with high sequence identity, presumably maintained by gene conversion. Many of the ancestral X-related genes previously reported in at least one mammalian Y Chromosome are represented either as active genes or partial sequences. This sequencing project has allowed us to identify genes—both single copy and amplified—on the pig Y Chromosome, to compare the pig X and Y Chromosomes for homologous sequences, and thereby to reveal mechanisms underlying pig X and Y Chromosome evolution.

Introduction
The therian (marsupial and placental mammal) sex chromosomes evolved originally from a homologous pair of autosomes about 170-180 million years ago (Cortez et al. 2014;Livernois et al. 2011), and have become extensively differentiated in terms of structure and sequence content. The gene content and organisation of the emergent X chromosome has been subject to strong conservation across different mammalian species with retention of much of the ancestral X (Ross et al. 2005;Bellott and Page 2010). In contrast, the acquisition of a dominant male sexdetermining function and accumulation of male benefit genes (e.g. genes involved in regulating male germ cell differentiation) on the Y chromosome has been accompanied by (a) the genetic isolation of much of the Y (suppression of recombination with the emergent X), (b) subsequent degeneration and loss of much of the ancestral Y gene content and (c) selection for dosage compensation in XX females to restore equivalence of gene expression between males and females for loci carried on the X that have degenerated or do not have a homologue on the Y (Bachtrog 2013;Graves 2010). Selection has also acted to retain a strictly X-Y homologous pseudoautosomal region (PAR) that permits X-Y pairing during meiosis and within which there is obligate recombination between the sex chromosomes. The gene and sequence content of the PAR varies between mammalian species and this reflects processes of expansion and pruning of the PAR in different mammalian lineages (Otto et al. 2011).
Comparisons of X chromosome sequences from several mammalian species have confirmed strong conservation of gene sequence and order (Chinwalla et al. 2002;Sandstedt and Tucker 2004). Groenen et al. (2012) published the first assembly of the porcine X chromosome as part of the initial description of the pig genome sequence, and again this demonstrated conservation of synteny across the chromosome. Nonetheless, sequence gaps and ambiguities remained within this first assembly, complicating genomic studies within pigs, and comparative studies between mammals.
In contrast to the broadly conserved X chromosomes, the hemizygous nature of the Y chromosome and suppression of recombination, in combination with normal processes of genome evolution, have led to a gradual degeneration of the chromosome over time, a large number of rearrangements, and colonisation by sequences from the X and autosomal chromosomes. Newly introduced genes will either drift or degenerate, or selection may act on variants to fix new genetic functions on the Y, particularly where these confer a benefit to the male. The haploid state of the sex chromosomes in males has generally led to the accumulation of male gametogenesis genes on both X and Y chromosomes (Vallender and Lahn 2004). A further consequence of the non-recombining status of the differential region of the Y is the relaxation of restraint on sequence amplification, leading to generation of ampliconic regions containing amplified gene and sequence families (Bellott and Page 2010).
The highly repetitive nature of many regions of mammalian Y chromosomes has impeded the generation of complete chromosome sequences; while there are tens of mammalian genomes sequenced, only a small fraction have a Y assembly. These few assemblies, plus several partial sequence assemblies have permitted the elucidation of chromosome topology and gene order in human and chimpanzee (Hughes et al. 2010), mouse (Soh et al. 2014), cattle (Elsik et al. 2009), horse (Paria et al. 2011), cat and dog (Li et al. 2013a). These works clearly show the great divergence in gene content, order, structure and sequence between Y chromosomes from different mammalian species. However, little data is available on the porcine Y chromosome sequence, gene order and their relationship to the X, despite the recent sequencing project for the pig genome (Groenen et al. 2012). Much of our knowledge of Y gene order comes from Quilter et al. (2002) who combined radiation hybrid mapping data with physical mapping of BAC clones to generate an ordered gene list. This paper presents a second-generation, much improved assembly and gene annotation of the porcine X chromosome. We present also a body of Y chromosome sequence, which has permitted a first-generation assessment of the Y chromosome short arm gene content and order, how this compares to other mammalian Y chromosomes, the evolutionary processes leading to the current Y organisation, and the structural relationships between the porcine sex chromosomes.

Results
A second-generation porcine X chromosome assembly

Organisation and evolution of the X chromosome
The pseudoautosomal region (PAR) in pig is of a similar size to the PAR in other closely-related mammals (e.g. cattle). We have previously discussed its gene content, its comparison to other mammalian species and the delineation of the boundary region (Skinner et al. 2013). Recently the precise location of the PAR boundary was confirmed to be within the gene SHROOM2 (Das et al. 2013). The X chromosome assembly (http://vega.sanger.ac.uk/Sus_scrofa/Location/Chromosome?r=X-WTSI) comprises 129,927,919bp of sequence in 5 contigs, with 13 gaps and an N50 length of 4,824,757bp.
Compared with the previous 10.2 build, many gaps have been filled and the order of sequences on the chromosome has been updated. Much of this improvement was aided by the use of optical mapping techniques, which helped resolve some of the more repetitive regions of the chromosome such as around the centromere; an example can be seen in the short arm clone CH242-202P13 (see Supplementary Methods for details on how the optical mapping approach was used here). Supplementary figure S1 shows a dot-plot alignment of the 10.2 X with our X assembly, highlighting the regions of the chromosome for which the sequence order has been corrected.
We aligned the X (and available Y) chromosome sequences of nine mammalian species ( Figure   1). The previously documented high level of conservation of X chromosome synteny is more apparent with the new pig X, as many of the reported breakpoints from cross-species comparisons to the build 10.2 X were due to errors that have been resolved in the new assembly. Specific lineages though, such as the rodents, have higher levels of chromosomal rearrangement. Li et al. (2013b) produced a genome assembly of a female Tibetan wild boar, and reported regions of the genome with apparent inversions with respect to the Duroc assembly. We compared the X inversions with our new assembly, and found that they lie outside the regions that have changed orientation from the 10.2 X. That is, these remain as potential inversions between Duroc and Tibetan wild boar.

Gene content of the X
Since the reference assembly of the X chromosome was high quality, the sequence underwent manual annotation, since the chromosome contains complex duplicated gene families such as olfactory receptors along with pseudogenes, which are hard to discriminate using only automatic pipelines. Table 1 shows the updated annotation compared to the 10.2 X build, plus comparable statistics for the annotation of the Y chromosome. The full gene annotation is provided in Supplementary Table S7 and is available through Vega browser. The majority (76%) of annotated loci in pig are shared with human. Again, the improved sequence assembly facilitated by optical mapping analysis assisted the new gene annotation, and many genes from the previous build that were disrupted by gaps have now been completed. The number of long noncoding RNA loci that have been identified on the new assembly has increased from 33 to 84, and, although the functionality of this category of loci is still under debate (Young and Ponting 2013), more and more evidence is coming to light suggesting that at least some of these loci have a functional role within the genome; for example, the lncRNA FIRRE has recently been identified in mouse and human and implicated in having a role in interfacing and modulating nuclear architecture across chromosomes (Hacisuleyman et al. 2014). Despite very low levels of sequence conservation between species, we find evidence for a FIRRE-like locus in the syntenic region of pig (OTTSUSG00000005757). Some genes with updated annotation from build 10.2 stand out as being of particular biological interest. Comparing the human and pig X chromosomes there are 38 coding loci in pig that are not found in human. Twenty-two of these coding genes, plus five novel pseudogenes, are in olfactory gene clusters. Pigs are known to have a large OR repertoire (Nguyen et al. 2012), and this adds to the reported collection. Figure 2 shows the improved assembly and annotation around one of the olfactory region clusters on Xq. The region lay within an inversion in the 10.2 assembly, corrected here and matching the gene order on the human X. The full list of genes present on the pig X chromosome, but not on the human X is provided in supplementary table

S9.
Besides the additional genes compared to human, there are also 11 protein coding genes which are present on the human X, which have been annotated as unitary pseudogenes (also known as Loss Of Function genes) on the pig X (Supplementary Table S10). These include GUCY2F (OTTSUSG00000005153) which in humans has been suggested as a possible candidate for involvement in X-linked retinitis pigmentosa (RP) (Yang et al. 1996). Pig models are now being developed for studying RP, taking advantage of the similarities between human and porcine retinal development (e.g. Ross et al. 2012;Fernandez de Castro et al. 2014). These models benefit from an improved understanding of the status of orthologues of disease-related genes: AWAT1 (OTTSUSG00000002936) is another unitary pseudogene. This gene is found on the X chromosome in a wide range of mammals, including cattle, humans and opossums (Holmes 2010). It encodes an acyl CoA wax alcohol acyltransferase involved in sebum production, and in humans is expressed in sebocytes, aiding in the prevention of surface desiccation of the skin (Turkish et al. 2005). The gene is a member of a family including DGAT1 and DGAT2, which have previously been associated with backfat thickness and intramuscular fat content respectively in pig breeds (Cui et al. 2011). It is likely that this pseudogenised AWAT1 in Duroc pigs represents a species-wide trait, as no transcripts for AWAT1 appear in the EST databases from any pig breeds.
Other identified unitary pseudogenes include ITIH6, a trypsin inhibitor, and RAB4, a member of the RAS oncogene family, as well as a number of zinc finger proteins and a transmembrane protein. Whether some of these loci were ever functional in pig, or merely reflect conserved regions that became functional at some point in the human lineage, is open for debate.
Other regions of difference lie in the cancer/testis (CT) antigen clusters found in humans and other primates but which are significantly reduced n pig. This is in line with evidence that enlarged primate CT antigen clusters arose due to a recent amplification in primates (Zhang and Su 2014), perhaps driven by a retrotransposition event. Their potential functions remain unknown, though they may have been involved in primate speciation.
Lastly it was found that INE1 (Inactivation-escape-1) is not present in pig. This is a non-coding transcript with unknown function within an intron of UBE1X (Thiselton et al. 2002). It appears to be unique to humans, where it escapes X inactivation in females.

A first generation porcine Y chromosome sequence assembly
Successfully assembled repetitive portions of Y chromosomes have been generated only for a limited number of species (see for example the human, mouse or chimpanzee Y; Hughes et al. 2010;Soh et al. 2014). Given the highly repetitive nature of the long arm of the pig Y chromosome (Quilter et al. 2002;Skinner et al. 2013), we targeted the short arm, which contains most, if not all, of the single-copy material.
The strategy employed was to prepare fosmid clones from flow-sorted Y chromosomes (origin Duroc) and initially to use fingerprint analysis to establish a framework of clones and to identify those inserts belonging to highly repetitive regions. Previously identified BAC clones (origin Meishan, from a Large White female/Meishan male cross: Anderson et al. 2000) containing Y chromosome genes conserved between several mammalian species served as anchor points for the elaboration of more complete sequence contigs. All of these BAC clones are known to map to the short arm of the Y (Quilter et al. 2002). Overlaps between BAC and fosmid clones were identified to extend the contigs, and selected clones within each contig were used for physical mapping (by metaphase-, interphase-and fibre-FISH). This, combined with existing data from linkage and radiation hybrid mapping, allowed us to produce an ordered map of the genes and contigs on the short arm of the Y. Selected repeat-containing clones were also sequenced to sample the diversity of repetitive sequence content of the Y chromosome. The final contigs we are predominantly of Duroc origin, with small regions of Meishan where no Duroc sequence was available. Supplementary table S1 summarises the ordered clones that we were able to incorporate into contigs including small contigs that could not be anchored to the physical map. The resulting male-specific Y sequence data is available at http://vega.sanger.ac.uk/Sus_scrofa/Location/Chromosome?r=Y-WTSI A two-pronged strategy was used to identify potential coding sequences in the body of fosmid sequence data. First, all fosmid sequence was subjected to gene prediction algorithms and potential coding sequences compared to the non-redundant protein database. Second, extensive BLAST analysis of known mammalian Y-linked genes and porcine X-linked homologues against the Y fosmid sequence dataset was used to identify porcine homologous sequences. The established intron/exon structure of these genes was used to drive the assembly of more extensive coverage of the porcine Y homologue. These sequences were used to identify ESTs in order to confirm expression status. Select genes were further tested by RT-PCR in a range of tissues. Figure 3 summarises the genes identified on the Y and their order, and the distribution of families of mammalian repetitive sequences across the short arm. The relative order and orientation of genes on the differential region of the Y short arm was established using a combination of the assembled contig sequences and high-resolution FISH with DNA fibres prepared by the Molecular Combing approach. The fibre-FISH also allowed an estimation of the gap sizes between adjacent contigs (see supplementary figure S2 for detailed presentation of supporting evidence). The genes are organised into two main blocks of low copy number sequences. These two blocks are separated by a region containing highly amplified sequences including the HSFY gene family (see Skinner et al, in submission). This amplified block maps to cytogenetic band Yp1.2 and is estimated from FISH analysis to be about 5Mb in length. Our final mapped sequence comprises seven contigs in the distal block and seven contigs in the proximal block. The final contig sequences were analysed for repeat content with RepeatMasker, and the densities of repetitive content found are shown in the column to the left of Figure 3. Two large (~400kb) contigs were assembled based on sequence overlaps and confirmed by fibre-FISH, but could not be assigned to the physical map (supplementary figure   S3). These unplaced contigs did not contain any identifiable genes, and are included in a separate assembly in VEGA.

Organisation of the porcine Y chromosome
Very few low-copy loci have been identified on the Y long arm. Indeed, FISH with clones containing male-specific repeats paints the entire long arm, indicating that it is mostly composed of repetitive sequences (e.g. Quilter et al. 2002). Although we attempted to sequence one of these clones, it was not possible to assemble a framework physical or sequence map from the repeats it contained, as most sequences collapsed into a single contig. The sequences we obtained from this clone are shown in the supplementary data and belong to previously published male-specific pig repeat families (Mileham et al. 1988;McGraw et al. 1988;Thomsen et al. 1992;Pérez-Pérez and Barragán 1998). Metaphase FISH did, however, reveal a small low-copy region at the q terminus (see supplementary figure S13).

Gene-related content of the Y chromosome
The Y chromosome sequence was run through the Otterlace/Zmap analysis pipeline which performs homology searches, de novo sequence analysis and gene predictions (Loveland et al. 2012). Repeat masking proved challenging due to the paucity of known pig-specific repeats.
Manual annotation resolved this as homology analysis is routinely run on-the-fly within the annotation tools, without repeat masking, to more accurately elucidate gene structures, especially using known Y chromosomal genes from other species as targets, and identifying the pig homologues where present. Many of the ancestral-X-related genes previously reported in other mammalian Y chromosomes are represented either as active genes or partial sequences, some of the latter with supporting ESTs. There is also evidence for rearrangement of a number of Y-linked genes which may have rendered them non-functional or modified their transcription products. RT-PCR analysis for several of these loci across a range of porcine tissues was performed to establish their transcriptional status. A list of the gene loci we identified in this and previous work is given in table 1. A fuller table including genes not tested here but with EST support is provided as supplementary table S3.

Ampliconic gene sequences
Although our sequence contigs are limited to the low copy regions of the chromosome, the FISH data show regions where gene-containing clones are present in multiple copies. Examples are the CUL4B genes, one of which is present as a partial copy in fosmid WTSI_1061-13A2. Probes against this sequence bind multiple targets proximal to SRY (and likely proximal to RBMY) as well as additional sites towards the centromeric end of the Y short arm (see supplementary figure S13). This indicates that CUL4B is likely to be an ampliconic gene. The sequence supports gene expression from a 'full-length' version of the sequence centromeric to SRY, and RT-PCR shows expression in testis, kidney and brain. A similar situation exists for clones containing TSPY.

Regions of X-Y Homology
We examined all the sequenced Y clones for homology with the X. Outside the PAR, two large (~50kb) regions of XY homology were seen, the remaining shorter regions comprising individual XY homologous genes and various classes of repetitive sequences. The overview X-Y comparison is given in supplementary figure S4). Most gene-related sequences show homology to the X.
We examined rates of synonymous substitution in open reading frames between the X and Y copies of genes conserved on Y chromosomes across mammalian species but within different evolutionary strata as defined for the human X and Y (Lahn and Page 1999) (Table 3). Gene pairs have similar substitution rates to their orthologues in other mammals.
A ~50kb region on the distal block encompassing the genes or gene fragments TRAPPC2-OFD1Y-GPM6B has high sequence identity to the X even in intronic and intergenic regions, suggesting a possible X-Y transposition (supplementary figure S6). OFD1Y is expressed highly in testis with lower expression in kidney and brain (supplementary figure S8), but there is no evidence for expression of the Y copies of TRAPPC2 or GPM6B. Furthermore, GPM6B is missing exons 1-3. GPM6B exon 1 lies near SMCY, possibly indicating an ancient inversion or other rearrangement. Recurrent X to Y transpositions of this region have been suggested before in cattle and dog (Li et al. 2013a;Chang et al. 2013); a gene tree constructed from dog, cattle, pig and horse OFD1X and Y sequences shows that the pig Y copy is more similar to the OFD1X than to any other species Y copies (supplementary figure S7), supporting the idea that this region arose on pig Y via a transposition from the X. A further ~33kb region of homology within the PAR (WTSI-X:1,840,693-1,874,130) can also be found on Xq 853,886,764), and potentially represents a region of duplication and transposition from the PAR onto Xq.

Inverted and duplicated blocks of sequence
Inversion and duplication to form palindromes has occurred around both SRY and CUL4BY ( Figure 4), evidenced by sequence and fibre-FISH data for fosmids. The SRY gene itself is present in two head to head copies, validated by fibre-FISH experiments using the gene sequence directly as a probe. There are unlikely to be more than two SRY copies on the chromosome (qPCR -supplementary figure S10). In both palindromes, the ancestral arm of the palindrome was the proximal copy, and the derived copy is distal. This is evidenced by a disrupted ancestral ERV element outside the SRY inversion, and a partially duplicated LINE element overlapping the CUL4BY inversion (see supplementary figure S9). The pattern of markers at the breakpoint regions reveals that the SRY duplication preceded the CUL4B duplication (see Figure 4 and supplementary figure S9). The transposable elements at both inversion boundaries are annotated as specific to the Sus lineage, suggesting these are relatively recent duplications. The copy number of SRY and CUL4B in closely-related suid species is therefore uncertain. The arms of the palindrome have high levels of sequence identity; we could not detect a difference in the SRY sequence from clones on one arm versus the other arm. Our sequence contigs do not cover the centres of the palindromes (about 20kb is missing in each), so we do not know if the arms abut -it is possible that there is a short stretch of unique sequence between the arms.
A further palindromic region lies at the proximal end of USP9Y, covering approximately 50kb, including the final exons 18-43 of the gene. The breakpoints lie well within known transposable elements; one end within a LINE, and the other end within a SINE, and these elements likely facilitated the original rearrangement (supplementary figure S11). Compared to SRY and CUL4BY, the breakpoints are less well defined, with sequence identity decreasing over some tens of base pairs; this probably reflects a more ancient event than SRY and CUL4BY.

Transposition within the Y chromosome
A distinctive pattern of repetitive elements covering 5kb is found in two Y contigs with 99% sequence identity. This likely resulted from a transposition of sequence at end ERV elements.
The sequence originated near SMCY, and transposed to OFD1; a cluster of ERVs can be seen at the distal end of OFD1Y in the distal gene block, and a further ERV cluster distal to TXLNGY in the proximal gene block (visible in Figure 3, higher resolution in Figure S2). The breakpoint regions are clearly visible within one ERV element (ERV1-4_SSC-I-int) upon alignment (see supplementary figure S5), while the breakpoint at the OFD1Y distal end has been further disrupted by an L1 LINE element insertion. These repetitive elements are all annotated as being Sus lineage specific, and show evidence for recent dynamic activity of transposable elements on the Y chromosome.

Structural rearrangements compared to other species
We generated likely pathways of gene-only rearrangements from the ancestral Y chromosome to pig, using data from (Li et al. 2013a) updated with pig and cattle (Elsik et al. 2009) Y chromosomes (supplementary figure S12). In pig, as in other species, the USP9Y-DBY-UTY genes are the only cluster to retain proximity from the ancestral Y chromosome. Global alignments of chromosome content with other available Y chromosome assemblies are also shown in Figure 1. This shows the striking rate of Y-chromosomal change across mammals.
Little of the Y sequences align, save for the genic regions, and these show gross levels of rearrangement. In particular, the comparison highlights some of the ancestral gene families that have become amplified in different lineages, e.g. HSFY in pigs and TSPY in dogs.

Discussion
We present here an updated and substantially improved assembly of the pig X chromosome, and a first-generation sequence and map of the pig Y chromosome with expression analysis of many of the genes it contains. This sequence has also allowed us to recover information on the evolution of the Y chromosome, and how this relates to sex chromosome evolution in other mammals.

An updated assembly and annotation of the porcine X chromosome
The picture of mammalian X chromosomes is one of general structural stability; see Figure 1.
Specific lineages, such as rodents, have many more X rearrangements than the others, but these species are characterised by a globally higher number of chromosomal rearrangements (Stanyon et al. 1999). Apparent inversions and translocations in the pig X, relative to the ancestral X, detected in previous builds are corrected here to an order more reflective of the inferred ancestral X chromosome. Similar findings may be seen with other mammalian X chromosomes as the quality of the assemblies improve. It paints a stark contrast to the dynamic and ever-changing mammalian Y chromosomes that we discuss below.

Optical mapping assists in chromosomal assembly
The assembly of the chromosome was assisted by the use of optical mapping technology to determine clone order and orientation, and to estimate gap sizes. Optical mapping has been used mainly for small genomes (e.g. Riley et al. 2011), and its use is growing for larger chromosome-level assemblies. Optical mapping assisted with the tomato genome assembly (Shearer et al. 2014), and the recent domestic goat (Capra hircus) genome (Dong et al. 2013) was assembled with a combination of short -read sequencing and automated optical mapping that produced super-scaffolds 5 times longer than those generated by conventional fosmid end sequencing. In this case, comparing the previous 10.2 X build with our assembly shows the regions of the chromosome that were particularly helped by the optical mapping approach (see supplementary figure S1). Unsurprisingly, these are the regions enriched for repetitive content; the pericentric heterochromatin especially contained a number of clones in which discrepancies in the sequence could be resolved.

Improved gene annotation of the porcine X
The revised gene annotation increases the number of protein-coding genes identified on the pig X to 689, bringing the reported gene content closer to that identified in the X chromosomes of well-studied species (i.e. humans and mice, with 813 and 940 protein-coding genes respectively). The majority of X-chromosome genes are shared between species (76% of annotated pig genes shared with human), in accordance with Ohno's law (Ohno 1967). We have highlighted some specific genes of interest in the Results section with an updated status from build 10.2 X. Beyond this, the improved gene annotation and ordering will be useful not just to the scientific community involved in mammalian sex chromosome evolution, but also to researchers using pigs as a model system for understanding disease. Many loci associated with mental retardation and brain function are found on the human X chromosome (Gécz et al. 2009). We have previously shown an association of X-linked QTLs associated with maternal aggression in sows towards their young with puerperal psychosis in humans (Quilter et al. 2007(Quilter et al. , 2012; this updated map of the pig X will make detecting significant loci and understanding their contribution to phenotypic traits far more robust in the future. These data also have relevance to the animal breeding industry; many of the QTLs mapped to the X chromosome are associated with production traits (e.g. average daily gain), meat quality traits (e.g. back fat), carcass traits (e.g. percentage lean meat) and reproduction traits (e.g. teat number).
Apart from the olfactory receptor gene clusters, we have not found evidence for widespread ampliconic gene families on the pig X. This contrasts with the X chromosomes of other species; both human and mouse X chromosomes contain independently amplified gene families, with little overlap between the species (Mueller et al. 2013). Human X chromosomes contain multiple inverted repeats with high sequence identity, enriched for testis-expressed genes (Warburton et al. 2004). Mice have a greater number of X-ampliconic genes than humans, apparently driven by a genomic conflict between X and Y sequences; the gene Sly on mouse Yq represses gene expression from sex chromosomes in spermatids, and the copy number of X genes has increased in response to maintain expression of key genes (Ellis et al. 2011). These examples led to an expectation that this might be a general feature of mammalian X chromosomes, and that the pig X would also contain unique ampliconic testis-expressed genes, but this appears not to be the case.

The porcine Y chromosome
One of the striking aspects of the Y chromosome organisation is that the single copy malespecific genes are found in tight clusters of contigs spanning only two or three megabases of sequence. This is a pattern observed in other mammalian Y chromosomes -for example humans (Skaletsky et al. 2003), mice (Soh et al. 2014), cattle (Elsik et al. 2009), cats and dogs (Li et al. 2013a). There is now considerable evidence that each lineage has preserved a small region of key ancestral X genes, and the remainder of the chromosome has evolved in a species-specific manner, often involving expansion of sequences either ancestral or newly acquired into ampliconic tracts.

Palindromic sequences
A recurring feature of the Y sequences we have assembled is the presence of palindromic regions, each on the order of 120kb end to end. These have resulted from duplication and inversion of sequences, and at least three such palindromes are visible. Two have very high levels of sequence identity; the inverted structure may facilitate the process of gene conversion, by allowing the formation of stem-loop structures as seen in the palindromes on the human Y chromosome (e.g. Hallast et al. 2013).
The first of these palindromes is in the proximal gene block, encompassing the two copies of the male-determining gene SRY. Multiple SRY copies are found in dog (Li et al. 2013a) and some rodent species (e.g. Prokop et al. 2013;Lundrigan and Tucker 1997;Murata et al. 2010), but there has previously been no suggestion of this being the case in the pig. Prior sequencing of the SRY gene and comparison across different pig breeds has given no evidence for polymorphisms in the recovered SRY sequence from any individual, despite there being breedspecific differences -i.e. there are no heterozygous pigs identified (see >300 sequences deposited in NCBI for S. scrofa alone). This suggests that this region undergoes frequent gene conversion that maintains the sequence identity between the copies. While most other species with multiple SRY copies have identifiable sequence differences, there is also a known example in rabbits of a palindrome encompassing SRY, with gene conversion maintaining sequence identity (Geraldes et al. 2010).
Close to SRY lies the second palindromic region, surrounding the CUL4B (cullin) fragments.
The palindromic region is of similar size to the SRY region, and in fact overlaps the SRY palindrome. This allowed us to determine that the duplication around SRY occurred first, bringing with it the first two cullin exons, after which the cullin fragment was duplicated again (Supplementary Figure S9).
The third example of this phenomenon is found in the distal gene cluster. A palindrome includes the latter 26 exons of the USP9Y gene, though without forming an open reading frame. Unlike the previous two palindromes, both breakpoint ends lie within complete transposable elements (TEs). These types of sequence have long been known to facilitate genomic rearrangements via processes including non-allelic recombination and non-homologous end joining (Stankiewicz and Lupski 2002;Hastings et al. 2009). In the case of USP9Y, sequence identity between palindrome arms is lower around the breakpoints, perhaps indicating the duplication results from an older event. In all three palindromes, the TEs are annotated as deriving from within the pig lineage -that is, these are not ancient repeat elements. This tells us that these palindromes have arisen relatively recently, and show the ongoing impact of repetitive content in the genome. However, we need data from more species to determine when the duplications occurred. Extant suids diverged since about 40Mya, and the copy number of the genes involved across these lineages remains to be identified.
The palindrome structures are reminiscent of ampliconic structures found on the mouse (Soh et al. 2014) and human Y chromosomes ( (Hughes et al. 2012). There are likely other similar palindromic sequences remaining to be discovered within the Y, and more examples of the impact of transposable elements that have promoted accelerated rearrangements will no doubt be found.

Ampliconic sequences
Most mapped mammalian Y chromosomes have been found to contain multi-copy gene families (e.g. Li et al. 2013a), and the pig is no exception. Besides the duplications due to palindromes, other sequences have amplified to a much greater extent. There are three gene families of note here.
The CUL4B fragments: in addition to the two fragments noted above, further cullin fragments are found proximal to SRY and the active cullin copy, and appear to form part of an ampliconic and fragmented region leading towards the centromere. Cullins are ubiquitin ligase genes, frequently found in multiple copies on mammalian Y chromosomes (Li et al. 2013a;Murphy et al. 2006), and are important for appropriate degradation of substrates during many process including gamete production (Sutovsky 2003).
The TSPY fragments: Other amplifications have been defined for TSPY, based on FISH mapping of clones. These again appear to be interspersed in the region leading towards the centromere but it not clear how these are arranged. TSPY is an ampliconic gene in many mammalian species, from artiodactyls to primates (Xue and Tyler-Smith 2011); the genes are involved spermatogenesis (Xue and Tyler-Smith 2011), and, in cattle, copy number variation of this gene has been linked with fertility in bulls (Hamilton et al. 2012).
The HSFY family: These genes are predominantly found in a block between the proximal and distal gene clusters. The organisation and content of this region is complex, and we have investigated it in more detail in a companion paper (Skinner et al, in submission). The HSFY genes show evidence for a recent amplification in the Sus lineage to ~100 copies, perhaps independent amplification in other suid species, and further independent amplification in cattle (Chang et al. 2013).
All of these pig ampliconic genes are involved in amplifications in other species, and the known functions suggest they are important in spermatogenesis.

Other amplified sequences
Yq is dominated by repeat sequences (as demonstrated by the painting pattern of FISH using certain BAC and fosmid clones). These clones are composed almost entirely of sequences related to male-specific repeats described previously for pig (e.g. Mileham et al. 1988), and thus a more detailed study is needed to understand the organisation of this arm of the chromosome.
At this time we can say only that patterns of hybridisation from FISH suggest to us that there is single copy sequence on Yq, including near/at the Yq terminus, but we were unable to isolate the sequences involved.

Comparative chromosome organisation and gene order between mammals
Previous work from primates, mouse, cat and dog has reconstructed a putative ancestral eutherian Y chromosome (Li et al. 2013a) based on gene order. We have incorporated our pig gene organisation into this, and added information available in the cow Y sequence assembly (Elsik et al. 2009) as shown in supplementary figure S12. One group of genes stands out from the comparison: USP9Y-DBY-UTY is the only ancestral cluster of genes that have retained their proximity to each other in all the studied species. There may be a selective disadvantage to disrupting this organisation, as has been proposed for other conserved syntenic blocks in general (Larkin et al. 2009), and for these genes in particular. Both USP9Y and DBY have been implicated as important in human spermatogenesis, though they may not be essential in all great apes (Tyler-Smith 2008).

TRAPPC2P-OFD1Y-GPM6B: A potential transposition from the X chromosome
Outside the PAR, there are regions of homologous sequence between the X and Y chromosomes. These can be as small as individual exons of a gene, or multi kilobase regions.
We found most of the XY homologies could be attributed to known XY homologous genes, or to otherwise unannotated repetitive sequences, such as ERV families enriched on the sex chromosomes. The large exception was the region TRAPPC2-OFD1-GPM6B, which has been subject to recurrent transposition onto the Y from the X chromosome in dog (Li et al. 2013a) and cattle (Chang et al. 2011)

lineages. A similar transposition affecting the RAB9A-SEDL-OFD1Y
genes has occurred in the human and chimp lineage (Hughes et al. 2010). This gene cluster is also at least partially present in the pig Y, with good sequence identity to the X across most of the region. Most discrepancies to the X sequence can be explained by recent transposable element activity (see Supplementary Figure S6 for X-Y alignment and content of the region).
The most likely explanations for the presence of this region in pigs, cattle and dogs are (1) a recurrent transposition of a similar region independently in dogs, pigs and cattle, or (2) in the ancestor of pigs and cattle and independently in dogs. More sequence information in pig is needed to be able to compare the region with cattle to determine if the breakpoints are shared.
Our analysis in pig is also complicated by an apparent further rearrangement affecting the region and splitting GPM6B, the first exon of GPM6B being proximal to UBE1Y. OFD1 is involved in cilia formation and the ubiquitin-proteosome degradation pathway. Ubiquitin degradation is an important part of gamete development, as seen in the CUL4B genes, and ciliopathies have been previously implicated in fertility issues (Fry et al. 2014). It is plausible that the testis-expressed OFD1Y has repeatedly acquired a function in sperm development in different mammalian lineages; certainly the X copies of OFD1 and also CUL4B have been found to be substantially downregulated in teratozoospermic men (Platts et al. 2007).

Conclusion
This work presents an improvement to the pig X chromosome assembly and gene annotation, plus the first sequence assembly for the pig Y chromosome, produced through a combination of fosmid sequencing, physical mapping by FISH and manual assembly. The assemblies we have generated have allowed new insights into the content and evolution of the pig sex chromosomes, and provide an important resource for the pig genomics community.

Acknowledgements
We wish to thank Genus for providing the Duroc boar from which a cell line was established and tissue samples obtained. We gratefully acknowledge the Wellcome Trust Sanger Institute core teams for fingerprinting, mapping, archiving, library construction, sequence improvement and sequencing which underpins this work. This work was funded by BBSRC grant BB/F021372/1.

Comprehensive methods are given in the supplementary material Library construction and sequencing
Phytohaemagglutinin-stimulated peripheral blood culture from a Duroc boar was used to prepare chromosomes for flow sorting. Flow-sorted Y chromosomes were purified, and 30-50kb sized fragments were cloned into the pCC1Fos vector (library WTSI_1061: http://www.ncbi.nlm.nih.gov/clone/library/genomic/330/). Clones for sequencing were targeted by minimal overlapping clones on a fingerprint contig (FPC) map. The targeted 897 clones for the Y chromosome were sequenced using a combination of 3 different sequencing platforms: Capillary, Illumina and 454. Raw sequence data is submitted to the public data repository, ENA http://www.ebi.ac.uk/ena/, under accession number ERP001277. Clones were assembled using a combination of four assembly scripts to produce de novo assemblies. All clone sequences are submitted to ENA. Manual alignment of clone sequences was used to build the clone map, expanding from clones containing known genes. These small contigs were oriented and ordered using fibre-FISH on single DNA-molecule fibres. Final sequence contigs were assembled based on this map using GAP5 (Bonfield and Whitwham 2010).

Molecular combing and FISH
Single-molecule DNA fibres were prepared by molecular combing (Michalet et al. 1997). Purified fosmid DNA was amplified using a GenomePlex® Whole genome Amplification (WGA) kit, and labelled using a WGA reamplification kit (Sigma-Aldrich) using custom-made dNTP mix as described before (Gribble et al. 2013). Fluorescence in-situ hybridisation followed standard protocols. Probes were detected with fluorescently conjugated antibodies. Slides were mounted with SlowFade Gold® mounting solution containing 4',6-diamidino-2-phenylindole (Molecular Probes/Invitrogen) and visualised on a Zeiss AxioImager D1 microscope. Digital image capture and processing were carried out using the SmartCapture software (Digital Scientific UK).

X and Y gene annotation, sequence content and chromosomal evolution
Manual annotation on the pig X and Y chromosomes was performed using the Otterlace/Zmap suite of annotation tools (Loveland et al. 2012) following previously established annotation protocols Dawson et al. 2013). The assembled chromosomes were run through an annotation pipeline (Searle et al. 2004), aligning EST, mRNA and protein libraries against the chromosomes with all annotated gene structures (transcripts) supported by at least one form of this transcriptional evidence. The HUGO Gene Nomenclature Committee (HGNC) (Seal et al. 2011) naming convention was used whenever possible for all pig genes, else HAVANA naming conventions (http://www.sanger.ac.uk/research/projects/vertebrategenome/havana) were followed.
RepeatMasker (Smit et al. 1996) was used to identify repetitive elements within Y contigs. Gene content and structure was identified both by gene prediction program GeneMark (http://exon.gatech.edu/GeneMark/) and manual comparison to known Y genes using BLAST.
Targeted resequencing was performed across specific genes to confirm their structure (primers given in supplementary table S4). Regions of XY homology were identified by comparing the repeatmasked X assembly to all sequenced repeatmasked Y clones (mapped and unmapped) using LASTZ (Harris 2007). Unannotated repetitive content enriched on the sex chromosomes was minimised by performing an X-X self alignment, and subtracting multiply-hit regions from the X-Y comparison. Candidate regions of 1kb or longer were further analysed using BLAST for similar sequences in the NCBI NT and EST databases, and the VEGA X gene annotations.
Evolutionary analyses between X and Y gene pairs were conducted using the Nei-Gojobori model (Nei and Gojobori 1986) in MEGA5 (Tamura et al. 2011). For each pair, positions containing gaps and missing data were eliminated. Reconstruction of ancestral Y chromosome organisations was performed using the Multiple Genomes Rearrangement (MGR) program (Bourque and Pevzner 2002) to calculate optimal rearrangement pathways between each species, as previously described (Skinner and Griffin 2012).

Gene expression
RT-PCR was used to confirm expression status of selected genes in five tissues (brain, liver, kidney, side muscle, testis), obtained from the same boar from which blood cultures were derived. Samples were taken from tissues stored in RNAlater (Qiagen) and homogenised in Trizol. Nucleic acids were extracted with phenol-chloroform and DNase treated. RNA was precipitated with isopropanol and stored at 1µg/µl in ddH 2 O at -80°C. RT-PCR was carried out using a OneStep RT-PCR kit (Qiagen) on 25 ng of total RNA. Primer sequences are given in supplementary table S5.

Copy number estimation of SRY by qPCR
Primers were designed to amplify a 1447bp region across the majority of the SRY ORF and UTRs (F:TAATGGCCGAAAGGAAAGG; R:TGGCTAATCACGGGAACAAC), and products were generated using a MyTaq Red kit (Bioline) using the profile 95°C for 3mins, 35 cycles of 95°C/53°C/72°C for 15s/15s/2min with a final extension of 72°C for 10min. Two female Duroc gDNAs were spiked with dilutions of the purified SRY product to give a standard curve of 4 copies SRY per genome to 0.25 copies per genome (assuming diploid genome size of 6Gb; Animal Genome Size Database; Gregory 2006). qPCR was performed using a SYPR-FAST qPCR kit (Kapa Biosystems) on the spiked females and on 5 male Duroc gDNAs with primers for SRY and the autosomal (SSC10) gene NEK7 (supplementary table S6). Annealing temperature was optimised at 57°C. Cycling conditions were 95°C for 3mins, followed by 40 cycles of 95°C/ 57°C/72°C for 10s/20s/30s. The fluorescent signal threshold crossing point (Ct) was normalized to the average signal from Nek7 to produce a normalised ∆C t . The data from spiked female gDNA was used to construct a standard curve relating SRY signal to Nek7 signal; from this, an estimate of the absolute SRY copy number in the male gDNA samples was produced.

Data access
All sequence and annotation is available via the Vega genome browser (http://vega.sanger.ac.uk/index.html). The pseudoautosomal region of X/Y homology between the X and Y chromosomes is represented on the X chromosome only in Vega and Ensembl. It is marked as an assembly exception in both chromosomes, but the underlying genomic sequence and annotation is that of the X chromosome. Only the unique regions of the Y chromosome are stored and annotated. The complete Y chromosome is represented by filling the 'gaps' with the PAR regions from the X chromosome. Pigs have an expansion of olfactory receptors throughout the genome; two clusters lie on the pig X, not found in humans. The region provides an example of the sequence improvement from optical mapping between the 10.2 assembly (top; image from Ensembl) and the current assembly (centre), with a comparison to the corresponding region on the human X (bottom; both images from Vega). The improvement in the current assembly allowed reorientation (Ensembl image has been flipped for clarity), the correct positioning of the VAMP7 gene which was previously annotated on an unplaced scaffold and improved annotation of the genes (blue) and pseudogenes (pink) within the region.

Figure 3 -The pig Y chromosome
The organisation of the pig Y chromosome. All identified male-specific single copy genes are on the short arm and split into two blocks by the ampliconic HSFY region. The genes (blue) and pseudogenes (grey) are shown within each block. To the left of the figure is the density of repetitive content based on a 10kb sliding window within the sequence contigs we obtained.

Figure 4 -The pig SRY region
The Yp proximal block of genes contains two overlapping palindromes of about 120kb each.
These surround the duplicated sequences CUL4BY exons 1-2 and SRY. A) FISH results from Y fosmid clones and probes for the SRY gene are shown with the BAC and fosmid clone sequences found mapping to the region. The inversion boundaries are both identifiable; the CUL4BY inversion runs from the last 3kb of 43B21 to within 72J17; the SRY inversion begins also within 72J17 and runs to 13A15. A schematic view is also shown of the regions surrounding B) the SRY and C) the CUL4B duplications. The SRY duplication disrupts and ERV element, revealing the proximal copy to be ancestral. The CUL4B duplication copies part of a LINE element, again revealing the proximal copy to be ancestral. The sequence alignments across the inversion breakpoints are shown in more detail for D) SRY and E) CUL4B. The order of events was therefore a duplication of SRY including the first two exons of CUL4B, followed by duplication of the region around the CUL4B copy. See also Figure S1 for broader context of the region. Tables   Table 1 -Comparison of gene content between porcine build 10.2 X, our updated build of the X (X-17), and the new Y annotation (Y-12). The number of identified coding genes has been substantially increased, and brings the porcine X closer in gene content to other well-sequenced mammalian X chromosomes (e.g. human and mouse).