Z Chromosome Detection from the Male Emu Genome (ZZ) Assembly

As shown in figure 1, 51 previously reported scaffolds from the ostrich genome, with a total length of 88.8 Mb, were used as scaffolds for the Z chromosome (Zhang et al. 2015; O’Connor et al. 2018; Xu, Wa Sin, et al. 2019). We constructed a draft of the entire Z chromosome for the cassowary that consisted of 106 scaffolds with a total length of 85.6 Mb (Sackton et al. 2019). Similarly, we constructed the counterpart region for the emu, which consisted of 13 scaffolds with a total length of 82.4 Mb. Using data from an emu male individual (ZZ) for the assembly, we succeeded in constructing six continuous scaffolds with a total of 55.3 Mb (PAR), one scaffold of 4.9 Mb (S1 region), and eight scaffolds of 22.2 Mb (S0 region). Moreover, one long scaffold, labeled scaffold1724899_cov79, covered part of the PAR (39.1 Mb), the entire S1 region (4.9 Mb), and part of the S0 region (7.7 Mb).

Open in a separate window

Fig. 1.

Overview of the Z chromosome with the boundary between a pseudoautosomal region (PAR) and stratum 1 (S1) region. The whole Z chromosomes of emu, cassowary, ostrich, and chicken are illustrated. Gray horizontal lines represent the scaffolds derived from the Z chromosome of each species. The order of the scaffolds in the ostrich is based on previous data (Zhang et al. 2015), and that in the cassowary and emu is based on an alignment with the chicken and ostrich Z chromosomes. The color map represents sequence alignments and their identities between the two species. The histograms shown on the top of each emu, cassowary, and ostrich scaffold represent normalized read coverage based on the sequence depth of the Illumina paired-end (PE) reads. Dark gray boxes on the histograms indicate the pseudoautosomal region (PAR) for each species. The red vertical line indicates the boundary between the PAR and S1 region of the emu and cassowary, and the blue line indicates the boundary between the PAR and S1 region in the ostrich. These graphs were plotted using an in-house python script with Matplotlib.

PAR Estimation from Read Mapping Data

The boundaries between the PARs and S1 regions were detected using the normalized read depth for each female (ZW) emu, cassowary, and ostrich. The normalized read depth along the Z chromosome clearly changed from two to one, suggesting that the upper region corresponded to the PAR (Z-W recombination); in contrast, the downstream region corresponded to the W-degenerated (Z-specific) region (fig. 1). Therefore, we defined the breakpoint of the read depth changes as the boundary between the PAR and S1 region in each species.

Comparisons of the loci of the breakpoints and synteny among the three species showed that the boundary was consistent between the emu and cassowary (Casuariiformes; indicated by the red vertical line in fig. 1), whereas the boundary in the ostrich chromosome (blue vertical line) was in a distinct position. We further confirmed that the ostrich Z chromosome had a relatively large S1 region that expanded approximately 2–3 Mb toward the PAR side compared with those of the emu and cassowary, consistent with previously studies (Zhou et al. 2014; Xu, Wa Sin, et al. 2019; Xu and Zhou 2020).

Phased De Novo Assembly of the Female Genomes

To construct phased S1 sequences, we initially applied Platanus-allee v2.2.2, a de novo haplotype assembler, to one individual genome for each female emu, cassowary, and ostrich (table 1). Read data from the female emu are shown in supplementary table 1, Supplementary Material online, and those of cassowary and ostrich were downloaded from the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database and are presented in supplementary table 3, Supplementary Material online. In the female emu, the total length of the assembled consensus sequences was 1.30 Gb, which was almost equal to the estimated genome size. The total length of the primary haplotype-phased sequences (referred to as “bubble” scaffolds in Platanus-allee) corresponded to the longer haplotype-phased sequences in a pair, which was 682 Mb with an N50 scaffold length of 128 kb. In contrast, the total length of the cassowary consensus sequences was 1.22 Gb, and the total length of the primary bubble sequences was shorter than that of the emu at 131 Mb, with an N50 scaffold length of 15 kb. The estimated heterozygosity rate in the cassowary genome was lower than that in the emu, suggesting a link between heterozygosity and total size of the phased regions (supplementary table 4, Supplementary Material online). Additionally, this indicated that the libraries used for de novo assembly contributed to the continuity of the scaffolds. Specifically, in the emu, four mate-pair (MP) libraries with an insert size of 3–15 kb were used for de novo assembly, whereas only the 3-kbp MP library was used for the cassowary assembly (supplementary tables 1 and 3, Supplementary Material online). By contrast, the total length of the ostrich consensus sequences was 1.29 Gb and that of the primary bubble sequences was 770 Mb, which were slightly longer than that in the emu and had an N50 scaffold length of 31 kb. Detailed assembly results and comparison of our results with other assemblies are shown in supplementary results and supplementary table 2, Supplementary Material online.

