Chromosome-level genome assembly of Korean native cattle and pangenome graph of 14 Bos taurus assemblies

This study presents the first chromosome-level genome assembly of Hanwoo, an indigenous Korean breed of Bos taurus taurus. This is the first genome assembly of Asian taurus breed. Also, we constructed a pangenome graph of 14 B. taurus genome assemblies. The contig N50 was over 55 Mb, the scaffold N50 was over 89 Mb and a genome completeness of 95.8%, as estimated by BUSCO using the mammalian set, indicated a high-quality assembly. 48.7% of the genome comprised various repetitive elements, including DNAs, tandem repeats, long interspersed nuclear elements, and simple repeats. A total of 27,314 protein-coding genes were identified, including 25,302 proteins with inferred gene names and 2,012 unknown proteins. The pangenome graph of 14 B. taurus autosomes revealed 528.47 Mb non-reference regions in total and 61.87 Mb Hanwoo-specific regions. Our Hanwoo assembly and pangenome graph provide valuable resources for studying B. taurus populations.


Background & Summary
Hanwoo is a native Korean taurine cattle breed with a 5000-year history as a draft animal for farming and transportation 1 .In a short period, Hanwoo underwent significant changes in its demographic history and selection.During the Korean war (1950-1953), the number of Hanwoo dropped to about 390,000, but recovered to 1.02 million by the late 1950s.With the development of the South Korean economy and agricultural industry, Hanwoo transitioned from a draft to a meat production breed in the 1960s.Modern breeding programs, including performance tests, artificial insemination and genomic selection were initiated by the South Korean government in the 1980s.These programs have improved carcass weight and meat quality of Hanwoo by increasing intramuscular fat (marbling).As a result of continuous artificial selection, Hanwoo has gained unique features both in genome and traits.
This study presents a high-quality assembly of Hanwoo which is the first chromosome-level genome assembly of Asian Bos taurus taurus using a combination of PacBio Hifi, Isoform and Illumina RNA sequencing, with scaffold N50 length of 89 Mb.The completeness of the genome was confirmed by the BUSCO score of 95.8%.The top 31 scaffolds are all greater than 17 Mb in size with a total length of 2.69 Gb. 48.7% of the Hanwoo genome is composed of various repetitive elements.The genome was annotated to contain 27,314 protein-coding genes, including 25,302 proteins with inferred gene names and 2,012 unknown proteins.
We generated a pangenome graph of 14 high-quality Bos taurus autosomes including high-quality genome assemblies of Hanwoo, Hereford, Angus, Brown Swiss, Highland, Holstein, Jersey, Original Braunvieh, Piedmontese, Simmental, Brahman, Nellore, N'Dama, and Ankole.We identified non-reference regions and breed-specific regions through the pangenome graph.In Hanwoo, 528.47 Mb of total non-reference nodes and 61.87 Mb of Hanwoo-specific nodes were identified.This pangenome graph would be used to extract structural variations and make insightful observations among various populations of Bos taurus.

Methods
Sample collection and extraction of genomic DNa and RNa.The samples used in the study of Hanwoo genome included blood, sirloin, liver, and subcutaneous fat from a steer named "bull 2050".The samples were collected from the Experimental farm of College of Agriculture and Life Sciences at Seoul National University, Pyeongchang-gun, Gangwon-do, Republic of South Korea (Fig. 1) and were approved by the Seoul National University Institutional Animal Care and Use Committee (SNU-201129-1-1).It was castrated in 9.4 months of age, slaughtered and sampled in 32 months of age.All blood sampling was carried out by trained veterinarians, according to the approved institutional protocols.Genomic DNA were extracted from whole blood using Wizard Genomic DNA Purification kit following the manufacturer's protocol.
Sirloin, liver and subcutaneous fat tissues of Hanwoo bull 2050 were collected immediately after slaughter and frozen using liquid nitrogen and stored in a deep freezer until RNA extraction.RNA was isolated using the RNeasy kits (Qiagen, Valencia, CA) following the manufacturer's protocol.
DNa library construction and sequencing.DNA sequencing libraries were prepared using SMRTbell Express Template Prep kit 2.0 (Pacific Biosciences, California, USA) and libraries larger than 20 kb were used Fig. 1 Picture of Hanwoo steer used in this study and a circos plot.Shown from the outer to inner circle are the following: gene density, with the intensity of color representing the number of genes in a 10,000 bp window; N (unknown base) ratio, with the height of the bar representing the percentage of bases that are N in a 1,000,000 bp window and the overall height of the track representing from the minimum to maximum value for the whole genome which are from 0% to 0.02%, respectively; GC content, with the height of the bar representing the percentage of GC in a 10,000 bp window and the overall height of the track representing from the minimum to maximum value for the whole genome which are from 27.07% to 74.80%, respectively; and the corresponding chromosome.
for next steps.HiFi reads were sequenced using 2 SMRT cells of 8 M Tray, Sequel II Sequencing Kit 2.0 in Pacific Biosciences (PacBio) Sequel IIe platform at NICEM in Seoul National University.Highly accurate consensus sequences were produced by PacBio CCS workflow (v 6.3.0),yielding a total of 3.5 M reads and 67.5Gbp corresponding to a genomic coverage of ~24.8X (Table 1).
RNa library construction and sequencing.For RNA-seq, paired-end libraries with insert size of 75 bp were prepared with TruSeq Stranded mRNA Sample Preparation kit (Illumina, San Diego CA USA) from total messenger RNA (mRNA) of sirloin, liver and subcutaneous fat tissues of a Hanwoo bull 2050.RNA of the three tissues were sequenced separately using Illumina NextSeq 500 with following adapters; liver: D701, D506; sirloin: D701, D507; subcutaneous fat: D701, D508.17.65 Gb of short paired-end RNA reads were sequenced using Illumina NextSeq 500 (Table 1).
For Iso-Seq, a total of 600 ng RNA from sirloin was used for full-length transcript sequencing with Pacbio Sequel system (Pacific Biosciences, CA, USA) according to the manufacturer's instructions.The Iso-Seq library was prepared according to the Isoform Sequencing (Iso-Seq) protocol using the NEBNext Single Cell/Low Input cDNA Synthesis & Amplification Module, PacBio SMRTbell Express Template Prep Kit 2.0 and ProNex ® Size-Selective Purification System.
Total 10 μL library was prepared using PacBio SMRTbell Express Template Prep Kit 2.0.SMRTbell templates were annealed using Sequel Binding and Internal Ctrl Kit 3.0.The Sequel Sequencing Kit 3.0 and SMRT cells 1 M v3 LR Tray was used for sequencing.SMRT cells (Pacific Biosciences) using 1200 min movies were captured for each SMRT cell using the PacBio Sequel System (Pacific Biosciences).
Genome size estimation and contig assembly.Hanwoo contigs were assembled using the HiFi consensus reads and validated following the VGP (Vertebrate Genomes Project) assembly pipeline 2 .Adapter sequences of HiFi reads (5′-ATCTCTCTCTTTTCCTCCTCCTCCGTTGTTGTTGTTGAGAGAGAT-3′) were removed by Cutadapt (v 4.0) 3 .Counting k-mer and generating histogram of the k-mer count were performed on adapter trimmed sequences with k = 21 by Meryl (v 1.3.0) 4 .Genome properties such as genome size, maximum read depth and transition parameter were inferred using GenomeScope (v 2.0) 5 from the 21-mer histogram generated by Meryl (v 1.3.0) 4 .Genome size of Hanwoo was estimated as 3.06 Gb based on the k-mer histogram (Fig. 2).Trimmed reads were assembled to contig level using Hifiasm (v 0.16.1) 6 , and the draft contig assembly consisted of 1311 contigs totaling 3.28 Gb with an N50 of 55.23 Mb (Table 2).Haplotypic duplication and low-coverage contigs of the draft contig assembly were removed using Purge_dups (v 1.2.5) 7 after self-alignment using Minimap2 8 .The primary contig assembly after removing haplotypic duplication included 603 contigs, with a size of 3.11 Gb and a contig N50 of 58.14 Mb.
Scaffolding and gap filling.The Hanwoo contigs after removing haplotypic duplication were scaffolded on autosome of ARS-UCD1.3,through reference-guided approach by RagTag (v 2.1.0) 9.Because the Y chromosome is absent in ARS-UCD1.3,autosome and X chromosome of ARS-UCD1.3,and Y chromosome of UOA_Angus_1 were used as reference genome for scaffolding.The reference-guided scaffolding using RagTag (v 2.1.0) 9consist of 'correct' and 'scaffold' steps.The 'correct' step identified and corrected potential misassembly based on alignment of contig assembly to the reference genome assembly.Part of contigs were broken at points of putative misassembly, and as a result, the number of contigs increased to 1915.In the 'scaffold' step, these RagTag 'corrected' contigs were aligned to the reference genome consist of autosome and X chromosome of ARS-UCD1.3,and Y chromosome of UOA_Angus_1.As a result, there were 1598 scaffolds including 31 chromosome-level scaffolds and 1567 unplaced scaffolds.
HiFi reads used in the Hanwoo assembly were aligned using Minimap2 8 to perform gap filling of the chromosome-level Hanwoo genome assembly using TGS-GapCloser (v 1.0.1) 10 .The final 31 chromosome-level scaffolds had a total size of 2.69 Gb, which was similar to chromosome size of ARS-UCD 1.3.(Tables 3, 4).These 31 chromosome-level scaffolds composed 86.66% of the assembly, with the remaining 414.6 Mb still unanchored and requiring further investigation.Further analysis including annotation and pangenome were performed on the chromosome-level scaffolds.
Circos plot denoting gene density, N ratio and GC content was generated with the advanced circos function from Java-based tool TBtools 11 .The gene density (number of genes), N ratio (%) and GC content (%) was calculated for every 10,000 bp increment of the genome and was visualized in a heatmap format for gene density and histogram format for N ratio and GC content using BIN size 100,000.