Construction of Phased S1 Regions

Following the haplotype assembly for female emu, cassowary, and ostrich, we detected bubble scaffolds derived from the S1 region and assigned each of the phased scaffolds to the Z or W chromosome (table 2 and supplementary fig. 1, Supplementary Material online). We verified the accuracy of our haplotype-resolved assembly in female individuals using two methods. First, we confirmed the high sequence identity between the male emu (ZZ) S1 region and the female emu Z-linked S1 region according to dot-plot alignments (supplementary fig. 2, Supplementary Material online). Second, we confirmed that both the Z-linked and W-linked S1 regions had high sequence identity with those of the PacBio-based female emu assembly (fig. 2A; ZJU1.0; Liu et al. 2021). The results of dot-plot with another Illumina-based study (ASM1339679v1; Feng et al. 2020; Xu and Zhou 2020) are shown in figure 2B. However, it is clear that this Illumina-based study did not provide sufficient phased assembly, as only a few regions were constructed as separate phased sequences into Z-linked and W-linked scaffolds (fig. 2B, enlarged image), and most regions were constructed as only one sequence. The sum length of the previously reported assembly was only approximately 5.5 Mb, whereas the sum length of our six scaffolds was approximately 9.9 Mb (4.9 Mb consisting of three Z-linked scaffolds and 5.0 Mb consisting of three W-linked scaffolds), which was the same as the length of the PacBio-based assembly results in the S1 region. In addition, the published Illumina-based assembly had approximately 96–98% homology with both of our W-linked and Z-linked sequences compared with PacBio assembly, which had almost 100% homology between the same phased sequences and approximately 94–96% homology between the opposite phased sequences. This result strongly indicated that both our assembly and PacBio-based assembly succeeded in obtaining fully phased S1 sequences, whereas the previous Illumina-based assembly consisted of a mosaic-like mixture of Z-linked and W-linked sequences. Similar results were obtained from analyses of the assembly results from cassowary compared with previous studies (supplementary fig. 3, Supplementary Material online). These results confirmed that our assembly has high accuracy and comprehensiveness (table 2 and fig. 2).

Open in a separate window

Fig. 2.

Dot-plot alignment of the female emu S1 region. (A) Dot-plot alignment of the female emu S1 regions between the Platanus-allee assembly constructed in this study and the PacBio-based assembly (ZJU1.0) (Liu et al. 2021). The horizontal axis shows the Z-linked and W-linked S1 region sequences from the female emu that was constructed by Platanus-allee assembler in this study. The vertical axis shows the Z-linked and W-linked S1 region sequences of the female emu constructed from the PacBio-based data (ZJU1.0). (B) Dot-plot alignment of female emu S1 regions between Platanus-allee and Illumina-based (droNov1) assemblies (Sackton et al. 2019). The horizontal axis shows the Z-linked and W-linked S1 region sequences in female emu constructed by the Platanus-allee assembler. The vertical axis shows the S1 region sequences in female emu constructed from Illumina-based data by the ALLPATHS-LG assembler (droNov1). Magnified views of the two regions are shown in the dot-plot alignments framed in pink and light blue, respectively.

Comparisons of Phased S1 Regions

Sequence alignments between the Z-linked S1 and W-linked S1 regions within species (gametologous sequences) are presented in figure 3A. We noted that these results were based on the alignment of corresponding regions to show the sequence homology between Z-linked S1 and W-linked S1 regions, which were not anchored to chromosomes and do not necessarily represent their positions on chromosomes. The detected S1 regions of the emu and cassowary were 4.9–5.1 Mb in length and nearly in the same loci in both species. Conversely, the S1 region constructed from phased scaffolds in the ostrich was approximately 7.2 Mb, with nine scaffolds in total that were approximately 2.2 Mb longer toward the PAR side than that in emu and cassowary (fig. 3B). Moreover, within the same species, average sequence differences between gametologous sequences in the S1 region were markedly higher than the heterozygosity rates (i.e., the average sequence difference rate of the whole genome), which were estimated with GenomeScope (supplementary table 4, Supplementary Material online). The sequence identities between gametologous sequences in the S1 regions of the emu and cassowary were approximately 95.9% and 96.2%, respectively, which were slightly higher than that reported previously (Xu and Zhou 2020). Gametologous sequence divergence, described as the density of heterozygous sites in the S1 region, slightly increased from the S0 region side toward the PAR side (fig. 3C). Whereas the sequence identity between gametologous sequences in the ostrich S1 region was 89.1%. These results suggested that there was a higher level of differentiation in the S1 region in the ostrich than in the emu and cassowary. Interspersed repetitive sequence landscapes in the Z-linked and W-linked S1 regions of emu and cassowary females showed that the expansion of the repetitive elements, especially for LTR families, occurred in the emu W-linked region after Z/W differentiation and emu/cassowary differentiation. The same insertion pattern of young repetitive sequences was also observed in the ostrich, however, no traces of recent repetitive sequence accumulation in the W-linked region of the cassowary were observed. (fig. 3D and supplementary figs. 4 and 5, Supplementary Material online).

Open in a separate window

Fig. 3.

Overview of the phased Z-linked and W-linked S1 regions and their features. (A) Phased S1 scaffolds and gene synteny. Gray horizontal lines represent phase scaffolds derived from the S1 region from each species. The color map indicates sequence identity between Z-linked and W-linked S1 gametologs. The arrows on the scaffolds indicate genes annotated in the S1 regions. (B) Schematic diagram of sex chromosomes from each species. The pseudoautosomal (PAR), Z-linked S1, and W-linked S1 regions are labeled in gray, light blue, and light orange, respectively. The Z-linked S0 and W-linked S0 regions are labeled in blue and orange, respectively. Red vertical lines show the boundaries between the PAR and S1 regions, and between the S1 and S0 regions. (C) Distribution of gametologous heterozygous sites in the S1 region of a female emu individual. (D) Interspersed repetitive sequence landscapes in the Z-linked S1 region (left) and W-linked S1 region (right) of a female emu individual. The horizontal axis represents the Kimura’s substitution rate from the consensus sequences of each repetitive unit. Older copies are located on the right of the graphs, whereas more recent copies are located on the left.

Emu Genome Sequencing

To obtain whole-genome sequencing (WGS) reads from two emu individuals (a male and a female), we constructed paired-end (PE) libraries and four mate-pair (MP) libraries for the Illumina platform (Hayward, CA, USA). PE libraries with an insert size of 600bp were prepared using a TrueSeq PCR-free DNA Sample Prep kit (Illumina). For MP libraries, we extracted approximately 3, 6, 10, and 15 kb DNA fragments from an agarose gel for each sample and prepared the libraries using a Nextera Mate Pair Library Prep kit (Illumina). We then confirmed the samples by sequencing a small part of the PE libraries with a read length of 250bp using an Illumina MiSeq platform. We then sequenced all the libraries with a read length of 150bp using an Illumina HiSeq X Ten platform. The sequence data of each sample are presented in supplementary table 1, Supplementary Material online. The basic statistical parameters of the emu genomes (e.g., estimated genome size and heterozygosity) from these sequence data were calculated using the method described in the supplementary methods, and the results are shown in supplementary table 4 and figure S7, Supplementary Material online.

Preparation of Raw Next-Generation Sequencing Data of Ostrich and Cassowary

Female ostrich and cassowary WGS data were downloaded from the NCBI SRA database. For the ostrich sequence data (NCBI BioProject number PRJNA212875), we used three Illumina PE libraries with insert sizes of 170, 500, and 800bp and four MP libraries with insert sizes of 2, 5, 10, and 20 kb. For the cassowary sequence data (NCBI BioProject number PRJNA433110), a PE library with an insert size of 220bp and an MP library with an insert size of 3 kb were used. The sequence data for each sample are presented in supplementary table 3, Supplementary Material online. The basic genomic features were estimated with the method described in the supplementary methods, and the results are shown in supplementary table 4 and figure 7, Supplementary Material online.

De Novo Assembly of the Male Emu Genome

The sequence reads of a male emu individual (ZZ) were assembled using Platanus v1.2.4 software (Kajitani et al. 2014). We used an assembler different from that used for the female individuals (ZW), prioritizing long connections because the goal was to not conduct phased assembly. Contig assembly was performed using only the PE library, following which scaffolding and gap closing were performed using all libraries. Additionally, we confirmed that scaffolds ≥500bp from the contamination were not included by performing alignments with the NCBI nr and nt databases using Diamond v0.9.22 (Buchfink et al. 2015) and Blastn 2.2.27+ (Altschul et al. 1997), respectively. For the subsequent analyses, we used a male emu draft genome with a complete mitochondrial genome. The assembly of the mitochondrial genome is described in the supplementary methods, Supplementary Material online

Z Chromosome Genome Detection