Genome annotation.
Illumina RNA-seq reads were trimmed to remove adapter sequences and low-quality bases using Trimmomatic (v 0.39) 13 .The BRAKER3 (v 3.0.3)pipeline was used for structural annotation of Hanwoo genome.The pipeline utilized three sources of extrinsic evidence; short-read RNA-seq (Illumina), protein sequences of Vertebrata in OrthoDB (v 11) 14 in addition to protein sequence of ARS-UCD1.3 to train Augustus (v 3.5.0) 15for gene prediction.The predicted gene sets were searched in 2 public functional databases, Swiss-Prot of UniProtKB 16 and Pfam (v 35.0) database 17 to identify the potential function with BLASTP (v 2.13.0+) 18 and functional domains with InterProScan (v 5.57) 19 .We used scripts included in MAKER (v 3.01.03) 20to integrate functional annotations into structural annotations.The protein annotation was evaluated by analyzing amino acid sequences of protein using BUSCO (v 5.3.2) 21 with the conserved core set of mammalian genes, yielding a completeness score of 87.9%.A total of 27,314 protein-coding genes were identified, including 25,302 genes with inferred names and 2,012 unknown proteins.assessment of the chromosome-level genome assembly.N50, L50 and lengths of the chromosome-level Hanwoo genome assembly was calculated by QUAST (v 5.0.2) 22 .Single copy gene completeness was assessed with BUSCO (v 5.3.2) 21 , using the metaeuk backend against 'mammalia_odb10' .Quality values (QV) was calculated with Merqury (v 1.3) 23 , with k-mer databases (k = 21) constructed by Meryl (v 1.3) 4 .
Non-reference nodes in pangenome graph.The multiple whole-genome alignments generated by CACTUS 24 were transformed into the Packed Graph (PG) format by chopping into 32 base pairs using 'hal2vg' with the options '-chop 32' and '-noAncestors' 32 .The reference nodes and non-reference nodes were separated using scripts from the Github repository (https://github.com/evotools/CattleGraphGenomePaper/tree/master/detectSequences/nf-GraphSeq) following previous research 29 .After excluding nodes flanking with gaps in 1 kb, the counts and lengths of the non-reference and breed-specific nodes were calculated (Table 6).Non-reference region and Hanwoo-specific regions longer and equal to 10 kb are marked in Hanwoo autosome using KaryoploteR 33 (Fig. 3).The Hanwoo-specific regions are encompassed within the non-reference region, with the majority of these regions being located in the telomeric and centromeric regions.Notably, the size of satellite repeats, as identified by RepeatMasker, amounted to 39.4 Mb (Table 5).The total size of the satellite repeat, a main component of the centromere, were similar to the differences in autosome length between Hanwoo and others.This finding implies that the larger genome and specific region of Hanwoo can be attributed to expansions within repeat-rich telomeric and centromeric regions.
Furthermore, HiFi-based assemblies generally have higher telomeric completeness than Oxford nanoporeor CLR-based assemblies 34 .The uniqueness of origin and evolution history also supported the larger and disctinct genome of Hanwoo compared to European taurine.Mitochondrial DNA haplogroup of Hanwoo is P, which is common in European aurochs but has not been detected in modern cattle in Europe 35 .The haplogroup P mtDNA in Hanwoo suggested the possibility of a minor and local event of domestication or introgression of Asian aurochs 36,37 .Furthermore, intensive inbreeding and small effective population size of Hanwoo might facilitate fixation of these distinctive regions in Hanwoo genome 38 .

Data Records
The final genome assembly was deposited at DDBJ/ENA/GenBank under the accession JARDUZ000000000 39 .This Whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession SRR23238456 40 .
The transcriptomic Illumina sequencing data of subcutaneous fat, liver and sirloin were deposited in the SRA at NCBI SRR23238453, SRR23238454 and SRR23238455, respectively 40 .
The transcriptomic PacBio sequencing data of sirloin were deposited in the SRA at NCBI SRR23238452 40 .
The Hanwoo genome assembly which were not processed by NCBI, genome annotation, transcript sequence and protein sequence are available in figshare 41 .
The pangenome graph in GFA format are also available in figshare 42 .

technical Validation
RNA degradation and contamination were monitored on Agilent RNA ScreenTape.The purity of RNA samples was checked using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA).The completeness of the Hanwoo genome assembly was evaluated using BUSCO 21 with the mammalian data set "mammalia_odb10." The evaluation found 95.8% (8842) of the core mammalian genes were present in the genome, including 93.9% single-copy, 1.9% duplicated, 1.9% fragmental, and 3.1% missing genes from the mammalian data set (Table 3).The k-mer databases (k = 21) constructed using HiFi reads by Meryl 4 , and the overall assembly quality was assessed using the k-mer databases using Merqury 23 .The assembly showed high quality values (QV > 63) with an error rate of 4.29 × 10 −7 (Table 3).The GC content of Hanwoo (43.44%) was slightly higher than that of ARS-UCD1.3(41.56%).These assessment results confirmed the completeness of Hanwoo genome assembly (Table 3).

Fig. 3
Fig.3Non-reference region and specific region in Hanwoo autosome.Non-reference regions and Hanwoospecific regions larger than or equal to 10 kb are visualized on Hanwoo autosomes.The Hanwoo-specific regions are marked in red, while the non-reference regions shared by other Bos taurus assemblies, excluding the Hanwoo-specific regions, are marked in blue.

Table 1 .
Statistics of sequencing data.

Table 2 .
Statistics of contig assembly before scaffolding.

Table 4 .
Length of Chromosome-level scaffolds.

Table 5 .
Statistics of repetitive elements.
The integrity of RNA was assessed using the RNA ScreenTape of the Agilent 2200 TapeStation System (Agilent Technologies, CA, USA).Only RNAs with an OD260/280 ratio of 2.0-2.2, an OD260/230 ratio of 1.8-2.1, and a RIN value of ≥9.0 were considered qualified for use.RNA concentration was measured using Quant-iT ™ RiboGreen ™ RNA Assay Kit in Victor Nivo (PerkinElmer, Waltham, MA, USA).