For ostrich, we first downloaded the whole-genome assembly (GCF_000698965.1) and then picked and ordered the scaffolds corresponding to the ostrich pseudo Z chromosome according to previously published reports (Zhang et al. 2015; O’Connor et al. 2018). For emu and cassowary, we initially detected Z chromosome scaffolds based on the similarity with the chicken Z chromosome (Gallus_gallus-5.0, NC_006127.4) independent of the ostrich Z chromosome and then ordered the scaffolds along the ostrich Z chromosome. In emu, the male (ZZ) draft genome was assembled with Platanus v1.2.4 and aligned to the chicken Z chromosome using LASTZ v.1.02.00 (Harris 2007). Then, scaffolds with ≥100 kb and more than 50% of the aligned length were detected as the emu Z chromosome. Finally, scaffolds were ordered and oriented according to the alignment with the ostrich Z chromosome. For the cassowary, we used the draft genome GCA_003342895.1 and subsequent methods described above for the emu genome construction.

PAR Estimation from Read Mapping Data

To estimate the PAR of each species, we calculated read depths normalized to one per haploid set for every 10-kb segment on the Z chromosome using female PE reads. Given that female individuals have Z and W chromosomes, the PAR is expected to be present in two copies (diploid), whereas the W degenerated region is expected to be present in one copy (haploid). The PAR was estimated by detecting changes in the copies along the Z chromosome. PE reads were mapped on each draft genome using BWA aligner v.0.7.12 (Li and Durbin 2009), and the sequence read depth was calculated using SAMtools v.1.8 (Li et al. 2009) based on the reads mapped with ≥98% identity and ≥80% alignment cover rates. Average read depth per haploid set was obtained from the distribution of whole-genome read depth. Lastly, normalized read depth for every 10 kb of scaffold was calculated by dividing the read depth of the 10 kb segment by the average depth per haploid set.

Phased De Novo Assembly of Female Individuals

To construct the phased S1 region, which was expected to have highly divergent W/Z gametologous sequences, we used Platanus-allee v2.2.2 as a de novo haplotype assembler (Kajitani et al. 2019) for female emu, cassowary, and ostrich individuals. The “assemble” command was performed using only the PE library, whereas the “phase” and “consensus” commands were performed using all libraries. For the subsequent comparative analyses, we used “allPhasedScaffold” FASTA files that contained phased (“bubble scaffolds”) and unphased (“nonbubble scaffolds”) sequences. In Platanus-allee, the longer sequence of a pair of phased bubble sequences was defined as the “primary bubble,” whereas the remaining sequence was defined as the “secondary bubble.” A schematic explanation of the Platanus-allee result is shown in supplementary figure 8, Supplementary Material online. Repetitive sequence analysis is described in the supplementary methods, Supplementary Material online.

Construction and Comparison of Phased S1 Regions

We obtained whole genome phased scaffolds that were expected to be haplotype sequences from female WGS data using Platanus-allee. Then, we detected bubble scaffolds derived from the S1 region and assigned each of the phased scaffolds to the original Z or W chromosome. In the pairs of bubbles derived from the S1 region, with each bubble derived from either the Z or W chromosome. We assigned each of the bubble sequences to either the Z or W chromosome according to sequence identity to the male-derived genome or transcriptome such that each of the bubble sequences with high identity to the male-derived sequence was attributed to the Z chromosome. Specifically, for emu and cassowary, we aligned bubble sequences to the male emu draft genome using LASTZ. For the ostrich female individual, we used male-derived transcript contigs for the assignment of the origin.

To obtain transcripts from a male ostrich individual, we assembled RNA-seq reads (SRR1619445) from the adult cerebellum using Trinity v2.2.2 (Haas et al. 2013). Additionally, we aligned the transcript contigs ≥500bp to bubble sequences using GMAP 2018-07-04 (Wu and Watanabe 2005). Then, average alignment identities for each phased S1 scaffold were calculated using transcript contigs with ≥80% alignment cover rate to both primary and secondary bubble scaffolds (supplementary table 9, Supplementary Material online). One bubble sequence of a pair was aligned with more than 99% sequence identity to male transcript contigs, whereas the other bubble sequence was aligned with 86.5–92.0% identity. According to the respective sequence identities of the nine pairs of phased scaffolds, the bubble with the higher identity to the male transcript contig was defined as Z-linked, whereas the other was defined as W-linked.

We calculated the sequence identity between phased scaffolds in the S1 region to define the sequence identity between gametologous sequences. First, secondary bubble scaffolds were divided into 1-kb segments and the segments were aligned to primary bubble scaffolds using LASTZ (–step = 5, –notransition, –nochain). Next we obtained each of the sequence identities of the alignment following synteny detection within the gametolous sequences and measured distributions of the sequence identity for each species.

This work was supported by JSPS KAKENHI (15H04401, 16H06279 [PAGS], 19H03267, and 19H03206). We thank the staff of the Itoh Laboratory at Tokyo Tech for supporting genome sequencing and their valuable participation in scientific discussions.