# Fourth report on chicken genes and chromosomes 2022

Jacqueline Smith, James Alfieri, Nick Anthony, Peter Arensburger, Giridhar Athrey, Jennifer Balacco, Adam Balic, Philippe Bardou, Paul Barela, Yves Bigot, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jacqueline Smith, James Alfieri, Nick Anthony, Peter Arensburger, Giridhar Athrey, et al. Fourth report on chicken genes and chromosomes 2022. Cytogenetic and Genome Research, 2023, 162 (8-9), pp.405-528. 10.1159/000529376. hal-04217636

HAL Id: hal-04217636

https://hal.inrae.fr/hal-04217636

Submitted on 26 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Fourth Report on Chicken Genes and Chromosomes 2022

Prepared by

Jacqueline Smithª James M. Alfieriᵇ, c, d Nick Anthonyᵉ Peter Arensburgerf Giridhar N. Athreyb, d

Jennifer Balaccog Adam Balicª Philippe Bardouh Paul Barelai Yves Bigotj Heath Blackmonb, c

Pavel M. Borodink Rachel Carrolll Meya C. Casonoi Mathieu Charlesm Hans Chengn Maddie Chiodi

Lacey Cigani Lyndon M. Coghillo Richard Crooijmansp Neelabja Dasi Sean Daveyi Asya Davidianq

Fabien Degalezr Jack M. Dekkerss, t Martijn DerksP Abigail B. Diackª Appolinaire Djikengu

Yvonne Drechslerv Alexander Dyominq Olivier Fedrigog Steven R. Fiddamanw Giulio Formentig

Laurent A.F. Frantzₓ, y Janet E. Fultonz Elena Gaginskayaq Svetlana Galkinaq Rodrigo A. Gallardo s, A

Johannes Geibelb, C Almas A. Gheyasª Cyrill John P. GodinezD Ashton Goodelli Jennifer A.M. GravesE, F

Darren K. GriffinG Bettina Haaseg Jian-Lin HanH, I Olivier Hanottel, J, K Lindsay J. Hendersona

Zhuo-Cheng HouL Kerstin HoweM Lan HuynhN Evans IlatsiaO Erich D. Jarvisg Sarah M. Johnsoni

Jim KaufmanN, P, Q Terra Kellys, A Steve KempR Colin KernS Jacob H. Keroacki Christophe KloppT

Sandrine Lagarriguer Susan J. Lamonts, t Margaret LangeU Anika LankeV Denis M. LarkinW

Greger LarsonX John King N. LayosY Ophélie LebrasseurZ, α Lyubov P. Malinovskayaβ

Rebecca J. MartinQ Maria Luisa Martin CerezoV Andrew S. Masonδ Fiona M. McCarthyi

Michael J. McGrewa, u Jacquelyn Mountcastleg Christine Kamidi MuhonjaO, R William Muirε Kévin Muretz

Terence D. Murphyη Ismael Ng'ang'ax Masahide Nishiboriθ Rebecca E. O'ConnorG Moses Ogugor

Ron Okimotoe Ochieng OukoO Hardip R. PatelI Francesco Periniª, K María Ines Pigozziλ Krista C. Potteri

Peter D. Priceμ Christian ReimerB Edward S. Ricev Nicolas RocosN Thea F. Rogersζ Perot Saelao s, S, o

Jens SchauerB Robert D. Schnabelπ Valerie A. Schneiderη Henner SimianerC Adrian Smithw

Mark P. Stevensª Kyle Stierso Christian Keambou TiamboR Michele Tixier-Boichardm

Anna A. TorgashevaK Alan TraceyM Clive A. TregaskesP, Q Lonneke Verveldeª Ying Wangss, S

Wesley C. Warrenv, π Paul D. WatersP David Webbη Steffen WeigendB, C Anna Wolct, z Alison E. Wrightμ

Dominic WrightY Zhou Wuª Masahito Yamagatao Chentao YangT Zhong-Tao YinL Michelle C. Youngi

Guojie Zhangu Bingru Zhaoφ Huaijun Zhou s, S

ªThe Roslin Institute and Royal (Dick) School of Veterinary Studies, University of Edinburgh, Easter Bush Campus, Edinburgh, UK;

ᵇInterdisciplinary Program in Ecology and Evolutionary Biology, Texas A&M University, College Station, TX, USA; ᵇDepartment of

Biology, Texas A&M University, College Station, TX, USA; ᶜDepartment of Poultry Science, Texas A&M University, College Station, TX,

USA; ᶜCobb-Vantress, Siloam Springs, AR, USA; ᵈBiological Sciences Department, California State Polytechnic University, Pomona,

CA, USA; ᵉThe Rockefeller University, New York, NY, USA; ᵇUniversité de Toulouse, INRAE, ENVT, GenPhySE, Sigenae, Castanet

Tolosan, France; ᵈUniversity of Arizona, Tucson, AZ, USA; ᵈPRC, UMR INRAE 0085, CNRS 7247, Centre INRAE Val de Loire, Nouzilly,

France; ᵇDepartment of Molecular Genetics, Cell Biology and Bioinformatics, Institute of Cytology and Genetics of Siberian Branch

of Russian Academy of Sciences, Novosibirsk, Russia; ᵈDepartment of Animal Sciences, Data Science and Informatics Institute,

University of Missouri, Columbia, MO, USA; ᵈUniversity Paris-Saclay, INRAE, AgroParisTech, GABI, Sigenae, Jouy-en-Josas, France;
ªUSDA, ARS, USNPRC, Avian Disease and Oncology Laboratory, East Lansing, MI, USA; ᵇDepartment of Veterinary Pathology, University of Missouri, Columbia, MO, USA; ᶜAnimal Breeding and Genomics, Wageningen University and Research, Wageningen, The Netherlands; ᵈSaint Petersburg State University, Saint Petersburg, Russia; ᵉINRAE, INSTITUT AGRO, PEGASE UMR 1348, Saint-Gilles, France; ᵈFeed the Future Innovation Lab for Genomics to Improve Poultry, University of California, Davis, CA, USA; ᵒDepartment of Animal Science, Iowa State University, Ames, IA, USA; ᵒCentre for Tropical Livestock Genetics and Health (CTLGH) – The Roslin Institute, Edinburgh, UK; ᵒCollege of Veterinary Medicine, Western University of Health Sciences, Pomona, CA, USA; ᵒDepartment of Zoology, University of Oxford, Oxford, UK; ᵒQueen Mary University of London, Bethnal Green, London, UK; ᵒPalaeogenomics Group, Department of Veterinary Sciences, LMU Munich, Munich, Germany; ᵒHy-Line International, Research and Development, Dallas Center, IA, USA; ᵒSchool of Veterinary Medicine, University of California, Davis, CA, USA; ᵒInstitute of Farm Animal Genetics, Friedrich-Loeffler-Institut, Neustadt, Germany; ᵒCenter for Integrated Breeding Research, University of Göttingen, Göttingen, Germany; ᵒDepartment of Animal Science, College of Agriculture and Food Science, Visayas State University, Baybay City, Philippines; ᵒDepartment of Environment and Genetics, La Trobe University, Melbourne, VIC, Australia; ᵒInstitute for Applied Ecology, University of Canberra, Canberra, ACT, Australia; ᵒSchool of Biosciences, University of Kent, Canterbury, UK; ᵒCAAS-ILRI Joint Laboratory on Livestock and Forage Genetic Resources, Institute of Animal Science, Chinese Academy of Agricultural Sciences (CAAS), Beijing, China; ᵒInternational Livestock Research Institute (ILRI), Addis Ababa, Ethiopia; ᵒCells, Organisms and Molecular Genetics, School of Life Sciences, University of Nottingham, Nottingham, UK; ᵒCentre for Tropical Livestock Genetics and Health, The Roslin Institute, Edinburgh, UK; ᵒNational Engineering Laboratory for Animal Breeding and Key Laboratory of Animal Genetics, Breeding and Reproduction, MARA, College of Animal Science and Technology, China Agricultural University, Beijing, China; ᵒWellcome Trust Sanger Institute, Hinxton, UK; ᵒInstitute for Immunology and Infection Research, University of Edinburgh, Edinburgh, UK; ᵒDairy Research Institute, Kenya Agricultural and Livestock Organization, Naivasha, Kenya; ᵒDepartment of Veterinary Medicine, University of Cambridge, Cambridge, UK; ᵒDepartment of Pathology, University of Cambridge, Cambridge, UK; ᵒCentre for Tropical Livestock Genetics and Health (CTLGH) – ILRI, Nairobi, Kenya; ᵒDepartment of Animal Science, University of California, Davis, CA, USA; ᵒINRAE, MIAT UR875, Sigenae, Castanet Tolosan, France; ᵒDepartment of Molecular Microbiology and Immunology, University of Missouri, Columbia, MO, USA; ᵒBASIS Chandler High School, Chandler, AZ, USA; ᵒDepartment of Comparative Biomedical Sciences, Royal Veterinary College, University of London, London, UK; ᵒThe Palaeogenomics and Bio-Archaeology Research Network, Research Laboratory for Archaeology and History of Art, The University of Oxford, Oxford, UK; ᵒCollege of Agriculture and Forestry, Capiz State University, Mambusao, Philippines; ᵒCentre d'Anthropobiologie et de Génomique de Toulouse (CAGT), CNRS UMR 5288, Université Toulouse III Paul Sabatier, Toulouse, France; ᵒInstituto Nacional de Antropología y Pensamiento Latinoamericano, Ciudad Autónoma de Buenos Aires, Argentina; ᵒDepartment of Cytology and Genetics, Novosibirsk State University, Novosibirsk, Russia; ᵒAVIAN Behavioural Genomics and Physiology, IFM Biology, Linköping University, Linköping, Sweden; ᵒDepartment of Biology, The University of York, York, UK; ᵒDepartment of Animal Sciences, Purdue University, West Lafayette, IN, USA; ᵒUniversité Paris-Saclay, Commissariat à l'Energie Atomique et aux Energies Alternatives, Centre National de Recherche en Génomique Humaine, Evry, France; ᵒNational Center for Biotechnology Information, National Library of Medicine, National Institutes of Health, Bethesda, MD, USA; ᵒLaboratory of Animal Genetics, Graduate School of Integrated Sciences for Life, Hiroshima University, Higashi-Hiroshima, Japan; ᵒThe John Curtin School of Medical Research, Australian National University, Canberra, ACT, Australia; ᵒDepartment of Agricultural, Food and Environmental Sciences, University of Perugia, Perugia, Italy; ᵒINBIOMED (CONICET-UBA), Facultad de Medicina, Universidad de Buenos Aires, Buenos Aires, Argentina; ᵒEcology and Evolutionary Biology, School of Biosciences, University of Sheffield, Sheffield, UK; ᵒDepartment of Animal Sciences, Bond Life Sciences Center, University of Missouri, Columbia, MO, USA; ᵒDepartment of Molecular Evolution and Development, University of Vienna, Vienna, Austria; ᵒVeterinary Pest Genetics Research Unit, USDA, Kerrville, TX, USA; ᵒDepartment of Animal Sciences, University of Missouri, Columbia, MO, USA; ᵒSchool of Biotechnology and Biomolecular Science, Faculty of Science, UNSW Sydney, Sydney, NSW, Australia; ᵒCenter for Brain Science, Department of Molecular and Cellular Biology, Harvard University, Cambridge, MA, USA; ᵒBGI-Shenzhen, Shenzhen, China; ᵒCenter for Evolutionary and Organismal Biology, Zhejiang University School of Medicine, Hangzhou, China; ᵒCollege of Animal Science and Technology, Nanjing Agricultural University, Nanjing, China

## Introduction

(Prepared by J. Smith)

The chicken continues to hold its position as a leading model organism within many areas of research, as well as being a major source of protein for human consumption. The First Report on Chicken Genes and Chromosomes [Schmid et al., 2000], which was published in 2000, was the brainchild of the late, and sadly missed, Prof. Michael Schmid of the University of Würzburg. It was a publica-

tion bringing together updates on the latest research and resources in chicken genomics and cytogenetics. The success of this first report led to the subsequent publication of the Second [Schmid et al., 2005] and Third Report on Chicken Genes and Chromosomes [Schmid et al., 2015] – each also proving popular references for the research community. It is now our pleasure to be able to introduce publication of the Fourth Report. Being 7 years since the last report, this publication captures the many advances that have taken place during that time. This includes pre-
sentation of the detailed genomic resources that are now
available, largely due to increasing capabilities of se-
quencing technologies and which herald the pangenomic
age, allowing for a much richer and more complete
knowledge of the avian genome. Ongoing cytogenetic
work also allows for examination of chromosomes, spe-
cific elements within chromosomes, and the evolutionary
history and comparison of karyotypes. We also examine
chicken research efforts with a much more “global” out-
look with a greater impact on food security and the im-
pact of climate change, and highlight the efforts of inter-
national consortia, such as the Chicken Diversity Consor-
tium. We dedicate this Report to Michael.

Multiple Chicken (Gallus gallus) Genome References
to Advance Genetic Variation Studies

(Prepared by W.C. Warren, O. Fedrigo, A. Tracey, A.S. Mason, G. Formenti, F. Perini, Z. Wu, T.D. Murphy, V.A. Schneider, K. Stiers, E.S. Rice, L.M. Coghill, N. Anthony, R. Okimoto, R. Carroll, J. Mountcastle, J. Bacacco, B. Haase, C. Yang, G. Zhang, J. Smith, Y. Drechsler, H. Cheng, K. Howe, and E.D. Jarvis)

We present two phased chromosome-scale assemblies
of chicken, a layer (GRCg7w) and broiler (GRCg7b), that
better meet research demands to characterize segregating
variation important for traits of interest. Annotation with
existing long- and short-read RNAseq data improved con-
tiguity, accuracy, and protein-coding and noncoding gene
counts, when compared to the existing Red Jungle Fowl
reference, GRCg6a. Most striking were the improvements
in placed telomeres, corrections for erroneous microchro-
mosome fusions, and gap reduction in these phased as-
semblies. We add 6 putative microchromosomes that
were previously missing in GRCg6a. Using a pairwise ge-
nome comparison of the parental genomes, and 2 inde-
pendent cohorts of sequenced chickens, we show small
discernible differences in mapping rates of whole genome
sequence (WGS) and RNAseq data, gene annotation, and
called single nucleotide variants (SNVs) or indels. Struc-
turally, some regional differences suggest future assembly
curation will further improve variant ascertainment.
These Gallus references also enabled a new genome-wide
review of endogenous Avian Leukosis Virus (ALVE) inte-
grations, exemplifying the improved representation of
chicken genomic diversity by these phased genomes. Our
genome references will collectively improve computa-
tional outcomes when testing multiple variant hypotheses
that are at the core of understanding avian biology.

Today, the poultry industry faces many challenges,
perhaps none more than the genetics underlying bird
health. A constant balance must be maintained to select
on several traits of immense economic impact, such as
fast growth in broilers and reproductive success in layers,
while not diminishing disease resistance. Genetic studies
offer promising avenues to selectively maintain this trait
balance with a new accounting of the most important
contributing factors: genes and environment [Wolc et al.,
2018]. More complete and accurate genomic resources to
support their continued discovery of these factors is par-
amount to generating the robust chicken germlines that
can meet a growing demand for this protein food source.

The chicken also supports a vast model organism com-
munity that uses a collateral source of scientific data to
comparatively inform developmental biology [see review
Cheng and Burt, 2018]. The chicken genome is also one
of the most frequently used resources for comparative ge-
nomic studies among vertebrates [Zhang G et al., 2014].
As the principal avian reference genome, it was used to
transfer gene annotation evidence to over 50 bird ge-
nomes, which were in turn used for clade- and species-
specific signals of genome evolution [Jarvis et al., 2014;
Zhang G et al., 2014]. Recent research attempts to deter-
mine the effect of structural variation (SV) on chicken
phenotypic differences, although resolution beyond
short-read mapping or hybridization methods must be
considered [Rao et al., 2016]. The site-directed gene
knockouts of chicken C2EIP [Zuo et al., 2016] and PAX
[Gandhi et al., 2017] genes are 2 functional instances of
gained insight into embryonic germ and satellite muscle
cell differentiation, respectively.

Since its first iteration in 2004 [International Chicken
Genome Sequencing Consortium, 2004], we have worked
to refine the assembly of the chicken genome as technol-
ogy progressed over the past 2 decades [Korlach et al.,
2017]. We recently summarized these advances [Rhie et
al., 2021], like the single haploid phasing of a diploid ge-
nome, i.e., trio binning [Koren et al., 2018], which sorts
and independently assembles divergent parental haplo-
types from F1 hybrids (inter- and intraspecies crosses) as
a highly efficient method for untangling complex se-
quence assembly graphs. This phasing strategy exploits
the higher heterozygosity observed in some F1 hybrids in
order to resolve diploid genomes more precisely and with
fewer gaps. Several successful haplotype-resolved de novo
assemblies for cats, cattle, zebra finches, and others have
been created using this method [Koren et al., 2018; Rice
et al., 2020; Bredemeyer et al., 2021; Rhie et al., 2021],
highlighting the astounding improvements in contiguity,
<table><caption>**Table 1.** Phased assembly comparisons of broiler and layer genomes to RJF for *Gallus gallus*. Each assembly contains the Z and W sex chromosomes despite there being only one copy of each from the parents</caption><thead><tr><th>Common name</th><th>Assembly version</th><th>N50 contig, Mb</th><th>Total size, Mb</th><th>Total contigs</th><th>N50 scaffold length, Mb</th><th>Unplaced sequences, Mb</th></tr></thead><tbody><tr><td>Red Jungle Fowl</td><td>GRCg6a</td><td>17.6</td><td>1,055</td><td>1,403</td><td>20</td><td>14.1</td></tr><tr><td>Broiler</td><td>GRCg7b</td><td>18.8</td><td>1,049</td><td>677</td><td>90</td><td>6.6</td></tr><tr><td>Layer</td><td>GRCg7w</td><td>17.7</td><td>1,046</td><td>685</td><td>90</td><td>7.1</td></tr></tbody><tfoot><tr><td colspan="7">Data based on NCBI assembly metrics.</td></tr></tfoot></table>

with some sequences spanning from telomere to telomere.

Recent characterization of different chicken genomes has demonstrated the necessity for pangenome resources in order to comprehend the comparative evolution of the *Gallus* genus [Li M et al., 2022]. To date, all chicken genetic studies have relied on the Red Jungle Fowl (RJF) genome, which portrays this diploid genome as a collapsed haploid genome containing a mixture of sequences from the two haplotypes. For the further investigation of trait selection indices, the adoption of additional high-quality chicken references, particularly those that better resemble commercial birds, is highly endorsed by the avian community and has broad applicability. In addition, these resources enable pangenome techniques that will provide higher resolution for discovering SVs that are exclusive to decades of artificial selection. Here, we used the trio binning approach to present 2 novel haploid de novo assemblies for chicken lines with extremely diverse genetic histories: one bred for muscle growth (broiler) and the other for egg production (layer). Large structural adjustments among microchromosomes, overall gap reduction, extension of the W chromosome by adding the pseudoautosomal region (PAR), better chromosome 16 (MHC region) representation, and enhanced telomere sequence placements are notable. We demonstrate these new assemblies' application in determining the extent of alignment, SNV identification, ALVE integration, and structural expansions and contractions in a small sample of chickens.

## Sequencing and Assembly

A parent-offspring trio composed of a paternal layer, a maternal broiler, and a female F1 offspring was sequenced to create these assemblies. Briefly, the parents were sequenced with low-coverage Illumina reads (150 bp), and the F1 was sequenced with 80× PacBio reads (12 kb on average), and all reads were used as input to TrioCa-

nu [see for review of methods, Rhie et al., 2021]. Similar to cross-species trio assembly of cattle and yak, the amount of haplotyped long reads phased from each paternal breed source was extremely similar (49.5 and 50.2%) with a low number of unknowns (0.16%) [Rice et al., 2020]. Broiler (n = 676) and layer (n = 688) birds had half as many constructed contigs as RJF (n = 1,403) birds, indicating that >53% of prior gaps have been bridged (currently 878 in RJF). Manuel curation with orthogonal evidence, including chromatin proximity (Hi-C) and Bisonano optical maps, delineated error locations that were fixed, e.g. 260 and 63 missed joins in GRCg7b and GRCg7w (online suppl. Material 1, Table 1; for all online suppl. material, see www.karger.com/doi/10.1159/000529376). Depending on the descriptive context, we use the assembled GenBank versions (GRCg6a, GRCg7b, and GRCg7w) and their common names (RJF, broiler, and layer) interchangeably throughout the remainder of this report.

While contig N50 length was comparable to GRCg6a (which also used PacBio long read data to fill gaps), phasing and revised mapping data led to a 4.5-fold increase in N50 scaffold length and a 2-fold decrease in the number of unplaced sequences in broiler and layer assemblies (Table 1). The paternal layer contributes Z to the ZW sex chromosomes, while the maternal broiler was particularly chosen for her haplogroup A mitochondrial genome and W. The female RJF reference is unique with a mitochondrial genome of haplogroup E. The layer Z chromosome is somewhat larger and contains more protein-coding genes than the BAC-curated GRCg6a version of Z (85.2 vs. 80 Mb; 1,492 vs. 1,345 genes) [Bellott et al., 2010]. Furthermore, the broiler W chromosome is more complete than the GRCg6a chromosome, which is 7.2 Mb in size, due in part to the insertion of the PAR that boosts its comparative utility (online suppl. Material 2, Fig. 1). The initial about 500 kb of the W chromosome show diploid coverage. We chose not to join this sequence to the begin-
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Assembled structural errors detected in RJF compared to broiler for chromosome 27 using Hi-C mapped data to the scaffolds. Genetic linkage map markers (n = 125) displayed as green tic marks below the x axis for the chromosome 27 heat map were mapped to each assembly to validate sequence order and orientation.</figcaption></figure>

ning of Z since we lack precise coordinates, and this portion of Z is partially collapsed. In each phased reference, for the sake of completeness, Z and W were incorporated notwithstanding their parental origins. By searching the NCBI assembly archive for “Gallus gallus,” you can find all fully annotated assemblies (see data availability).

### Assembly Accuracy Benchmarking

There are inherent assembly artifacts present in all reference genomes, including the human genome. With this knowledge, we wanted to estimate the detected errors in GRCg6a, given its extensive use in chicken genetic studies, and repair them in our new phased assemblies using a previously established iterative procedure [Howe et al., 2021]. A greater number of GRCg7b and GRCg7w chromosomes exhibited telomere ends (24 and 13, respectively), than GRCg6a (just 3), demonstrating the much-improved completeness of the new assemblies. Among mi-

crochromosomes, we discovered many instances in which GRCg6a chromosomes were wrongly fused into a single chromosome instead of two distinct ones (online suppl. Material 1, Table 1). The first 2 Mb of GRCg6a chr27 is not associated with chr27, but rather a variety of alignments to other chromosomes, including W and chr2, which, upon curation, accurately sizes this chromosome; 8 Mb as opposed to 5.2 Mb in GRCg7b (Fig. 1). Other errors include the fusion of chr31 and chr29 in GRCg6a, which is likely due to repeat sequences identified on the Hi-C heat map (online suppl. Material 2, Fig. 2).

Avian microchromosomes show more frequent recombination, and thus the positive correlation between recombination and interspecies divergence observed in mammals is not seen in birds, at least at the resolution of whole chromosomes [International Chicken Genome Sequencing Consortium, 2004]. If the origin of the microchromosomes was initiated by a number of random fis-
sion events that were channeled towards the present day macro/microchromosome arrangement, it was clear from earlier evaluations that there were not sufficiently long compositionally uniform regions in any of the sequenced avian genomes, that would satisfy some classifications, e.g., the classical isochore definition within microchromosomes [Waters et al., 2021]. We now have corrected several assembly errors, mostly among the microchromosomes, to test more accurately these and other hypotheses regarding their evolution.

The chicken karyotype has a diploid number of 78 chromosomes, classified as a haploid autosome count of 10 macrochromosomes and 28 microchromosomes [Burt et al., 1999]. In earlier chicken assemblies, microchromosomes 29 and 34–38 were absent, primarily due to the absence of linkage groups or physical maps that might assign missing scaffolds to any of these smaller chromosomes [Groenen et al., 2000], as well as difficulty in sequencing through high-GC rich microchromosomes. In both GRCg7b and GRCg7w, using long reads that get through GC-rich regions, as with zebra finch [Kim et al., 2022], we identify these missing microchromosomes (online suppl. Material 2, Fig. 3) and an additional microchromosome to a final total of 39 autosomes. Future cytogenetic evaluations or new combinatorial approaches that can yield telomere-to-telomere stepwise assembly of more complete chromosomes [Logsdon et al., 2021] will be necessary to rule out the possibility these nominated microchromosomes are not affiliated with other macro- or microchromosomes. Moreover, the availability of almost complete genome copies of these uniquely selected lines and others will drive reevaluations of all types of segregating variation in a pangenome-dependent manner [Siren et al., 2021].

## Structural Differences

To estimate the major structural differences among these phased references, we employed 2 methods: high resolution alignments to reveal major synteny differences using SyRi [Goel et al., 2019] and the predicted contractions and expansions of deletions, insertions, and repeat elements with different size distributions using Assemblytics [Nattestad and Schatz, 2016]. Across the chicken genome, differences in local chromosomal synteny were predominately one-to-one; however, when we discover discrepancies, they frequently occur towards chromosome ends, highlighting the difficult nature of placing sequences in these repetitive telomeric regions (Fig. 2). Regardless of their length distribution, Assemblytics alignment results show comparable total base size differences

(Fig. 2; online suppl. Material 1, Table 2, online suppl. Material 2, Fig. 4). However, these differences vary by type, such as deletion versus insertion, which may be the result of numerous factors, including genetic diversity and assembly completeness and accuracy of each reference. When employing a phased assembly for pairwise broiler versus layer alignments, the total number and base sizes of discovered deletions and insertions drop relative to RJF, suggesting the more diversified origins of RJF and its mixed haplotype assembly architecture are the cause (online suppl. Material 1, Table 2). A genome-wide perspective of SVs in RJF, layer, and broiler genomes, including overall alterations in repeat content, will require additional research to validate their patterns of segregation in larger populations of chickens and their accuracy of ascertainment. Overall, we observe structural differences, despite the fact that the percent masked sequence, a measure of all repeat types, is comparable between these references (20.5, 20.2, and 20.3) using default WindowMasker output [Morgulis et al., 2006].

## Gene Annotation

First, protein-coding gene representation was evaluated with BUSCO, which demonstrated an average 54% reduction in the number of missing universal single-copy orthologs in both GRCg7 haplotypes compared to GRCg6a (online suppl. Material 1, Table 3). Automated gene annotation of GRCg7b and GRCg7w using the NCBI workflow [Sayers et al., 2021] reveals an increase in the overall number of protein-coding and noncoding genes (online suppl. Material 1, Table 4). Recent gene annotation of multiple chicken genomes revealed 1,335 more protein-coding genes relative to GRCg6a [Li M et al., 2022]. The addition of at most 546 genes to the GRCg7 phased assemblies is a modest increase, but not unexpected given the NCBI annotation methods are likely more conservative and do not rely on many varied genome annotations as in Li M et al. [2022]. We also analyzed the differences in gene set ontology between GRCg7b and GRCg7w given their distinct selection histories. As determined by enrichment analysis, we find 82.8% overlap between GRCg7b and GRCg7w, where uniqueness is most often in large gene families, e.g., immune genes (online suppl. Material 2, Fig. 5; online suppl. Material 1, Table 5). In the broiler annotated set, there are 154 genes not found in the layer set, but there were only 44 unique to the layer set (online suppl. Material 1, Table 5). This disparity may reflect the slightly higher contiguity of the broiler reference (Table 1). Future study will be required to provide a precise accounting of the genes that are
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Sequenced differences in the phased broiler and layer genomes for macro- (a) and micro-autosomes (b). From the inside out SNV density (red), window size of 500 kb, range of 0 to 2.5%, indels <50 bp (coral), 500 kb window size and 0–0.8%; large indels</figcaption></figure>

(blue) per Mb, range of 0 to 60; CNV count per Mb (green); highlighted inversions (black dashes); chicken karyotype (varied color); ideograms of GRCg7b and GRCg7w chromosomes (varied colors).

unique to breeds, commercial or research lines, and wild strains of *Gallus gallus*. In addition, solutions to the question of whether avian genes were genuinely lost during the ancient divergence of avian and mammalian lineages will begin with the availability of a full genome, as recently completed for a human [Nurk et al., 2022].

## WGS Mapping and SNV Analysis

It is probable that the choice of broiler, layer or RJF as a reference for alignment of various resequenced chicken populations could contribute to SNV ascertainment bias as has been shown in human [Schneider et al., 2017]. We examined the mapping rates of WGS short-read data of 6 genetically diverse chicken samples (online suppl. Material 1, Table 6): a male and female for layer and broiler chicken, as well as an Ethiopian indigenous chicken breed. For various mapping metrics, regardless of the reference, we find no large differences, indicating that for measures of genetic diversity, all 3 references have comparable initial abilities to call SNVs or indels (online suppl. Material 1, Table 7).

Despite similar WGS mapping rates across references, the optimal SNVs set for the experimental purpose intended is not certain. Next, we mapped WGS data from a separate cohort of solely broilers (n = 10) to all three genome assemblies and called SNVs using GATK version 4.2.0. Although we found differences in total SNVs, these were not large, suggesting that SNV detection will be comparable when beginning with any of the three assemblies (online suppl. Material 1, Table 8). However, regional variations may be encountered and must be addressed if certain loci are of great experimental relevance (Fig. 3).

## RNAseq Mapping

The mapping of RNAseq data to estimate transcriptome changes between samples for biological interpretation is a crucial reference usage. To address this application, we first analyzed a large number of diverse tissues where total percent mapping is available in the NCBI gene annotation report (https://www.ncbi.nlm.nih.gov/genome/annotation_euk/Gallus_gallus/106/) and found very little difference between haplotypes (GRCg7b and
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** The distribution of called heterozygous SNVs across chick- en macrochromosome 7 (a) and microchromosome 20 (b) in the three assemblies. Rainfall plots of heterozygous variants depict their location, and each unique color indicates a different type of base substitution. We only include variants that passed all filters and were heterozygous in either reference source.</figcaption></figure>

GRCg7w). However, since RNAseq alignment in this ap-
plication is optimized for verifying gene model predic-
tions, we also tested the STAR aligner, which is typically
considered best practice for bulk RNAseq studies [Dobin
et al., 2013]. Using a small number (n = 8) of RNAseq
samples from diverse tissue origins, including ileum,
bone-derived macrophages, and uterus, we observe a
small average range of 0.7–1.7% differences among 6
samples in the total percentage of reads uniquely map-
ping to each of the three references (online suppl. Mate-
rial 1, Table 9). However, it is unclear why GRCg6a has a
somewhat greater percentage of uniquely mapped reads
across all samples (online suppl. Material 1, Table 9). We
also highlight the two individual female and male layer
muscle RNAseq samples with vastly different average
rates of unique mapping between assemblies, 70.6 and
68.5% to each GRCg7 haplotype, respectively, compared
to 88% in GRCg6a (online suppl. Material 1, Table 9). In
the female layer sample, after analyzing all secondary
alignment counts (not read counts), the olfactory recep-
tor 14C36-like gene (>10 M) is most abundant in GRC-
g6a, whereas for GRCg7 haplotypes it is ribosomal RNAs
(Fig. 4). Although the number of individual rRNAs in
each reference is comparable (~90), the total base size is
notably different. The total assembled lengths of rRNA in
each genome are GRCg6a (20,568), GRCg7b (113,817),
and GRCg7w (55,916) (Fig. 4), suggesting that when eval-
uating unique versus multi-mapping events for any of
these references, the type of library sequenced, i.e., ribo-
some depleted or polyA selected should be considered

[Zhao S et al., 2018]. Overall, for the majority of RNAseq
samples studied, we observe minimal differences in
unique mapping rates across all references; nonetheless,
prior to conducting RNAseq mapping experiments, the
reference choice should be considered.

Contrasting ALVE Diversity in Varied Genome
Assemblies

To show an additional advantage of these phased ref-
erences, we revisit the question of the RJF reference be-
ing unrepresentative of ancestral *Gallus gallus* ALVE
diversity. ALVEs are species-specific retroviral integra-
tions which retain the potential for retrotransposition
and retroviral expression (Fig. 5a). The previous RJF
reference assembly contained 2 Avian Leukosis Virus
subgroup E (ALVE) integrations: ALVE6 (ALVE-JFe-
VA), a truncated ALVE widespread across many breeds;
and ALVE-JFevB, an intact integration found in no oth-
er chicken to date [Mason et al., 2020a]. The new GRCg7
phased assemblies contain a total of 11 ALVEs: 5 from
the maternal broiler and 6 from the paternal layer (sum-
marized in Fig. 5b; detailed locations in online suppl.
Material 1, Table 10). Leghorn layers typically have few-
er than 6 ALVEs [Mason et al., 2020b, c], but the identi-
fied ALVE1, ALVE3, ALVE15, ALVE_ros034 and the
slow feathering-associated ALVE21 are representative
of this Leghorn layer breed. ALVE_ros005, however,
has previously only been identified in brown-egg layers
and Ethiopian indigenous birds [Mason et al., 2020b].
The ALVEs of the broiler haplotype are widely found
across brown-egg commercial layers and broiler lines,
representing their recent shared ancestry [Muir et al.,
2008b], and the presence of ALVE-TYR supports the
observed recessive white phenotype [Fox and Smyth,
1985; Chang et al., 2006].

Seven of the ALVEs are full-length (Fig. 5b; online
suppl. Material 1, Table 10), and 5 have completely in-
tact retroviral ORFs, accounting for the –1 ribosomal
frameshift between gag and pol [Nikolic et al., 2012].
Despite this, ALVE transmission between cells is un-
likely, as both parental haplotypes exhibit ALVE-resis-
tance at the TVB receptor (TNFRSF10B): maternal
Q58* (TVB<sup>R</sup>; rs736008824) and paternal P61L
(rs318006572). This perhaps represents the effects of
selection against P27 expression in commercial birds.
Additionally, the ALVE_ros034 gag ORF truncates
within P27, and similar mutations have been observed
in ALVEs in other commercial backgrounds (Fig. 5b)
[Mason et al., 2020b]. Additional high quality chicken
genome references of diverse genetic backgrounds in-
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** The RNAseq alignment detection of multimapping events and rRNA number and size distributions by reference source.</figcaption></figure>

interpreted in pangenome visualization modes will continue to resolve the evolution of ALVEs and their role in trait presentation.

## Future Chicken References

Pangenomic starting points as opposed to single linear representations have been proposed in humans [Siren et al., 2021] to overcome reference bias in genotyping. In Siren et al. [2021], the utility of a human pangenome reference in variant ascertainment demonstrates that this is the optimal course of action for future chicken genetic studies, particularly structural analyses. As a result, we are generating the requisite read-types to follow this same de novo assembly process in building multiple telomere-to-telomere single haplotype reference sources. Using these individual linear genome graphs to construct pangenome references will ensure the availability of the next generation of computational resources for optimally estimating segregating variation for significant genotype-to-phenotype connections in poultry production. The phased assemblies of the broiler and layer genomes as well as the RJF reference provide new insights into their general structure. In addition, we believe a new era in the use of avian genome references has already begun due to the rapid development of methods to build full genome copies.

## Bird Husbandry

The parent-offspring trio of this study is composed of a male White Leghorn and female broiler, the parents, each raised at the University of Arkansas avian housing facilities. A female F1 offspring was chosen from this cross for sequencing. DNA for each parent and the F1 was extracted from white blood cells using standard practices for each intended use.

## Sequencing and Primary Genome Assembly

We followed the workflow established by the Vertebrate Genomes Project (VGP) to create the haplotype-phased chicken assembly [Rhie et al., 2021]. Libraries were sequenced on the PacBio Sequel II instrument with the sequencing kit 2.1 (#101-310-500) and 10 h movie time to a total of ~98 GB. Because sequence coverage is lowered when phasing a diploid genome, we targeted a high read coverage of ~80×, to attempt the accurate assembly of repetitive microchromosomal regions and ZW sex chromosomes. TrioCanu (v1.8+287) was used to bin Consensus Long Reads (PacBio) of the F1 female into maternal and paternal haplotypes using haplotype-specific 21-mer markers derived from the Illumina short reads of the mother and father. Following binning, TrioCanu independently generated contigs for each haplotype (haplotigs). From this point, the maternal and paternal haplotigs independently underwent the same
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** ALVE integration, propagation, and degradation within the chicken genome. **a** Retroviral genomic lifecycle. Retroviral positive sense, single-stranded RNA is reverse transcribed into cDNA and associates with the retroviral integrase integration complex, which primes the cDNA 3' ends and initiates strand transfer with genetic DNA. Integration creates overhangs which are repaired by host machinery, creating target site duplications (TSDs; grey). Following integration, retroviral expression and retrotransposition is possible. Over evolutionary timescales integrated ERVs degrade,</figcaption></figure>

either by nonhomologous recombination events (I, II) or internal LTR recombination leaving solo LTRs (III). **b** Schematic indicates an intact ALVE with putative transcripts, with the ribosomal –1 frame slip and recognition site for miR-155 indicated. Phased chicken genome ALVE content and integrity is shown, with likely transcript and regulatory implications. CA, capsid; INT, integrase; LTR, long terminal repeat; MA, matrix; NC, nucleocapsid; PR, protease; RH, RNaseH; RT, reverse transcriptase; SU, surface; TM, transmembrane.

steps. Separately, we assembled the mitochondrial (MT) genome with the mitoVGP pipeline (v2.2) [Formenti et al., 2020] and added it to the haplotigs to keep any raw MT reads from being mapped to nuclear sequences preventing conversion of possible mitochondrial nuclear integrations into MT sequence during the polishing steps. We used Arrow from smrtlink (v6.0.0.47841) to improve base calling accuracy and purge_dups (v1.0.0) [Guan et al., 2020] in an adapted trio mode to remove erroneous duplications.

The median insert sizes of WGS libraries were approximately 400 bp and individual libraries were tagged with unique dual index DNA barcodes to allow pooling and minimize the impact of barcode hopping. Libraries

were pooled for sequencing on the NovaSeq 6000 (Illumina) to obtain at least 750 million 151-bp reads per individual.

### Assembly Scaffolding and Curation

Various maps were constructed to facilitate scaffolding of the phased contigs [Rhie et al., 2021]. Briefly, long linked read libraries were generated from unfragmented high molecular weight DNA on the 10X Genomics Chromium instrument (Genome Library Kit & Gel Bead Kit v2 PN-120258, Genome Chip Kit v2 PN-120257, i7 Multiplex Kit PN-120262). We sequenced this 10× library on an Illumina HiSeq X instrument with 150 bp read length to ~60× coverage. For optical mapping, the extracted
DNA (~750 µg) was labeled with a direct labeling enzyme (DLE-1) following the BioNano Prep Direct Label and Stain (DLS) Protocol (Document Number 30206). Labeled samples were imaged on the BioNano Saphyr instrument. Finally, Hi-C crosslinks were generated by Arima Genomics (https://arimagenomics.com/) using the Arima-Hi-C kit (P/N: A510008). From size selected fragments, Illumina-compatible libraries were generated using the KAPA Hyper Prep kit (P/N: KK8504). The resulting libraries were sequenced on an Illumina HiSeq X instrument to ~70× coverage.

With 10× long-linked reads, BioNano, and Hi-C maps in hand, the earlier polished and purged haplotigs were scaffolded in 3 stages according to Rhie et al. [2021]: first, we used the 10× linked-reads in 2 rounds of scaff10x (v4.1.0) (https://github.com/wtsi-hpag/Scaff10X) to generate the primary scaffolds. Second, we generated BioNano cmaps and used BioNano Solve (v3.2.1_04122018) [Lam et al., 2012] for hybrid scaffolding and to break misassemblies. Third, we used Salsa2 (v2.2) [Ghurye et al., 2019] to generate chromosomal-level scaffolds using the molecular contact information from Hi-C linked reads. Finally, we performed a second round of Arrow polishing on the maternal and paternal scaffolds with the binned long reads. During this round of polishing, gaps between contigs were closed by the gap-filling function of Arrow. The two haplotypes were then combined in a single assembly and underwent 2 rounds of short-read polishing using longranger (v2.2.2) [Bishara et al., 2015] and freebayes (v1.3.3) [Garrison and Marh, 2017]. After separating the scaffolds back into their respective haplotypes and removing the MT genome from each assembly, the two phased assemblies underwent manual curation using gEVAL as described previously [Chow et al., 2016; Howe et al., 2021], particularly to correct structural assembly errors.

## Assembly Statistics and Evaluation

Following each stage of the assembly, we calculated various metrics of assembly quality, for example, N50 contig length, number of contigs, and quality value (QV) scores for each base call to assess progress. We used Mercury (v1.0) for overall assembly evaluations (including k-mer completeness and spectra copy number analysis) as well as phasing assessment with hap-mers. We first generated 21-mer databases (dbs) from the raw F1 10× data and the parental Illumina data using meryl. We then built inherited hap-mer dbs by taking the difference between the maternal and paternal k-mer dbs, filtering according to the filter level used by TrioCanu for binning,

intersecting both with the F1 dbs, and filtering again, as below (steps 1–4). For evaluation of genome completeness and protein-coding gene representation, we ran BUSCO v4.0.2 [Manni et al., 2021] on our phased assemblies to determine the representation of near-universal single-copy orthologs in the vertebrate avian lineage (n = 8,338); aves_odb10 (online suppl. Material 1, Table 3).

## Genome Synteny and Structural Variation

To estimate sequence structural changes between assemblies for synteny, structural variation, and repeat expansion and contractions we used SyRi [Goel et al., 2019] with default parameters or Assemblytics v1.2.1 [Nattes-tad and Schatz, 2016] with a unique sequence length requirement of 10,000 on nucmer alignments between GRCg6a, GRCg7b, and GRCg7w assemblies.

## Gene Annotation

Both assemblies, GRCg7b, and GRCg7w were gene annotated using the standard NCBI pipeline [Pruitt et al., 2014], including masking of repeats prior to ab initio gene predictions, for evidence-supported gene-model building. All annotation processes used publicly available RNA-seq and Iso-Seq data from diverse tissue sources. We relied on the NCBI gene annotation report release 106 to compare the outcomes for each assembly. GRCg6a gene annotation data were reported earlier in NCBI release 104 using the same process as above.

## Interspersed Repeat Estimation

Two independent assessments were made to estimate the percentage of repeats to confirm their similarities between assemblies. RepeatMasker v4.0.9 [Smit, 2013] with -excln and -species chicken was used to identify and annotate repetitive regions of each genome while ignoring gap sequence; then WindowMasker analysis was carried out using default parameters [Morgulis et al., 2006].

## WGS and RNAseq Mapping

WGS short-read data of 6 chicken samples; a male and female layer, broiler and Ethiopian indigenous chicken breed were used to compare the mapping rate across the 3 genome assemblies (online suppl. Material 1, Table 6). WGS were first checked for quality using Fastqc [Andrews, 2010]. Trimmomatic [Bolger et al., 2014] was used to remove the remaining Illumina adapter sequences and low-quality bases with default parameters. Clean reads were mapped to the 3 reference assemblies (GRCg6a, GRCg7b, and GRCg7w) using bwa-mem with default parameters [Li and Durbin, 2009]. Picard (http://broadin-
stitute.github.io/picard/) was used to sort the mapped
files and merge files from multiple sequencing runs and
to mark duplicate reads. Finally, SAMtools [Li et al., 2009]
was used to assess mapping quality.

RNAseq data from 8 chicken samples were used to
compare the mapping rate across the 3 genome assem-
blies. We retrieved all sequence data from the NCBI Se-
quence Read Archive that included the diverse tissue
sources of ileum, bone-derived macrophages, uterus and
muscle from a male and female layer (online suppl. Mate-
rial 1, Table 9). Sequencing quality was checked by FastQC
software (v 0.11.7), qualifying reads were mapped using
STAR software (v 2.5.3a) with default parameters to all
assemblies (GRCg6a, GRCg7b, and GRCg7w), and the
percentage of uniquely mapped reads, multiple mapped
reads, and reads mapped to too many loci were taken for
the comparison. Moreover, the resulting bam files were
used for assessing the mapping rate for each sample with
the Samtools (v 1.9) "flagstat" command [Li et al., 2009].
Percentages of correctly paired reads were used for com-
parison.

SNV Analysis

To estimate SNV differences in the starting reference
alignments we used short-read sequences from a small
cohort of broilers (n = 10) representing commercial birds
generated by Cobb-Vantress (available upon request). All
samples attained genome coverage depth greater than
20× and individual reads were aligned to each reference
with the Nvidia Clara Parabricks (version 3.6) implemen-
tation of the BWA algorithm. Variants were called in
GVCF mode with Nvidia Parabricks HaplotypeCaller,
and GVCF files were loaded into GenomicsDB using
GATK 4.2.0 [Poplin et al., 2017], GenomicsDBImport,
and joint-genotyped with GATK's GenotypeGVCFs.
Hard-filtering was performed on the resulting raw VCF
using GATK's current best-practices for filtering.

BCFTools 1.12 was used to extract statistics on SNVs
and insertions and deletions (indels) per chromosome.
Variants that did not pass the filtering criteria were re-
moved, and mapping data for chromosomes were com-
pared between all assemblies using command-line tools
and then plotted in R. Ideograms were generated using
karyoplotER (v1.16.0) in R. Colored regions of the chro-
mosome denote annotated feature regions for that chro-
mosome. Rainfall plots of variants depict where variants
were found in the analysis along each chromosome. Each
unique color indicates a different type of substitution. We
only include variants that passed all filters and were het-
erozygous in the reference source.

ALVE Annotation

Assembled Avian Leukosis Virus subgroup E (ALVE) inte-
grations were identified by BLAST v2.10.0 [Altschul et al., 1990]
using the ALVE1 reference sequence (GenBank: AY013303.1)
and annotated for ORFs and miR-155 recognition sites [Hu et
al., 2016]. Analogous GRCg6a locations were identified using
flanking sequences, then compared with known ALVE integra-
tion sites and target site duplications (TSD) [Mason et al., 2020b,
c]. ALVE susceptibility was assessed by identifying the TVB re-
ceptor (TNFRSF10B) genotype.

Acknowledgments

All figures were generated with the BioRender software.

Statement of Ethics

For these experiments we used CO₂ gas euthanasia, following
the current standards for poultry euthanasia provided in AVMA
Guidelines for the Euthanasia of Animals (2020 Edition). All ex-
periments presented herein were carried out in accordance with
the approval of the Institutional Animal Care and Use Committee,
University of Arkansas, Fayetteville, AR. Moreover, all methods
were performed in accordance with the ARRIVE guidelines.

Conflict of Interest Statement

All authors have no conflicts of interest.

Funding Sources

This work was supported by USDA NIFA 2020-67015-31574 to Y.D. at the Western University of Health Sciences, and USDA NIFA 2022-67015-36218 to W.C.W., at the University of Missouri. The work of T.D.M. and V.A.S. was supported by the National Center for Biotechnology Information of the National Library of Medicine (NLM), National Institutes of Health.

Data Availability Statement

We have deposited the primary data underlying these analyses as follows: genome assemblies are deposited in the NCBI assembly archive (GRCg7b - PRJNA660757 and GRCg7w - PRJNA660758), PacBio SMRT reads associated with each reference are found in the SRA under the Bioproject number PRJNA673216. In addition, all sequence types and files are available in the VGP database (https://genomeark.s3.amazonaws.com/index.html?prefix=species/Gallus_gallus/bGalGal1/). The mitochondrial genome is available under NCBI accession number NC_053523.1. All Illumina data used in evaluating mapping rates, WGS or RNAseq, are described in various supplemental tables.

Karyotypes and Chromosome-Level Assemblies in
Chickens and Other Dinosaurs

(Prepared by D.K. Griffin, D.M. Larkin, and R.E.
O'Connor)

Any article about chicken, or avian generally, genomics
is inherently about dinosaur genomics. In the light of re-
cent paleontological evidence, statements such as “birds
evolved from dinosaurs” and “birds are *related to* dinosaurs” require substantial revision. Rather, it is very clear that birds, including our humble chicken, are in fact ex- tant dinosaurs. Omnipresent in literature, film, television, popular culture and the media since the original fossil dis- coveries, notions that dinosaurs were obliterated entirely by the latest mass extinction event have undergone radical revision. On the contrary, dinosaurs are the great survi- vors of mass extinction events and we suggest, this may be due, at least in part, to their unique genome organization (in other words, their karyotype). Studies of chromosome- level genome assemblies guided us to this conclusion.

*Chromosome-Level Assemblies*

Genomic approaches such as array comparative ge-
nomic hybridization (CGH) and next-generation se-
quencing (NGS) for cytogenetic application would not be
possible for many clinical and veterinary uses were it not
for “chromosome-level” assemblies (CLAs). In other
words, one ultimate objective of any genome assembly is
that all the sequences are correctly aligned and assigned
to their place on the appropriate chromosome, chromo-
some arm, and chromosome band.

Similarly, genome assembly needs cytogenetic analy-
sis. As Lewin et al. [2010] put it “every genome needs a
good map.” Thankfully, the outline of that map for any
given species is provided in the shape of a karyotype.
While karyotypes typically are used to detect chromo-
somal disease in humans, the karyotype can be consid-
ered the most basic low-resolution genomic map of an
organism. When whole-genome assemblies are yet to
reach the heights of a CLA (probably the case for most
sequenced genomes), their applicability for evolutionary
genomics is impeded. To give one of many examples,
CLAs are crucial for furthering applied agricultural re-
search in that an established order of sequences is an es-
sential pre-requisite to establish genotype-phenotype as-
sociations, e.g., in genome-wide association studies
(GWAS). CLAs have facilitated these in chicken, turkey,
and duck as well as in mammals such as pig, cattle, and
sheep. Other advantages of CLAs include discovery of
mendelian traits, spotting cryptic chromosome translo-
cations, as well as isolating quantitative trait nucleotides
(QTNs), expression quantitative trait loci (eQTLs), and
long-range regulatory interactions. The ultimate goal of
these studies is increased food production efficiency as
well as global food security. With the knowledge at our
fingertips for a number of species' CLAs, comparative ge-
nomics is much more practicable in silico, and intractable
karyotyping of chromosome rearrangements can be de-

tected by fluorescence in situ hybridization (FISH). Com-
parative genomics, thereafter, allows us to describe the
genome structure of more obscure species (by compari-
son to a standard, such as chicken) and the identification
of chromosome rearrangements that led to each species'
unique karyotype. CLAs also facilitate addressing funda-
mental biological hypotheses pertaining to genome evo-
lution, e.g., the mechanisms of chromosome breakage
and fusion, as well as the significance and genomic cor-
relates of evolutionary breakpoint regions (EBRs) and
homologous synteny blocks (HSBs). In the era of genom-
ics, cytogenetics (or, more specifically, cytogenomics) is
not just a descriptive science but it provides a framework
for the conceptualization of the structure of any genome.
It gives us a starting point through which we can under-
stand genome-phenome correlations.

The particular challenge with establishing genome
structure in birds is that while, in most species, correlat-
ing genome assembly with karyotype is a bit like identify-
ing cities and towns in the major landmasses of Europe
(some big, some small), doing it in birds is more like do-
ing it for Polynesia where there are a lot of tiny islands,
i.e., the microchromosomes. In our research, we have
looked at CLAs in a multitude of bird species, providing
novel insight into the genome organization of the avian
forebears – the extinct dinosaurs.

What Is a Dinosaur?

The official definition of a dinosaur is “Triceratops,
modern birds, their most recent common ancestor and all
their descendants.” As biologists however, perhaps an
easier way to visualise them is as reptiles (in this defini-
tion we class birds as reptiles) with hind limbs held erect
beneath the body (like mammals). This distinguishes
them from most other reptiles such as lizards and croco-
dilians where the legs are positioned to the side. If we then
eliminate the sister group, the pterosaurs, which are eas-
ily distinguishable, we can usually spot a dinosaur with
relative ease from another organism.

Accepting then that modern paleontological evidence
is explicit that birds are dinosaurs, rather than consider-
ing them as a group of animals that the Chicxulub meteor
wiped out, dinosaurs are, in point of fact, survivors of nu-
merous extinction events including the latest – the K-Pg
(Cretaceous-Paleogene) extinction event. This segment
describes how groups mostly in Kent and the Royal Vet-
erinary College combined molecular cytogenetics and
bioinformatics to establish that this resilience and ability
to recover from extinction event may be due, at least part-
ly, to their karyotype.
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Evolution of major groups of reptiles (including birds) with major extinction events noted. Timelines given but scale is not linear.</figcaption></figure>

*The Dawn and Dusk and Dawn of the Dinosaurs*

As illustrated in Figure 6, around 325 million years ago, the amniote lineage split into synapsids (ultimately becoming mammals, amongst other groups) and the rep- tile/bird lineage – diapsids. There are over 17,500 living diapsid species, the majority of which (around 11,000) are birds. Dinosaurs (including birds), pterosaurs, turtles, and crocodilians all share a common ancestor that lived 275 million years ago [Shedlock and Edwards, 2009; Hedges et al., 2015], with the turtles (Testudines) diverg- ing first (around 255 million years ago), the crocodilians around 252 million years ago, the pterosaurs about 245 million years ago, and the true dinosaurs about 240 mil- lion years ago. For the next 30 million years, dinosaur species were few in number, but amid the Jurassic period, this number, their geographical spread, and their body size all dramatically increased [Benton et al., 2014]. The proceeding 135 million years of dinosaur evolution is spectacular for being a time when they had an incredible range of species diversity and thus became the dominant vertebrates on the planet. Once human beings evolved scientific investigation, popular culture and media how- ever, then the dinosaurs' legacy was complete. Incredibly, dinosaurs survived the Carnian-Norian and End-Triassic mass extinction events (228 and 201 million years, re- spectively) and now there are >1,000 discovered fossil species. About 30 more appear annually, not including birds, in the fossil record [Weishampel, 2004]. The dev- astation of the K-Pg extinction event 66 million years ago nearly wiped them out, but they bounced back again as

modern birds, with more species than any other terres- trial vertebrate (Fig. 6). The cytogenomic study of birds is an independent line of enquiry to paleontology and cir- cumvents some of the problems associated with fossil dat- ing.

The remarkable diversity and species abundance that is observed in dinosaurs is often put down to the fact that competitor species were wiped out, thereby allowing the dinosaurs to thrive. It has nonetheless also been proposed that such impressive levels of abundance and diversity reflect genomic adaptations that are particular to dino- saurs, facilitating their survival over other species in harsh environments. Some examples include extraordinary bone growth rates and highly evolved respiration systems [Farmer and Sanders, 2010], including unidirectional respiration [O'Connor and Claessens, 2005]. These sorts of evolutionary adaptations could have led to the evolu- tionary success of avian species, clues of which may be found in their genome structure and organization.

*Dinosaur (Bird) Evolution Just before, during, and after the Last Mass Extinction Event as Revealed by Multiple Genome Sequencing Efforts*

As a result of a multiple genome sequencing effort in birds [Jarvis et al., 2014; Zhang G et al., 2014], a revised avian phylogeny based on genome assemblies (some of which were CLAs) corrected the timing of avian diversi- fication. We can now consider the first avian evolutionary divergence at about 100 million years ago with the Pale- ognathae (Ratites/Tinamous) branching from the Neog-
nathae (Galloanseres/Neoaves). The second is where Gall-
loanseres (Galliformes and Anseriformes) and Neoaves
then diverged 80 million years ago, with Galliformes
(landfowl, e.g., chicken, turkey, quail, pheasant) and An-
seriformes (waterfowl, e.g., geese, ducks, swans) diverg-
ing about 66 million years ago. Another major divergence
of the Neoaves into Columbea (e.g., pigeons) and Pas-
sarea (e.g., songbirds) is dated slightly earlier (67–69 mil-
lion years ago). Data from the Jarvis et al. [2014] analysis
plus that of Prum et al. [2015] suggest that, following the
K-Pg mass extinction event [Schulte et al., 2010], around
the same time of these two divergences, there was a rapid
period of diversification, with 36 lineages evolving in a
very short evolutionary period of 10–15 million years
[Jarvis et al., 2014]. Genomics has thus updated our un-
derstanding of dinosaurs through comparative studies,
proving intriguing insights into relationships with diver-
sity and phenotype [Jarvis et al., 2014; Zhang G et al.,
2014]. The karyotype of dinosaurs was something there-
fore that warranted deeper investigation.

Dinosaur Karyotype Evolution

The extraction of intact DNA from Jurassic blood-
sucking insects is an outstanding plot device for novelists
and film-makers but, alas, not a feasible means of facilitat-
ing the making of metaphase preparations. We can, how-
ever, with enough avian CLAs, glean insight into extinct
dinosaur karyotypes by inference. Analysis of (near)
CLAs from 6 living birds plus an *Anolis* lizard outgroup
allowed us to infer the most likely ancestral karyotype of
all birds [Romanov et al., 2014]. Generally speaking, the
common avian ancestor was probably a bipedal, terres-
trial, chicken-sized small Jurassic dinosaur with some fly-
ing ability [Witmer, 2002], and our studies established
that its karyotype was very similar to that of a chicken or
a ratite bird. We then went on to reconstruct the most
parsimonius sequence of events that led to the karyotypes
we typically observe in avian species. The humble chicken
(*Gallus gallus*) is in fact the closest, karyotypically, to the
reconstructed ancestral pattern among the birds that we
studied, whereas zebra finch (*Taeniopygia guttata*) and
budgerigar (*Melopsittacus undulatus*) appear to have un-
dergone the most intra- and inter-chromosomal changes,
respectively [Romanov et al., 2014]. We later reconstruct-
ed the ancestral avian karyotype using an algorithmic ap-
proach, DESCHRAMBLER, on fragmented genome as-
semblies. There [Damas et al., 2018], we performed large-
scale analysis of ancestral avian chromosome structure
around 14 key nodes of bird evolution. Our results pro-
vided insight into the variability in the rates of rearrange-

ment that occurred during avian evolution. It also al-
lowed us to detect patterns related to chromosome distri-
bution of EBRs and microchromosomes.

In the same year, we applied a comparable approach to
recreate the most likely ancestral karyotype of diapsids
[O'Connor et al., 2019]. We developed a universally hy-
bridizing bacterial artificial chromosome (BAC) FISH
probe set that was able to hybridize directly across species
that diverged hundreds of millions of years ago [Damas
et al., 2017]. The BAC probes employed in FISH experi-
ments gave strong signals on lizard (*Anolis carolinensis*)
chromosomes and more so on chromosomes of the red-
eared slider (*Trachemys scripta*) and spiny soft-shelled
turtle (*Apalone spinifera*). FISH experiments then al-
lowed us to anchor the series of events from the perspec-
tive of an archelosaur (bird-turtle) ancestor. By combin-
ing molecular cytogenetics and bioinformatics, we then
recreated the cytogenetic changes that occurred from the
ancestral diapsid ancestor, through to the archelosaur an-
cestor [Benton et al., 2015], through the theropod lineage,
through to birds – including of course chicken. Screen-
shots of these events are depicted in Figure 7 in the panels
that represent an animation (the video of which is avail-
able in online suppl. Material 3).

Hybridization of BACs to *T. scripta* (2n = 50) and *A. carolinensis* metaphases (2n = 36) also unveiled larger chromosomes with microchromosomal homologues attached, giving clues to the ancestral pattern of the diapsid ancestor (see Fig. 6). Figure 7 thus depicts a diapsid ancestral karyotype (275 million years ago) with 2n = 36–46 (50:50 macro- and microchromosomes) [Beçak et al., 1964; Alföldi et al., 2011]. This underwent rapid change over ~20 million years and we established that most of the major features associated with a typical bird karyotype were already laid down in the archelosaur ancestor 255 million years ago. This is because most chicken (for this, read ancestral avian) chromosomes (numbered 1–28 + Z) are precisely syntenic to those of *A. spinifera* (2n = 66). Analogous studies using chicken chromosome painting on Chinese soft-shelled turtle (*Pelodiscus sinensis*) (2n = 66) [Matsuda et al., 2005], *T. scripta* [Kasai et al., 2012], and painted turtle (*Chrysemys picta*) chromosomes (both 2n = 50) [Badenhorst et al., 2015] further point to the notion that bird and turtle macrochromosomes are precise counterparts of one another. From this basic pattern that was laid down 255 million years ago in the archelosaur ancestor, only about 7 fissions would be needed to establish the familiar karyotype pattern directly observed in most of the major groups of birds including the Ratites, Galliformes, Anseriformes, Columbaea, Pas-
seriformes, and others. Determining the precise timings that these “final” chromosome fissions occurred cannot be easily achieved with the available evidence. However, if the same fission rate that had occurred for the previous 20 million years continued for another 15 million years, a complete bird-like karyotype would have emerged before the appearance of the earliest dinosaurs and pterosaurs 240 million years ago [Baron et al., 2017].

To paraphrase, our evidence strongly suggests that *the chromosomal pattern that we see in the majority of bird species that we choose to karyotype has remained mostly unchanged, not only in most birds, but with a reasonable degree of certainty, in many, if not most, extinct dinosaurs too* [O'Connor et al., 2019 and unpublished results]. We thus go so far as suggesting that if we had the opportunity to make chromosome preparations “a la Jurassic Park/World” from the blood of extinct dinosaurs, then karyotype analysis and zoo-FISH results would differ very little from that of a modern chicken (Fig. 8). Those most closely related to modern birds such as the theropods (Spielberg favorites *Tyrannosaurus rex* and *Velociraptor* are representatives) are the most likely to have a very avian-like karyotype.

*Karyotype, the Reduction in Genome Size, and the Evolution of Flight*

Dave Burt [Burt, 2002] suggested that some avian microchromosomes were present in the avian ancestor >80 million years ago [Cracraft et al., 2015], purporting that it maybe had a karyotype of around $2n = 60$. A very worthy summation of the evidence at the time of writing. As detailed above, however, on the basis of our new evidence over 15 years later, we challenged this notion. We also challenged the notion that this fragmented genome organization (i.e., a karyotype with $2n = 80$ chromosomes) accompanied the genome size reduction in birds that has sometimes been associated with the evolution of flight [O'Connor et al., 2018b]. That is, it had previously been suggested that there was some evidence of a correlation between genomes with fewer chromosomes (and no microchromosomes) and larger genome sizes (2.5–3 Gb).

**Fig. 7.** Chromosome evolution from the diapsid ancestor, to the archelosaur ancestor, via the theropod dinosaur lineage, to modern birds and finally to chicken. A recognisable avian pattern had evolved just before the dinosaurs emerged 240 million years ago. After this time, chromosome inversions were the principal mechanisms of change. **a** 275 million years ago. **b** 255 million years ago. **c** 240 million years ago at the dawn of the dinosaurs. **d** Snapshot 100–150 years ago. **e** Modern chicken.

<figure><img src="image_9.png" /></figure>
such as is seen in mammals and crocodilians [St John et al., 2012; Kapusta et al., 2017]. We however suggest that the avian karyotype was in place first, and that this was followed by a reduction in genome size, which was then followed by the evolution of flight.

## Why Not Change?

Why then is a near identical karyotype pattern present in most birds, and why has it been in place for around 255 million years? Conceptually there could be two reasons: that is, either there is little opportunity for change and/or the configuration is so evolutionarily successful that there is no pressure to change. In the case of the former, repetitive elements provide substrates for interchromosomal rearrangement, often seen in mammals but hardly even seen in avian species, suggesting that the bird karyotypes provide fewer opportunities for interchromosomal rearrangement because there are fewer recombination hotspots [Kawakami et al., 2014; Smeds et al., 2016], repeat structures [Mason et al., 2016; Gao et al., 2017; Warren et al., 2017], or endogenous retroviruses [Cui et al., 2014; Romanov et al., 2014; Farré et al., 2016]. We also provided evidence of purifying selection acting on some of the smallest microchromosomes [Damas et al., 2018]. A karyotype that hardly changes for 255 million years also however suggests that it is an evolutionarily successful one. The large chromosome number (especially microchromosomes) with high rates of recombination, could, we hypothesise, be the cause of the great variation we see in dinosaurs (including birds) mediated through random chromosome segregation and increased genetic recombination. Phenotypic variation is, of course, the driver of evolution and, although having many chromosomes is not the only mechanism through which variation can be generated, it may nonetheless explain this apparent paradox in dinosaurs of remarkable phenotypic diversity but very little karyotypic diversity. We should nonetheless recognise that it is possible (indeed likely) that some dinosaurs underwent a lot of interchromosomal change. Kingfishers [Christidis, 1990] (many fissions), parrots

[Nanda et al., 2007; O'Connor et al., 2018a], and falcons [Damas et al., 2017; Joseph et al., 2018] (many fusions) are modern examples of where this has occurred. Which specific extinct dinosaur groups did this, may however always remain a mystery.

## Chromosome Inversion, the Role of Gene Ontology Analysis, and Dinosaur Phenotype

In the absence of interchromosomal change, the principal mechanism for chromosomal change in dinosaur genome evolution was probably chromosome inversion (also depicted in Fig. 7). Using the ancestral genome re-construction tool Multiple Genome Rearrangement and Analysis (MGRA) [Avdeyev et al., 2016], we generated contiguous ancestral regions likely to represent the chromosomes of the diapsid ancestor. Compared to extant birds, we identified inversions along the path from the diapsid ancestor to the modern chicken, probably underestimating this number of inversions at 49 chromosome inversions. We believe that the rate of intrachromosomal change increased in modern times, even in the chicken [Romanov et al., 2014]. An even greater degree of change was seen however in some bird clades, particularly the songbirds [Skinner and Griffin, 2012; Zhang G et al., 2014; Farré et al., 2016], the group with the most species. It seems reasonable to hypothesise therefore that periods of faster speciation may have also been accompanied by increased chromosome inversion rates in other dinosaur groups [Skinner and Griffin, 2012; Romanov et al., 2014; O'Connor et al., 2018b].

We [O'Connor et al., 2018b], identified around 400 HSBs delineated by EBRs that characterize dinosaur genome evolution. Other genomic studies in other species (mostly mammals) established that EBRs commonly occur in gene-dense loci, with genes related to lineage-specific biology, transposable elements and other repetitive sequences [Pevzner and Tesler, 2003; International Chicken Genome Sequencing Consortium, 2004; Larkin et al., 2009; Rao et al., 2012]. HSBs on the other hand contain more developmental genes and regulatory elements

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Imagined dinosaur karyotype (deliberately manipulated image based on chicken and spiny soft-shelled turtle chromosomes).</figcaption></figure>
[Larkin et al., 2009; Warren et al., 2017]. Regions more likely to break, e.g., open chromatin areas or recombination hotspots, and chromosome breaks that do not disrupt key genes or provide a selective advantage, are more likely to be fixed in populations [Farré et al., 2016].

HSB analysis [O'Connor et al., 2018b] using gene ontology (GO) tools, established significant enrichments relevant to amino acid transmembrane transport and signalling plus synapse/neurotransmitter transport, nucleoside metabolism, cell morphogenesis, and cytoskeleton and sensory organ development. Former studies established that HSBs are enriched for GO terms related to evolutionary constant phenotypic features [Larkin et al., 2009], and our dinosaur results support this hypothesis. EBRs on the other hand, are often proposed to be where the "action" in genome evolution occurs [Sankoff, 2009]. We initially found GO terms in avian EBRs that were associated with specific adaptive features, e.g., enrichment for forebrain development in the budgerigar EBRs (consistent with vocal-learning) [Farré et al., 2016]. Later, we identified significant enrichments in genes and single GO terms pertaining to chromatin modification, chromosome organization, and proteasome/signalosome structure [O'Connor et al., 2018b].

The discovery that the avian karyotype likely dates back before the dawn of the dinosaurs complements paleontological research that demonstrates that features such as feathers and pneumatised skeletons arose first among more ancient dinosaur or archosaurian ancestors [Zhou, 2004; Baron et al., 2017]. Dinosaurs were the dominant group of animals for around 200 million years, with significant radiations occurring in response to two mass extinction events and, despite being almost wiped out by a third (K-Pg), their resilience as a highly diverse and speciose clade (extant birds) [Barrowclough et al., 2016] is evident.

## Conclusions

Far from being simply a curating exercise (and a speculative one at that), the study of likely dinosaur karyotype sheds new light into genome evolution, with distinct clues about phenotype and an alternative line of enquiry compared to more established methods. A highly fragmented genome that appears like it has been hit by a meteor is ironic since its pattern was most likely established around 200 million years before Chicxulub hit the earth. Our studies reveal an apparent paradox of an organization that is remarkably unchanging karyotypically during evolution, yet, quite possibly, the driver of such phenotypic evolutionary change. Did the fictional creators of Jurassic

Park ever attempt karyotyping? Well, apparently, they did: Check out the last of the Jurassic World movies: at one point, an historical recording of the (now deceased) character Charlotte Lockwood appears, observed by her young daughter Maisie. In the clip, there is a screen with chromosomes on them in the background. There are however no microchromosomes! Unless Spielberg knows something that we don't, these chromosomes were unlikely to be dinosaur in origin.

## Conflict of Interest Statement

The authors have no conflicts of interest to declare.

## Funding Sources

This research was funded by the Biotechnology and Biological Sciences Research Council (BB/K008226/1 and BB/J010170/1 to D.M.L, and BB/K008161/1 to D.K.G).

## Chicken as a Model for Bird Genome Alignments

(Prepared by P.D. Waters, J.A.M. Graves, and H.R. Patel)

Chicken has long been a model for bird genomes [International Chicken Genome Sequencing Consortium, 2004], but it is only in recent years that a chromosome-level assembly has been available that includes the gene-dense microchromosomes.

The chicken karyotype is characterised by 9 macrochromosomes (defined as >35 Mb in the genome assembly), which include the Z, and 30 microchromosomes (defined as <35 Mb) that are characteristic of bird genomes [Mendonça et al., 2016]. Like all birds, chicken has a ZW female/ZZ male sex chromosome system, in which sex is determined by dosage of the gene *DMRT1* on the Z chromosome [Smith et al., 2009].

Both the macro- and microchromosomes are conserved across the most distantly related bird lineages [Waters et al., 2021], with the exception of some clades that have undergone rearrangement [Nanda et al., 2006; Huang et al., 2022] (Fig. 9). In fact, strong homologies are conserved with reptiles (lizards, snakes, turtles, and alligators), with differences mostly explained by chromosome fusions (less often fissions) so that large regions of synteny are retained [Waters et al., 2021].

Strikingly, it was demonstrated that almost all amphioxus chromosomes shared homology with one or two bird microchromosomes, in addition to two or three regions of macrochromosomes [Waters et al., 2021]. These multiple regions of amphioxus homology to the bird genome are likely to represent the four copies resulting
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Bird whole-genome alignments. **a** Alignment of chicken to zebra finch [Rhie et al., 2021], Anna's hummingbird [Rhie et al., 2021], superb fairywren [Penalba et al., 2020], and the ancestral emu genome [Liu J et al., 2021]. **b** Alignment of the chicken genome to species that have undergone significant chromosome rearrangement: saker falcon, California condor [Robinson et al.,</figcaption></figure>

2021], and golden eagle [Mead et al., 2021]. Macrochromosomes are green, microchromosomes purple, and Z chromosomes orange. Green ribbons show macro-to-macro homologies, purple show micro-to-micro homologies, blue show macro-to-micro homologies, and grey show Z-to-Z homologies.

from two rounds of whole-genome duplication in vertebrates [Simakov et al., 2020]. The fact that each amphioxus chromosome shared striking homology with one or two bird microchromosomes suggests that each bird microchromosome represents one copy of an ancestral chromosome.

Microchromosomes are CG-rich and gene-dense [Waters et al., 2021], as are amphioxus chromosomes. Additionally, microchromosomes cluster together in the nucleus and have many more interactions with each other than with the macrochromosomes [Liu J et al., 2021; Waters et al., 2021]. Therefore, they are physically isolated from macrochromosomes, so are less likely to fuse with them.

It is simplest to envisage that after genome duplication one complete copy was rapidly packaged up in the interior of the nucleus and isolated from the other genome duplicate, rather than a mix of chromosomes from duplications being isolated. Given their tight association in the nucleus, it is somewhat surprising that microchromosomes have not all fused with each other in birds, as they have in crocodiles and mammals. This may have resulted

from the selective advantage of higher recombination due to obligate cross overs on each chromosome and their independent segregation at meiosis.

Here we report whole-genome alignments of the chicken genome with those of other bird representatives with chromosome-level assemblies. The chicken genome aligns very well with the emu genome, which is considered to represent the ancestral bird genome [Liu J et al., 2021]. The only whole-chromosome difference is the fusion of a microchromosome to chromosome 4, uniquely in chicken [Shetty et al., 1999; Burt, 2002]. There have been some small internal rearrangements at the termini of macrochromosomes, two that are chick-specific and three that are shared with other carinate clades so presumably occurred in the ancestor of all carinates. The Z chromosome has two large internal rearrangements with respect to the emu Z, but these are shared by other carinates.

As a gold standard genome, we propose that the chicken assembly is the best model for the bird genome. It has few autosomal rearrangements relative to the ancestral genome, and four of these are shared by other carinates.
All microchromosomes are represented in the most recent assembly. Therefore, whole-genome alignments to it can be used to assess assembly errors and/or real rearrangements in new bird assemblies for which karyotype information is available. This presents a more robust comparison than another popular bird model, the zebra finch, which has undergone macrochromosome fission and internal rearrangement (Fig. 9). It is also more useful for aligning the genomes of birds with very rearranged chromosomes, since the multiple rearrangements are independent in different clades (e.g., golden eagle and falcon/parrot, Fig. 9).

*Conflict of Interest Statement*

There are no conflicts of interest to declare.

*Funding Sources*

P.D.W. is supported by Australian Research Council Discovery Projects (DP180100931, DP210103512, and DP220101429). J.M.G. is supported by Australian Research Council Discovery Projects (DP210103512 and DP220101429). H.R.P. is supported by the Australian National University and NHMRC Synergy Grant (APP2011277). The authors acknowledge the provision of computing resources provided by the Australian BioCommons Leadership Share (ABLeS) program. This program is co-funded by Bioplatforms Australia (enabled by NCRIS), the National Computational Infrastructure and Pawsey Supercomputing Centre.

**The Avian W Chromosome: Simply a Y with a Different Name?**

(Prepared by P.D. Price, T.F. Rogers, and A.E. Wright)

Sex chromosomes have long fascinated biologists due to their unique gene content and evolutionary trajectories relative to the rest of the genome [Furman et al., 2020]. In particular, the halting of recombination between sex chromosome pairs has resulted in the evolution of highly degenerate sex-limited W and Y chromosomes in many species [Charlesworth, 1991]. Identifying the function of these chromosomes and understanding if and how they can resist the degenerative forces arising from reduced recombination has been the focus of numerous studies [Bachtrog et al., 2011].

We now know a considerable amount about Y chromosomes, despite the difficulties in sequencing highly heterochromatic and repetitive genomic regions [Tomaszkiewicz et al., 2017]. Their evolution is typically characterised by the accumulation of genes with male-specific functions, large-scale gene amplification, and rapid turnover of gene content across lineages [Bachtrog, 2013; Subrini and Turner, 2021]. In contrast, our under-

standing of the W chromosome has lagged. However, the last decade has seen an explosion in the number of W-linked genes sequenced across birds [Zhou Q et al., 2014; Bravo et al., 2021], ranging from songbirds [Smeds et al., 2015; Xu et al., 2019; Sigeman et al., 2021; Huang et al., 2022; Warmuth et al., 2022] to fowl [Moghadam et al., 2012; Ayers et al., 2013; Wright et al., 2014] to paleognaths [Xu and Zhou, 2020; Liu J et al., 2021], and a reference assembly of the chicken W chromosome [Warren et al., 2017]. In theory, W chromosomes are in many ways comparable to Y chromosomes, as both are sex-limited and often don't recombine, and so they might be expected to share similar evolutionary fates. However, there are key differences, most notably that the W chromosome is limited to females whereas the Y chromosome is only present in males [Bachtrog et al., 2011; Mank, 2012]. Below, we outline new insights into avian W chromosomal evolution and ask whether W and Y chromosomes are really that different.

*What Are the Evolutionary Dynamics of W Chromosomes across Birds?*

It has been known for decades that the chicken W chromosome is a degenerated version of the Z, with the most recent build of the W reference (GRCg7b) identifying only ~80 protein-coding genes across ~9 Mb [Warren et al., 2017]. However, establishing whether the chicken W chromosome is representative of the avian W more generally has only recently been possible due to the plethora of W-linked sequences now available across the avian phylogeny [Wang Z et al., 2014; Zhou Q et al., 2014; Xu et al., 2019; Sigeman et al., 2020].

Sex chromosomes diverge as recombination is suppressed between them, typically assumed to occur in a stepwise process through sequential inversions [Charlesworth et al., 2005]. Consistent with this, 'strata' of different ages can be detected on avian Z and W chromosomes [Wang Z et al., 2014; Wright et al., 2014]. These strata are thought to reflect the halting of recombination through large-scale chromosomal rearrangements, such as inversions [Wright et al., 2016]. However, recent evidence suggests that recombination suppression at the earliest stages of avian sex chromosome divergence is a more mosaic and gradual process [Sigeman et al., 2021]. This cessation of recombination has occurred independently in different avian lineages [Zhou Q et al., 2014] and many species exhibit a heavily degraded W chromosome, similar to the chicken.

Despite degeneration proceeding independently across birds, the set of ancestral genes retained on the W chromosome is remarkably conserved, suggesting that
decay is non-random [Xu and Zhou, 2020]. For instance, over 80% of W-linked genes in the oldest stratum are conserved across chicken, songbirds and tinamous [Xu and Zhou, 2020]. This is in stark contrast to the Y chromosome, where frequent gene movement onto and off the chromosome is common [Hughes et al., 2015; Mahajan and Bachtrog, 2017]. This has been interpreted as a product of differing selective pressures acting on Y versus W chromosomes, with the W chromosome subject to stronger purifying selection compared to the Y due to its higher effective population size relative to the autosomes [Wright and Mank, 2013].

However, there is still remarkable variation in the extent of Z-W divergence across birds [Zhou Q et al., 2014]. In contrast to the chicken, the paleognath W chromosome recombines along a large proportion of its length and so has experienced limited decay, although this varies across species [Zhou Q et al., 2014; Yazdi and Ellegren, 2018; Liu J et al., 2021] with greater recombination suppression in ostrich and emu than tinamou. The growing amount of long-read sequencing data for birds has also revealed that fusions between sex chromosomes and autosomes to create neo-sex chromosomes are not uncommon, with two independent origins across Psittaciformes [Huang et al., 2022], four across Sylvioidea [Pala et al., 2012; Sigeman et al., 2020, 2021], one in the eastern yellow robin (*Eopsaltria australis*) [Gan et al., 2019], cuckoo (*Crotophaga ani*) [Kretschmer et al., 2020] and Raso lark (*Alauda razae*) [Dierickx et al., 2020] identified to date. No doubt this number will increase as more bird genomes are probed for sex chromosomes, making it possible to test the evolutionary pressures responsible for driving these fusions. Together, this challenges the traditional view that the avian W chromosome is genetically inert and highly conserved across species.

*Is the Avian W Chromosome Selected for Female-Specific Functions?*

Given the sex-limited inheritance pattern of Y and W chromosomes, theory predicts that they should be subject to sex-specific selection and accumulate genes with sex-specific functions [Rice, 1984]. Indeed, the Y chromosome in many species is enriched with genes predominantly expressed in testes that function in spermatogenesis [Bachtrog, 2013; Subrini and Turner, 2021], although there is a growing awareness of its role in non-reproductive traits [Cirulis et al., 2022]. It follows that we might expect the W to be subject to female-specific selection to retain genes with female fitness benefits.

The avian W chromosome lacks a candidate sex-determining gene [Schmid et al., 2015]. Instead, sex in birds is determined by dosage of the Z-linked gene, *DMRT1* [Hirst et al., 2017]. However, there are several lines of evidence implicating the avian W chromosome in female fertility, although the precise functions of genes on the W have yet to be defined. First, W-linked genes are highly expressed in developing chicken ovaries [Moghadam et al., 2012; Ayers et al., 2013; Xu and Zhou, 2020]. This is consistent with feminization of the W [Mank et al., 2010], as key stages of oogenesis are restricted to embryogenesis unlike spermatogenesis, which is a continuous process throughout adult life. Second, W-linked genes expressed during late female development are convergently upregulated in the ovaries of chicken layer breeds subject to artificial selection for fecundity relative to their modern ancestor, the Red Jungle Fowl, and other chicken breeds [Moghadam et al., 2012].

However, unlike most Y-linked genes, which typically exhibit testes-specific expression [Subrini and Turner, 2021], expression of genes on the avian W chromosome is not limited to the ovary. Instead, studies from chicken and collared flycatcher show that many W genes are active in both somatic and reproductive tissue [Smeds et al., 2015; Bellott et al., 2017; Xu and Zhou, 2020]. While this does not preclude a specific role of the W chromosome in oogenesis, it has led to suggestions that this chromosome has instead been selected to maintain gene dosage and ancestral expression levels of essential genes. Consistent with this, many avian W-linked genes are subject to purifying selection [Wright et al., 2014; Sigeman et al., 2021], exhibit a high degree of sequence conservation as well as similar expression patterns to their Z-linked partner [Ayers et al., 2013; Smeds et al., 2015; Xu and Zhou, 2020], and have human orthologs that exhibit detrimental effects when haploid [Bellott et al., 2017; Xu et al., 2019; Bellott and Page, 2021; Sigeman et al., 2021].

It is plausible that apparent differences in the function of Y and W chromosomes could arise from their contrasting inheritance patterns. For instance, W-linked genes, which only pass through the female germ line, are not exposed to sperm competition and so might be subject to weaker sex-specific selection than genes on the Y chromosome. However, it is worth noting that our understanding of the function of the avian W is based on expression data from a limited number of species (chicken and flycatcher) taken across whole, heterogeneous, adult tissue. This precludes accurate contrasts of expression between Z and W orthologs [Price et al., 2022] and so could lead to false inferences of selection to maintain gene dos-
age between gametologs. Further expression analyses, in-
corporating a broader taxonomic range and data for in-
dividual cell types throughout female development, are
essential to ascertain why specific genes have been re-
tained on the avian W chromosome.

*How Do Multi-Copy Gene Families Evolve on the W Chromosome?*

Y chromosome degeneration is frequently character-
ised by massive gene amplification where many remain-
ing Y-linked genes persist as members of multi-copy gene
families [Skaletsky et al., 2003; Soh et al., 2014; Bachtrog
et al., 2019; Vegesna et al., 2020]. However, until recently,
gene amplification on the W chromosome has received
comparatively less attention and it remained unclear
whether large-scale gene amplification is a general feature
of sex chromosome evolution or a peculiar quirk of the Y.

A handful of W-linked multi-copy gene families have
been identified in a limited number of avian species [Hori
et al., 2000; Backström et al., 2005; Davis et al., 2010;
Smeds et al., 2015; Rogers et al., 2021]. The most compre-
hensively studied of these is histidine triad nucleotide-
binding protein (HINTW), an ampliconic gene family
that is hypothesized to play a role in female reproduction
and oogenesis [O'Neill et al., 2000; Ceplitis and Ellegren,
2004]. At present, approximately 10 different copies of
HINTW are annotated in the most recent chicken W
chromosome assembly (GRCg7b), however, this is likely
an underestimation with previous studies estimating over
40 copies [Hori et al., 2000; Backström et al., 2005]. Fur-
thermore, large-scale amplification of HINTW is con-
served across avian non-ratites [Hori et al., 2000]. Cur-
rently, evidence for the functionality of HINTW is lack-
ing. However, it is known that HINT can form a
heterodimer and the amino acid residues that form the
dimer binding site are conserved in many HINTW cop-
ies, although many copies are nonfunctional [Hori et al.,
2000; O'Neill et al., 2000]. Therefore, HINTW might act
to disrupt the function of its Z-linked ortholog (HINTZ)
by forming a heterodimer. Interestingly, the size of the
HINTW gene family varies between chicken layer breeds
subject to artificial selection for fecundity relative to oth-
er chicken breeds [Rogers et al., 2021], potentially sug-
gesting a role of female-specific selection in driving gene
amplification, although this relationship was absent
across duck breeds.

The paucity of multi-copy gene families on the avian
W chromosome is in stark contrast to the abundance of
ampliconic genes often present on the Y chromosome.
Several mechanisms have been proposed to drive the evo-

lution of multi-copy gene families on the Y, including
meiotic drive, sperm competition, genetic drift, and gene
conversion [Skaletsky et al., 2003; Ellis et al., 2011; Coc-
quet et al., 2012; Larson et al., 2018; Bachtrog et al., 2019;
Vegesna et al., 2020]. In theory, the strength of these pro-
cesses might differ between the Y and W due to their con-
trasting inheritance patterns [Wright and Mank, 2013].
Notably, the Y chromosome is exposed to spermatogen-
esis, whereas the W is subject to oogenesis, and this likely
leads to differences in the potential for antagonistic co-
evolution between the X and Y versus the Z and W. An-
tagonistic co-evolution is predicted to drive the co-ampli-
fication of genes on sex chromosomes but should be
weaker during oogenesis than spermatogenesis, poten-
tially explaining the limited number of W-linked multi-
copy gene families [Bachtrog, 2020]. Targeted avian gene
knockouts [Ioannidis et al., 2021] provide an exciting op-
portunity to elucidate the functionality of HINTW cop-
ies, whether this varies across avian species, and the po-
tential for antagonism between W and Z orthologs.

*Is There a “Toxic W” Effect?*

There appears to be a cost for males to carrying a de-
generated Y chromosome [Brown et al., 2020; Xirocostas
et al., 2020; Nguyen and Bachtrog, 2021; Connallon et al.,
2022], where males in species with XY chromosomes tend
to die earlier [Xirocostas et al., 2020]. Several hypotheses
have been put forward to explain this phenomenon, in-
cluding the presence of deleterious recessive mutations
on the single X in males that would otherwise be shielded
in females (“unguarded X”) or the accumulation of muta-
tions and repetitive elements on the Y chromosome
(“toxic Y”). There is also evidence that the Y chromosome
acts as a heterochromatin sink, reducing the efficiency of
heterochromatin maintenance across the rest of the male
genome [Francisco and Lemos, 2014; Brown et al., 2020].

Similar processes may operate on the W chromosome,
where females exhibit a shorter lifespan than males across
a range of species [Xirocostas et al., 2020]. Consistent
with a “toxic W” effect, the avian W chromosome is a ha-
ven for repetitive material and transposable elements in
several species. For instance, females in species with a de-
generate W carry between 20 and 90% more endogenous
retroviruses than males [Peona et al., 2021]. Furthermore,
transposable element suppression is less effective on the
crow W chromosome than the rest of the genome, leading
to higher expression of transposable elements in females
[Warmuth et al., 2022]. Although transposable elements
can facilitate adaptive evolution, they also have the poten-
tial to reduce fitness through the disruption of gene activ-
ity and the promotion of chromosomal rearrangements [McDonald, 1993]. In theory, they may also contribute to an increased chance of female sterility in hybrids, where mechanistic mismatches between transposable repressor mechanisms and the W chromosome lead to reduced female fertility. This would provide further support for Haldane's rule where the heterogametic sex is more likely sterile in hybrids [Haldane, 1922].

## Final Remarks

Recent studies have provided new insight into avian W chromosome evolution, challenging the traditional view that the avian W is genetically inert and highly conserved across species. There are clear parallels with Y chromosomal evolution but also key differences, primarily regarding the relative importance of the W in reproduction and fertility traits. Recent technological advances offer new potential to resolve uncertainty over the functionality of the W, for instance by using single-cell RNA-seq to establish fine-scale expression patterns of Z- and W-linked genes through development and across species [e.g., Estermann et al., 2020] and targeted gene knockouts to test gene function [e.g., Ioannidis et al., 2021]. Therefore, the next couple of years hold much promise for disentangling the function and evolution of the W chromosome in birds.

## Conflict of Interest Statement

The authors declare no conflicts of interest.

## Funding Sources

This work was funded by a NERC Independent Research Fellowship to A.E.W. (NE/N013948/1) and a NERC ACCE DTP to P.D.P.

# Recovering the “Missing” Avian Genes Using Multi-Omics Data

(Prepared by Z.-T. Yin, J. Smith, and Z.-C. Hou)

Gene gain and loss are common events in the evolution of species, especially of birds, which have evolved many unique characteristics such as feathers, wings and flight capabilities, strong and lightweight skeletons, toothless beaks, high metabolic rates and heat absorption, sex, and unique respiratory and excretory systems [Kennedy and Vevers, 1976; Blomme et al., 2006]. The release of the first chicken genome provided the basis for systematic analysis of the similarities and differences between vertebrate and avian genomes [International Chicken Genome Sequencing Consortium, 2004]. In comparison

with other amniotes, bird genomes are more compact, and this difference may be related to the overall smaller cell size [Hughes and Hughes, 1995; Hughes and Friedman, 2008]. The reductions in genome size may be the result of the loss of noncoding DNA sequences, with bird genomes having less repetitive DNA, fewer pseudogenes, and shorter introns than mammalian genomes [International Chicken Genome Sequencing Consortium, 2004; Hughes and Piontkivska, 2005]. Importantly, the evolution of avian genomes also appears to involve the loss of protein-coding genes, as the total number of uniquely identified avian coding genes is much smaller than in other tetrapods (i.e., 23,294 in humans, GRCh38.p14; 19,404 in lizards, AnoCar2.0; 17,007 in chickens, GRCg7b). Paralog analysis revealed a higher overall incidence of gene families with fewer members in birds compared to other vertebrates [Hughes and Friedman, 2008]. Likewise, birds have a high rate of chromosomal rearrangements compared to other organisms, all of which may result in the deletion of protein-coding genes [Backström et al., 2010]. In recent years, the genomes of a large number of birds and lizards have been assembled and annotated, including zebra finches [Warren et al., 2010], chickens [International Chicken Genome Sequencing Consortium, 2004], turkeys [Dalloul et al., 2010], and duck [Zhu et al., 2021]. Moreover, large-scale bird genome projects [Jarvis et al., 2014; Zhang G et al., 2014], and chicken pan-genomes [Wang K et al., 2021; Li M et al., 2022] have also generated considerable genomic data. These large comparative genomic datasets identified hundreds of lost genomic-blocks in the bird genomes, and also suggested that hundreds of genes are missing in birds [Lovell et al., 2014; Zhang G et al., 2014].

The missing genes seem to be directly related to the unique physiological phenomena of birds. Several functionally important genes in mammals are supposed “missing” in chickens and have caused long-debated questions in bird biology. Spurious discovery of the missing/hidden genes in the bird genome has continued for decades. Previously, *BGN* [Blaschke et al., 1996], *CORO1A* [Xavier et al., 2008], *MAPK3* [Lemoine et al., 2009], *MMP14* [Simsa et al., 2007], *TBX6* [Lardelli, 2003; Ahn et al., 2012], *TSSK4* [Shang et al., 2013], and five adipokine genes [Dakovic et al., 2014] were reported to be missing in birds, however, several long-debated genes including *TNF-α*, and *leptin* have been cloned in birds [Prokop et al., 2014; Seroussi et al., 2016; Rohde et al., 2018]. This hide-and-seek game still continues, and does not appear to be ending anytime soon [Elleder and Kaspers, 2019].
Here we summarize recent efforts using multi-omics data to probe those genes missing/hidden in avian genomes.

## Reconstruction of Missing Genes in the Chicken Genome

While the hypothesis of missing genes in birds has been proposed for decades, researchers have found that some of the missing genes were, in fact, present in chickens or other birds. In the presence of large gaps and imperfect gene annotation in the genome, the de novo assembly of gene sequences using RNA-seq is considered to be an efficient way to identify unannotated genes in the genome. Attempts that only used a few tissues/organs have identified many missing genes in birds [Hron et al., 2015; Bornelov et al., 2017; Botero-Castro et al., 2017]. Recently, we used the raw data from 26 chicken tissues downloaded from the GenBank database to assemble and obtain 2,048,631 transcripts and identified 589 missing genes in birds [Yin et al., 2019b].

At the same time, the continuity and integrity of chicken genome assemblies have been rapidly improving. The chicken genome released in 2017 was assembled by third-generation sequencing technology, and the number of annotated genes increased significantly (2,768 noncoding and 1,911 protein-coding genes) [Warren et al., 2017]. In the Gallus_gallus-5.0 genome, 442 (77.41%, from a total of 571) genes thought to be missing in chickens [in Lovell et al., 2014, see Table S1, Table S6, plus select entries in Table S4 and Table S18] were annotated, indicating that there is no systematic deletion of genes in birds. With the development of sequencing and hybrid assembly technology, the genomes of different chicken breeds continue to be assembled and another 136 missing genes were further annotated in our recently assembled *Silkie* genome (unpublished). To date, it has now been shown that 528 (92.47%) genes that were thought to be missing, actually exist in chickens. This has been made possible by exploiting a large amount of multi-omics data available in chicken and has led to the revelation of genes with important functions such as *TNF-α* and *leptin* [Seroussi et al., 2016; Rohde et al., 2018]. Recent large-scale chicken pan-genome data have also identified thousands of genes that are not presented in the current chicken reference genome (Li M et al., 2022).

## Reconstruction of Missing Genes from Other Birds

In addition to chicken, researchers have reconstructed many genes thought to be missing from other birds. We collected data from various important tissues from duck

(24), pigeon (11), goose (8), and zebra finch (22) [Yin et al., 2019b], and an avian transcriptomic database containing a total of 9,296,247 transcripts was constructed by de novo transcriptome assembly. From this, we identified several genes in duck (583), pigeon (558), goose (537), and zebra finch (543) from 806 genes that were thought to be missing in birds [in Lovell et al., 2014, see Table S1; in Zhang G et al., 2014 see Table S10]. Only 135 genes were not found in this bird transcriptome database. The number of missing genes reconstructed in different birds by de novo assembly of large transcriptome data is similar, indicating that these genes thought to be missing exist across different bird species.

In recent years, duck functional genomics has developed rapidly. We have assembled the Mallard, Pekin duck, and Shaoxing laying duck genomes using a combination of third-generation sequencing, Bionano, and Hi-C sequencing technologies. These have proved to be a rich source of genetic information [Zhu et al., 2021]. In the Mallard duck the CAU_wild 1.0 genome has 1,872 more protein-coding genes annotated than the previous CAU 1.0 genome, including 89 genes previously thought to be missing in birds. Among these 89 genes, 5 genes have become pseudogenes, losing part of their gene function, 3 genes have been annotated as lncRNAs, and the remaining 81 genes remain as protein-coding genes. In addition, 240 genes were annotated as paralogous genes, and 108 genes had similar segments in the genome. Mining large multi-omics data assemblies and annotations now reveals that only 10 genes (from a total of 806 missing genes), to date, have not been reconstructed in birds, with the rest of the genes thought to be missing in birds having been shown to actually exist. The recovered gene list is shown in online supplementary Material 4.

## Development of New Methods to Identify More Missing Genes

Summarizing the characteristics of these reconstructed missing genes in birds and the reasons why they are thought to be missing can provide insights and methods for us to identify more missing genes. First, these reconstructed gene sequences have high GC content and length in many birds. The GC content of most of these “missing” genes is more than 60%, and few genes even have over 80% (the median GC content of the chicken genome is 42.22% and the median GC content of the duck genome is 41.99%) [Hron et al., 2015; Bornelov et al., 2017; Botero-Castro et al., 2017; Yin et al., 2019b]. At the same time, the multi-tissue transcriptome expression profiles of birds
showed that most of the reconstructed genes usually have strong tissue-specific expression. These genes are gener- ally expressed predominantly in one tissue and are rarely expressed in the other tissues [Yin ZT et al., 2019b]. High- throughput transcriptome-based assembly approaches have limitations for fully recovering missing genes due to technical factors such as the PCR amplification bias against GC-rich fragments [Beauclair et al., 2019]. Ex- pression patterns, i.e., tissue-specific expression patterns, and low expression, also limit the ability for full transcrip- tome assembly. Now, the third-generation sequencing technologies, which have less GC bias, such as single- molecule real-time (SMRT) and nanopore sequencing technologies, can obtain full-length transcripts directly, without assembly [Yin et al., 2019a; Kuo et al., 2020]. The missing genes will continue to be discovered with the ac- cumulation of full-length transcriptome data from more avian tissues from different physiological conditions.

Furthermore, the missing genes annotated in the chicken and duck genomes are mainly distributed on the microchromosomes, the ends of the chromosomes, and within regions showing a high content of tandem repeats clustering with non-canonical DNA structures [Zhu et al., 2021; Li M et al., 2022]. Long repetitive regions [Tre- angen and Salzberg, 2011], regions of high GC content [Chen et al., 2013], telomeric regions, fragmented micro- chromosomes [O'Connor et al., 2019], and adaptive as- sembly strategies have always proved problematic for en- abling complete bird genome assembly. To fully resolve the whole chicken gene sets, a telomere-to-telomere (T2T) genome is necessary. The recently completed hu- man T2T genome has now paved the way for the finished bird genome assembly [Miga et al., 2020; Hoyt et al., 2022; Mao and Zhang, 2022; Nurk et al., 2022]. Ultra-long ONT sequencing, high-precision HiFi sequencing data, multi- type auxiliary assembly data, and hybrid assembly using multiple strategies will greatly promote the quality of bird genome assembly [Sohn and Nam, 2018]. For large pres- ence/absence variations within species, we can enrich ge- nomic information by constructing high-quality multi- breed pan-genomes [Vernikos et al., 2015]. The Bird 10,000 Genomes (B10K) Project [Zhang et al., 2015] has generated insightful results and the future bird T2T ge- nome and pan-genome will undoubtedly reveal more genes. This complete gene map of birds will be critical for the further understanding of the biology and evolution of birds.

Finally, precise genome annotation will also provide
the necessary sequence and structural information for
mining more genes in birds. Annotation errors are un-

avoidable in genome annotation using automated pro-
cesses, especially for some protein-coding genes that can-
not be annotated in complex and high GC regions [Salz-
berg, 2019]. While applying full-length transcriptomic
data for genome annotation [Nudelman et al., 2018;
Wang X et al., 2019; Kuo et al., 2020], the use of novel an-
notation methods developed based on machine learning
can further improve the accuracy of annotation [Mahood
et al., 2020; Stiehler et al., 2020]. More accurate manual
annotation of important genome regions is also necessary
for novel gene identification [Dunn et al., 2019]. It can be
seen that, with the continuous development of omics
technology and analysis methods, the genome informa-
tion will be more complete, the annotation will be more
accurate, and the genes that were previously thought to
be missing in birds will continue to be discovered.

Statement of Ethics

All experiments with birds were performed under the guidance
of ethical regulations from the Animal Care and Use Committee
of China Agricultural University, Beijing, China.

Conflict of Interest Statement

The authors declare no competing interests.

Funding Sources

The work was supported by the National Waterfowl-Industry Technology Research System (CARS-42), the National Nature Science Foundation of China (31972525, 31572388), Beijing Joint Research Program for Germplasm Innovation and New Variety Breeding (G20220628007), Beijing Municipal Science & Technology Commission (Z211100004621007).

Data Availability Statement

Data have been submitted to the public databases under the fol-
lowing accession numbers: The raw data for assembling the tran-
scriptome database of chicken, duck, goose, pigeon, and zebra
finch were deposited in Sequence Read Archive (SRA) database
under the accession number SRP141084. The Mallard genome is
stored in NCBI under accession number PRJNA554956. The raw
data for the Silkie genome assembly can be found in the SRA da-
tabase under the accession number PRJNA805080 (unpublished
data).

Evolution of Protein-Coding and Long Noncoding
Genes of the Chicken Genome through the
Different Genome Assemblies and Their Associated
Annotations

(Prepared by F. Degalez, K. Muret, and S. Lagarrigue)

The chicken genome was the first avian genome se-
quenced because of its importance in human food pro-
duction, in fundamental biology like the study of devel-
opment or gene function conservation across evolution
[International Chicken Genome Sequencing Consor-
tium, 2004]. Since its first version in February 2004 (Gal-
gal2/WUGSC1.0), five new genome assemblies have been
released, each improving the genome sequences' accura-
cy. Along with these genome assemblies, numerous ge-
nome annotations were released, providing at least mod-
els for gene loci and transcripts supporting them. Since
the first annotated version (Ensembl v22 - May 25, 2004)
associated with the Galgal2 assembly, the number of
genes and the diversity of their biotypes have increased,
especially in 2015 with the introduction of long noncod-
ing RNAs (lncRNAs), which is concurrent with the first
initiatives of lncRNA annotation [Chodroff et al., 2010;
Necsulea et al., 2014; Li et al., 2015; Muret et al., 2017].

LncRNAs represent a large and heterogeneous class of
genes defined by transcripts longer than 200 nucleotides
without coding-potential capabilities [Derrien et al.,
2012]. They represent a variety of regulatory elements im-
plied in gene expression and can act at different levels by
using diverse biological mechanisms based on DNA,
RNA, or protein interactions [Guh et al., 2020]. As illus-
trated in Figure 10, lncRNAs can interact with DNA,
RNA, and proteins and act at different molecular levels:
nuclear organization (e.g., MALAT1 [Wang X et al.,
2021b]/NEAT1 [Yamazaki et al., 2018]) (Fig. 10A), ge-
nome integrity (e.g., TERRA [Barral and Déjardin, 2020])
(Fig. 10B), histone marks modification for silencing (e.g.,
Fendrr [Grote et al., 2013]) or activating (e.g., GATA3-
AS1 [Gibbons et al., 2018]) gene transcription (Fig. 10C),
loop formation to connect enhancers to promoter regions
(e.g., MYMLR [Kajino et al., 2019]) (Fig. 10D). LncRNAs
can modulate RNA splicing (e.g., linc-HELLP [van Dijk
et al., 2015]) (Fig. 10E), miRNA maturation (e.g., CCAT2
[Yu et al., 2017]/uc.372 [Guo et al., 2018]) (Fig. 10F), and
protein translation (e.g., BC1 [Wang et al., 2002]/MC-
M3AP-AS1 [Guo C et al., 2020]) (Fig. 10I) or their activ-
ity (e.g., NORAD [Munschauer et al., 2018]) (Fig. 10K).
They can also control the stabilization or the degradation
of molecules as miRNAs (e.g., ROR [Li C et al., 2017]/
DSCR8 [Wang Y et al., 2018b]) (Fig. 10G), mRNAs (e.g.,
PTB-AS [Zhu L et al., 2019]/TINCR [Xu et al., 2015])
(Fig. 10H), and proteins (e.g., PiHL [Deng et al., 2020]/
MALAT1 [Yan et al., 2016]) (Fig. 10J). LncRNAs can host
small ORFs [Choi et al., 2019] which code for peptides
(e.g., CASIMO1 [Polycarpou-Schwarz et al., 2018]/
DWORF [Nelson et al., 2016]) (Fig. 10M) or host in their
introns small RNAs [Sun Q et al., 2021] (e.g., MCM7
[Agranat-Tamir et al., 2014]/DLEU2 [Morenos et al.,
2014]) (Fig. 10N). They can control protein transfers

from cytoplasm to nucleus (e.g., NRON [Willingham et
al., 2005]) or from nucleus to cytoplasm (e.g., Discn
[Wang L et al., 2021]) (Fig. 10L). Finally, they can migrate
to other cells with exosomes (e.g., ZFAS1 [Pan et al.,
2017]/GAS5 [Chen et al., 2017]) (Fig. 10O).

Through their key roles in gene regulation, lncRNAs
are consequently involved in diverse biological and
pathophysiological processes [Ponting et al., 2009;
Muret et al., 2019; Gil and Ulitsky, 2020; Statello et al.,
2021]. Moreover, since most of the trait-associated
variations identified by genome-wide association stud-
ies (GWAS) concerned noncoding intervals of the ge-
nome [Manolio et al., 2009; Bouwman et al., 2018], this
reinforces the need to characterize the regulatory re-
gions of domesticated species such as lncRNA genes.
lncRNA genes have different characteristics compared
to protein-coding genes (PCGs). They are less expressed
[Derrien et al., 2012; Muret et al., 2017; Le Béguec et al.,
2018; Jehl et al., 2020], explaining why they have been
detected only recently – i.e., this last decade – by high-
throughput transcriptome sequencing technologies
(RNA-seq). Furthermore, lncRNA expression is more
specific to tissues, life stages, and conditions than that
of PCGs [Cabili et al., 2011; Derrien et al., 2012; Jehl et
al., 2020]. The identification of these genic entities is
therefore dependent on the variety of RNA-seq data
available to detect them.

After presenting the different chicken genome assem-
bles developed over the last 2 decades, we discuss the as-
sociated genome annotations provided by NCBI's RefSeq
and EMBL-EBI's Ensembl, the two reference annotation
databases. We characterize them in terms of number of
gene and transcript models, variety of biotypes, or in
terms of models that are shared by the two reference da-
tabases. We show that lncRNA loci are even less well-
known than PCG ones, although, for the latter, knowl-
edge of their transcripts can be further improved. Finally,
we discuss the impacts of these weaknesses and the value
of gathering different genome annotation resources, in
particular, for a better description of lncRNA loci, and
then present two initiatives. The MANE project yet lim-
ited to the human genome aims to synergize the NCBI's
RefSeq and EMBL-EBI's Ensembl "gene" databases to es-
tablish a consensus annotation. The second project, spe-
cifically realized for the chicken, is to provide a "gene"
database built from various resources including the NC-
BI's RefSeq and EMBL-EBI's Ensembl databases and oth-
er resources such as FAANG multi-tissue resources and
NONCODE database. This gene catalog is maintained at
each significant update in chicken genome assembly and
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Different mechanisms of lncRNA roles. Effects at the nuclear and telomere (A, B), transcriptional (C, D), post-transcriptional (E–H), translational (I), and post-translational levels (J–L). Role as small ORF host (M) and small noncoding RNA host (N). Implication in the exosome-mediated transfer (O). In purple, lncRNA;</figcaption></figure>

in blue, DNA; in green, other RNAs; in dark red, proteins. For more examples, other genes are presented in Muret et al. [2019], specifically genes involved in the regulation of lipid metabolism and their regulatory mechanisms.

genome annotation, the last version of June 2022, associated to the GRCg7b assembly, being composed of 23,926 PCGs and 44,428 lncRNA genes (available at http://www.fragencode.org).

## Evolution of the Reference Sequence of the Chicken Genome

As illustrated in Figure 11a, while the overall coverage of Galgal2/WUGSC1.0 was 6.63× [International Chicken Genome Sequencing Consortium, 2004], this parameter is
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Assembly versions associated with the chicken genome (a) and the number of publications associated with them (b). a RJF, Red Jungle fowl; W-L, White Leghorn; Cov., coverage; Scaf, scaffold; N50, scaffold N50. For more details, see online supplementary Material 5, Table 1.</figcaption></figure>

b Blue numbers, articles published using the corresponding assembly during the year 2020. Identification was made on PubMed Central by searching for the assembly name in different formats (e.g., GRCg6a or GRCg6 or Galgal6).

doubled for the two successive assemblies, i.e., Galgal3/WUGSC2.1 [NCBI RefSeq, 2006] and Galgal4 [NCBI RefSeq, 2011], which were released in November 2006 and 2011, respectively. Compared to the first version, these include 33 chromosomes (1–28; 32; W/Z; mitochondrial chromosome [MT]) but the number of scaffolds remained very high (~17,000), with 915 gaps between the scaffolds and a scaffold N50 which was quite low (~12 Mb) showing the incompleteness of the chicken genome sequence. Note that scaffold N50 is defined as the sequence length of the shortest scaffold at 50% of the total genome length. With the release of the Galgal5 assembly in December 2015 and with the improvement of sequencing technologies [NCBI RefSeq, 2015; Warren et al., 2017], the average coverage exploded and reached a global depth of 70×, leading to a better knowledge of the genome sequence. A few chromosomes were newly defined (addition of chromosomes 29–31). However, the quality of this genome sequence remained low with a high number of scaffolds (~24,000) and a lower scaffold N50 (~6.4 Mb) than before. These weak performances are likely due to long-read sequencing, which improved the detection of smaller scaffolds thus decreasing the N50 value [Warren et al., 2017].

In March 2018, a new assembly called GRCg6a was released by the Genome Reference Consortium, which has

taken the lead concerning the chicken genome assembly previously managed by the International Chicken Genome Consortium [NCBI RefSeq, 2018]. Tremendous progress — due to the addition of long read sequences, improved de novo assembly algorithms, manual annotation of contigs, and integration of finished BAC clone sequences — was made regarding the genome accuracy, as indicated by the drop in the number of scaffolds (from ~24,000 to ~500) with only 68 gaps between the scaffolds and an increase in the scaffold N50 (from ~6.4 M to ~20 M).

Furthermore, with the latest GRCg7 genome assemblies [NCBI RefSeq, 2021a, b], the knowledge of the chicken genome sequence improved even more in two main ways. First, the accuracy of the genome sequence increased due to improvements in sequencing and especially assembly technologies. The chicken genome is now composed of 42 chromosomes (1–39; W/Z; MT), reaching the number observed in the chicken karyotype. This assembly includes more microchromosomes with ~250 scaffolds with no gap between scaffolds, and the scaffold N50 reaching 90 Mb. Second, whereas a Red Jungle Fowl breed (known as RJF #256) was always used in previous assembly versions, a trio of chickens from diverse breeds was used for GRCg7. The new reference sequence was
generated from a female offspring from a cross between a broiler female and a white leghorn laying male. The RJF breed, considered to be the descendant of domestic chickens, was used as a good representation of broiler and layer chicken breeds; however, such a choice has a significant impact on the detection of variants leading to the identification of false positives. Since actual breeds have diverged from the RJF, a variation (e.g., SNP) at one position may be detected according to the RJF genome sequence, whereas this position was fixed in the population of interest. As an illustration, we have previously shown, using RNA-seq data from the liver of 11 different breeds (~750 RNA-seq of ~400 birds) aligned on the Galgal5 assembly, that the SNP number with reliable genotypes was on average 549,634 per population, but this number dropped to 339,539 (-38.2%) with a minor allele frequency ≥10% [Jehl et al., 2021]. This drop is mainly due to fixed variants in the populations since the number decreased to 438,837 (-20.2%) after only excluding the fixed variations.

Consequently, two genome assemblies were released in January 2021: GRCg7b representing the broiler breed and considered as the new genome reference, and GRCg7w representing the laying breed and considered as an alternative.

Because of the quite frequent change of the genome assembly compared to the time needed to conduct and publish a scientific study, a lot of works are published in outdated versions which can lead to the publication of misleading results and in disagreement with more recent versions (Fig. 11b). For example, in 2020, two years after the release of GRCg6a, 79 studies using this genome reference were published against 46, 41, and 5 published with the Galgal5, Galgal4, and Galgal3 versions, respectively. Some tools such as Liftoff [Shumate and Salzberg, 2020] or LiftOver [Kuhn et al., 2013] can be used to convert coordinates from one version to another. The first tool is based on the alignments of the gene features from one annotation to another, whereas the second tool is based on alignments of the best/longest syntenic regions for each region of the genome between assemblies (chain files). However, the use of these tools must be done with caution, especially for remote versions, because of important changes in the genome sequence.

## Evolution of the Two Reference Genome Annotations: A Breakthrough in 2015 with the Apparition of the lncRNA Gene Biotype

Genome annotation is not only evolving according to the version of the genome assembly but also to the evolu-

tion of annotation bioinformatics pipelines and data resources, mainly composed nowadays of RNA-seq data. In Figure 12a, genome annotations, from 2004 with the Galgal2/WUGSC1.0 assembly through 2022 with the GRC7b assembly, produced by the reference centers, NCBI's RefSeq and EMBL-EBI's Ensembl, have been analyzed. As illustrated, the gene number has increased, particularly due to the apparition in 2015 of lncRNAs, with 5,763 and 4,641 lncRNAs modeled by NCBI's RefSeq (v103) and EMBL-EBI's Ensembl (v94). This increase continues with the last genome annotation v107 (associated with the GRCg7b assembly) provided by EMBL-EBI's Ensembl with 11,944 lncRNAs compared with 5,504 for v106 (associated with GRCg6a). The number of PCGs remains constant at around 17,000 (see further for more explanation regarding such evolution).

In parallel to the gene number, it is important to make some comments about the transcript models that support these genes. As observed in Figure 12b, the transcript models can still be improved, as illustrated by the numerous changes observed between the two versions v105 and v106 of NCBI's RefSeq [NCBI RefSeq, 2022]. Only 6.4% of 93,980 transcripts identified in the 106 version are identical to those found in version v105. Such results can also be observed between genome annotations of different genome assemblies (not shown), or between genome annotations from the two reference centers, NCBI and EMBL-EBI, as illustrated in the next section.

## Differences between the Latest NCBI RefSeq and EMBL-EBI Ensembl Genome Annotations

For the same genome assembly, the two genome annotation bioinformatics centers, EMBL-EBI and NCBI, do not provide the same annotations, as illustrated in Figure 13. First, EMBL-EBI's Ensembl provides twice the number of lncRNA gene models compared to NCBI's RefSeq (shown in Fig. 13a) resulting in a total of 30,108 gene models (associated to 72,689 transcripts) including 17,007 PCGs, 11,944 lncRNAs, and 674 miRNAs compared with 25,638 gene models (associated with 85,704 transcripts) including 18,024 PCGs, 5,791 lncRNAs, and 799 miRNAs for RefSeq. These differences can be explained by the sample datasets and the annotation pipeline thresholds used specifically by the two bioinformatics centers. For example, NCBI's RefSeq does not consider lncRNAs supported by a single mono-exonic transcript in contrast to EMBL-EBI's Ensembl (with 1,157 lncRNA loci). Second, using the "GffCompare" software [Pertea and Pertea, 2020], we observed that most of the transcript models are different between the two ge-
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Gene numbers provided by NCBI's RefSeq and EMBL-EBI's Ensembl according to the genome annotation and genome assembly versions (a) and transcript model changes between two genome annotation versions from NCBI RefSeq for the same assembly – GRCg7b (b). a PCG, protein-coding gene; lncRNA: long noncoding RNA. b Comparison between versions 105 and 106 provided by NCBI [NCBI RefSeq, 2022]. Briefly, a score (between 0 and 1) for current and previous transcript features is calculated</figcaption></figure>

based on overlap in exon sequence and matches in exon boundaries. Pairs of current and previous features were categorized based on these scores and considering changes in attributes. New, new transcript models; Deprecated, transcripts removed or merged in the new version; Major changes, changes with great impact on the sequence or on the transcript attributes; Minor changes, minimal change ensuring similarity.

genome annotations, as shown in Figure 13b. Among the 72,579 transcripts from EMBL-EBI's Ensembl considered in the analysis, only 17.8% are strictly equal in the NCBI's RefSeq annotation. More than half (55.9%) are identified as new isoforms of an existing locus and 26.1% (18,922) transcripts are associated with 9,958 new gene loci resulting in more than one-third of the 30,108 gene models from EMBL-EBI's Ensembl not being known in NCBI's RefSeq.

Important differences exist between PCG and lncRNA transcripts. For PCG, most of transcripts from EMBL-EBI's Ensembl (70.8%) are new isoforms of the same gene loci existing in the two databases. These results show that the transcript isoforms are not well described with current RNA-seq resources. Indeed, most of RNA-seq data available in the public database are short-read RNA-seq; the long-read RNA-seq studies using the new technologies such as ONT or PacBio are still very limited [Kuo et al., 2017; Guan et al., 2022] due to the cost of these technologies and their low sequencing depth. For lncRNA, most of the lncRNA transcripts from EMBL-EBI's Ensembl (77.4%) are considered as new loci compared to NCBI's RefSeq. The main cause of this very low gene

overlap between the two genome annotations is the difficulty in capturing and therefore modeling lncRNAs compared to PCGs, due to specific features of lncRNA. First, lncRNAs are characterized by a global low expression; around less than 10% of the total reads of a sample analyzed by common technologies support lncRNA transcripts [Lagarrigue et al., 2022]. Second, they are tissue-, developmental stage-, and condition-specific [Cabili et al., 2011; Derrien et al., 2012; Jehl et al., 2020], conditions which are not covered by the limited number of RNA-seq samples used by the reference genome annotation centers compared to the tens of thousands of short-read RNA-seq generated by the avian scientific community which are available in the public database.

Moreover, the transcript models from NCBI's RefSeq are significantly longer than those of EMBL-EBI's Ensembl, as shown in Figure 13c, particularly for lncRNAs (almost twice the length), with nearly two supplementary exons by transcript (resp. a median of 5 vs. 10 exons/transcript, $p < 10^{-16}$) for both PCG and lncRNA models, with median exon sizes which remain similar (~250 bp). Moreover, NCBI's RefSeq provides a higher extreme distribution of transcripts per gene for PCGs compared to
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Features of the current NCBI RefSeq (v106) and EMBL-EBI Ensembl (v107) genome annotations based on the latest GRCh38 genome assembly. **a** Number of genes and transcripts according to gene biotypes for the two genome annotations. **b** The transcript models were compared between the two annotations according to 4 main classes (Equal isoform, New isoform, New loci, and Artifacts) according to the software “GffCompare” (op-</figcaption></figure>

tions: -S --no-merge) [Pertea and Pertea, 2020]. **c** #Tr/gn, number of transcripts per gene; #Ex/tr, number of exons per transcript; Tr. Size, transcript size considering only exonic regions; Ex. Size, exon size; In. Size, intron size. The median transcript sizes between RefSeq and Ensembl are 3,465 bp versus 2,317 bp, respectively, *p* < 10⁻¹⁶ (Wilcoxon rank sum test); for PCG 3,634 bp versus 2,870 bp, *p* < 10⁻¹⁶; for lncRNAs 2,952 bp versus 1,487 bp, *p* < 10⁻¹⁶.

EMBL-EBI's Ensembl (resp. 5 vs. 3 for the third quartile and resp. 9 and 5 for the last decile). Note that these numbers are far below what is described in human (resp. 6, 11, 18 transcripts per gene for the median, the third quartile, and the last decile (EMBL-EBI's Ensembl v107 with the GRCh38.13 assembly). This discrepancy can be explained in part by the variety of samples each reference used. EMBL-EBI's Ensembl combines a short-read RNA-seq dataset of 21 tissues from the Roslin Institute (1 stage-condition of a same breed per tissue, 21 samples of individual pool) with a short-read dataset of 7 tissues from the GENE-SWitCH project (3 stages of a same breed, 84 samples) and a long-read dataset of 6 tissues (7 samples) [EMBL EBI's Ensembl, 2022]. NCBI's RefSeq integrates data from various projects representing more than 20 tissues, different development stages and breeds for a total of 100 and 89 samples for short-read and long-read RNA-seq, respectively, in addition to “Cap Analysis Gene Expression” (CAGE) data including those from the FANTOM project [Lizio et al., 2017] for improving the annotation of transcription start sites [NCBI RefSeq, 2022].

For lncRNAs, the pattern in the distribution of the number of transcripts per gene is inverted between NCBI's RefSeq and EMBL-EBI's Ensembl with respectively 2 versus 3 for the third quartile and 3 and 5 for the last decile. Interestingly, these numbers are of the same order of magnitude in human (resp. 1, 2, and 5 transcripts per gene for the median, the third quartile, and the last decile), highlighting the general difficulty in capturing the transcript models associated with lncRNA genes.

## Interest in an Annotation Combining NCBI's RefSeq and EMBL-EBI's Ensembl

In summary, different genome annotations coexist with important differences in transcript models for PCGs and gene models for lncRNAs. Initiatives like the MANE project [Morales et al., 2022] for the human genome aim to synergize the NCBI's RefSeq and EMBL-EBI's Ensembl reference genome annotations to establish a consensus, although, so far, these efforts have focused only on PCGs. Such initiatives have yet to exist for livestock species, especially chicken. So far, most RNA-seq studies have ana-
lyzed gene expression and focus only on PCGs, using only one of these two reference annotations. As previously reported, the last two chicken reference genome annotations are quite similar in terms of PCG loci. Indeed, 18,024 and 17,007 PCG loci are respectively annotated for NCBI's RefSeq and EMBL-EBI's Ensembl; 15,711 (87.2%) loci from NCBI's RefSeq are shared with 15,848 (93.8%) loci from EMBL-EBI's Ensembl, even if most of the transcript models supporting these PCGs are different. However, these numbers drop for the 5,791 and 11,944 lncRNA loci respectively from NCBI's RefSeq and EMBL-EBI's Ensembl where 2,008 (34.7%) loci from NCBI's RefSeq are shared with only 2,118 (17.7%) loci from EMBL-EBI's Ensembl. Therefore, the use of only one of the two reference annotations enables the investigation of most PCG loci but can bias the study of lncRNA loci. Moreover, even when the expression is quantified at the gene level and not the transcript level, the high difference of transcript models previously reported — even for PCG loci — can have an impact. Thus, in the context of gene expression studies, results could differ depending on the annotation used [Zhao and Zhang, 2015]. Furthermore, the difference between transcript models, especially for PCGs, may have an important impact on variant prediction [McCarthy et al., 2014].

Concerning the recent studies interested in lncRNA gene expression, most of them have not used a reference genome annotation because of the very limited number of lncRNA loci represented in the versions – before the latest EMBL-EBI's Ensembl v107 – and produce a de novo annotation from investigators' own samples [for review, see Lagarrigue et al., 2022]. Such an approach is due to the recent democratization of RNA-seq data and the RNA-seq processing, gene modeling, and lncRNA prediction pipelines. However, these genome annotations are specific to one tissue or a set of tissues and characterized by their own gene identifiers, making result comparison difficult from one study to another. As reported in recent reviews [Kosinska-Selbi et al., 2020; Lagarrigue et al., 2022], the number of such publications has been constantly growing since 2015, with most of them focusing on the tissue-specific expression of lncRNAs or their differential expression in a given tissue between breeds or animal groups contrasted for an economically important trait in the species of interest. In most of these studies, a few lncRNAs have been highlighted as associated with the trait or tissue of interest whereas the lncRNA catalogues are not really exploited by the scientific community.

In parallel to these tissue-specific studies, a few multi-tissue studies have been performed in order to provide a

more comprehensive annotation of lncRNAs and considering their high tissue specificity. We can point two studies, that are part of the Functional Annotation of ANimal Genome consortium (FAANG) [Andersson et al., 2015], which have provided a multispecies lncRNA annotation: the first, Foissac et al. [2019], used 3 tissues of 4 female and male biological replicates of 4 farm species including the chicken updated to 12 tissues (personal communication); the second one, Kern et al. [2018], used 8 tissues of 2 female and male biological replicates of chicken, pig and cattle and was recently updated to 19 tissues for the chicken species [Guan et al., 2022]. Nevertheless, for a given species such as chicken, these studies remain limited due to the range of tissues, stage of development, condition that may exist.

In this context, we proposed since 2020 to provide a comprehensive gene catalog for chicken by gathering different resources, including EMBL-EBI's Ensembl, NCBI's RefSeq, and other multi-tissue databases, that we update at each important change of chicken genome assembly and annotation. Since the release of the new GRCg7b chicken genome sequence, we have recently updated the gene catalogue of 52,075 genes published in 2020 [Jehl et al., 2020], considering the last NCBI's RefSeq and EMBL-EBI's Ensembl annotations available in June 2022. First, we gathered the two genome annotation references, i.e., the v106 of RefSeq and the v107 of Ensembl resources. In addition to these two references, we chose to gather the two updated FAANG multi-tissue resources described above [Foissac et al., 2019; Guan et al., 2022], in which lncRNAs have been modeled in parallel with the PCG loci. The NONCODE resource composed only of lncRNA loci has also been used, even if this resource has not been updated since 2014 for the chicken [Zhao et al., 2021]. As a result, the EMBL-EBI's Ensembl and NCBI's RefSeq references grew respectively from 17,007 to 18,024 to 23,926 PCGs and from 11,944 to 5,791 to 44,428 lncRNA genes. This atlas associated to GRCg7b assembly is publicly available at http://www.fragencode.org [Degalez et al., in preparation] as the previous ones published in 2020 and associated to the Galgal5 and GRCg6a genome assemblies. In addition to the gene atlas (i.e., gtf file), a functional annotation of the genes across 40 tissues using different public resources is also provided as well as the lncRNA gene naming according to the official HUGO gene nomenclature committee (HGNC). Briefly, for the lncRNAs with an unknown function (frequent cases), the lncRNA adopts the symbol gene name of the gene harboring it, enriched by a suffix describing its genom-
<table><caption>**Table 2.** Chicken cell atlas: scRNA-seq data</caption><thead><tr><th>Focused tissue</th><th>Method<sup>a</sup></th><th>Reference</th><th>doi</th><th>GEO accession</th></tr></thead><tbody><tr><td>Embryo</td><td>SMART-seq, dissociation, FACS, LCM</td><td>Morrison et al. [2015]</td><td>10.1002/dvdy.24274</td><td>–</td></tr><tr><td>Neural crest</td><td>SMART-seq, dissociation, FACS</td><td>Morrison et al. [2017]</td><td>10.7554/eLife.28415</td><td>GSE108230</td></tr><tr><td>Primitive streak</td><td>SMART-seq Fluidigm C1</td><td>Vermillion et al. [2018]</td><td>10.1016/j.ydbio.2018.04.007</td><td>GSE89910</td></tr><tr><td>Developing limb</td><td>10x, dissociation</td><td>Feregrino et al. [2019]</td><td>10.1186/s12864-019-5802-2</td><td>GSE130439</td></tr><tr><td>Gonads</td><td>10x, dissociation</td><td>Estermann et al. [2020]</td><td>10.1016/j.celrep.2020.03.055</td><td>GSE143337</td></tr><tr><td>Neural crest (hindbrain)</td><td>10x, dissociation</td><td>Gandhi et al. [2020]</td><td>10.7554/eLife.57779</td><td>PRJNA624258</td></tr><tr><td>Retina</td><td>10x, dissociation</td><td>Hoang et al. [2020]</td><td>10.1126/science.abb8598</td><td>https://github.com/jiewwwang/Single-cell-retinal-regeneration</td></tr><tr><td>Skeletal muscle and fat (breast)</td><td>10x, dissociation</td><td>Li et al. [2020]</td><td>10.1186/s12864-020-07136-2</td><td>CRA002353</td></tr><tr><td>Developing retina</td><td>10x, dissociation</td><td>Tegla et al. [2020]</td><td>10.7554/eLife.54279</td><td>GSE142244</td></tr><tr><td>Limb muscle</td><td>10x, dissociation</td><td>De Lima et al. [2021]</td><td>10.1038/s41467-021-24157-x</td><td>GSE166981</td></tr><tr><td>Primitive streak</td><td>inDrops, dissociation</td><td>Guillot et al. [2021]</td><td>10.7554/eLife.64819</td><td>GSE161905</td></tr><tr><td>Cochlea</td><td>SMART-seq, dissociation</td><td>Janesick et al. [2021]</td><td>10.1016/j.celrep.2021.108900</td><td>https://umgear.org/dataset_explorer.html?search_terms = 33,761,346</td></tr><tr><td>Developing heart</td><td>10x, dissociation</td><td>Mantri et al. [2021]</td><td>10.1038/s41467-021-21892-z</td><td>GSE149457</td></tr><tr><td>Retina</td><td>10x, dissociation</td><td>Yamagata et al. [2021]</td><td>10.7554/eLife.63907</td><td>GSE159107</td></tr><tr><td>Pituitary</td><td>10x, dissociation</td><td>Zhang et al. [2021]</td><td>10.3389/fphys.2021.562817</td><td>CRA003604</td></tr><tr><td>Developing hypothalamus</td><td>10x, methanol</td><td>Kim et al. [2022]</td><td>10.1016/j.celrep.2021.110251</td><td>GSE171649</td></tr><tr><td>Lymphocyte</td><td>10x, suspension</td><td>Qu et al. [2022]</td><td>10.3389/fmicb.2022.800618</td><td>PRJNA687808</td></tr><tr><td>Germ cell</td><td>10x, dissociation, FACS</td><td>Rengaraj et al. [2022]</td><td>10.1016/j.csbj.2022.03.040</td><td>PRJNA761874</td></tr><tr><td>Neural crest</td><td>10x, dissociation</td><td>Williams et al. [2022]</td><td>10.7554/eLife.74464</td><td>GSE131688, GSE181577</td></tr></tbody></table>

Tabula Gallus) has been proposed to create a cell atlas of all tissues in the mature and developing chicken [Yamagata, 2022]. Like several ongoing cell atlas projects (see below), the chicken project will collect scRNA-seq data for chickens, characterize each cell type, and eventually make available information that integrates diverse modalities.

The chicken cell atlas is still in its infancy (Table 2) [Yamagata, 2022; Liu Y et al., 2022]. Nonetheless, a couple of pioneering studies have realized this endeavor, re-

vealing various cell types and their states in the chick limb buds [Feregrino et al., 2019; De Lima et al., 2021], the early primitive streak stage [Vermillion et al., 2018; Guillot et al., 2021], the neural crests [Morrison et al., 2017; Gandhi et al., 2020; Williams et al., 2022], and the neural retina [Hoang et al., 2020; Tegla et al., 2020; Yamagata et al., 2021]. Most studies have used embryonic or juvenile tissues, likely due to accessibility. It is practicable because hatched chicks are highly active and generally considered mature. It takes 21 days, on average, for an egg to hatch

<footnote><sup>a</sup> 10x, 10x chromium system from 10x Genomics (Pleasanton, CA, USA). Dissociation indicates enzymatic dissociation.</footnote>
once incubation begins. However, embryos are often not consistently developed, thus staged according to the Hamburger and Hamilton series [Hamburger and Hamilton, 1951]. Interestingly, in single-cell analysis, one embryo consists of a series of cell types and their states covering different developmental stages. Therefore, making multiple developmental atlases at close time points is not essential. In contrast, auxiliary atlases must be generated reflecting different factors such as sex [Clinton et al., 2001] and variable genetic background [Núñez-León et al., 2019].

The raw data from and references to scRNA-seq studies are searchable at NCBI's GEO database (https://www.ncbi.nlm.nih.gov/geo/) and several single-cell reference sites (e.g., https://panglaodb.se/papers.html, https://www.nxn.se/single-cell-studies/). To explore the original datasets, some interactive viewers for single-cell data are available (Single Cell Portal, https://singlecell.broadinstitute.org; UCSC Cell Browser, https://cells.ucsc.edu; EMBL-EBI, https://www.ebi.ac.uk/gxa/sc/home). Thus, the first endeavor toward the chick cell atlas project will be to create an armamentarium, incorporate multimodal data, display those datasets and assist users. Tabula Muris, Tabula Sapiens/Human Cell Atlas (HCA)/HuBMAP, Tabula Drosophilae, and C. elegans atlas are examples of other species [Tabula Muris Consortium, 2018; HuBMAP Consortium, 2019; Haniffa et al., 2021; Taylor et al., 2021; Lindeboom et al., 2021; Eraslan et al., 2022; Li H et al., 2022; Tabula Sapiens Consortium, 2022]. Species-specific community websites like GEISHA (http://geisha.ari-zona.edu/geisha) and Chickspress (https://geneatlas.arl.arizona.edu/) may provide vital starting points as a resource for chicken as in the cases for other model species. However, in the coming years, emerging crypto-technologies such as peer-to-peer data storage and smart contract protocols could transform data sharing methods (see below).

## Next Steps: Multimodal, Spatial, and Temporal Atlases

In addition to scRNA-seq and relevant single-nucleus RNA sequencing, other multimodal single-cell technologies, which simultaneously profile multiple data types in the same cell, represent a frontier for discovering new cell types and characterizing cell states [Stuart and Satija, 2019]. The additional modalities include epigenome, proteome, glycome, metabolome, electrophysiology, morphology, and connectome (Guo S et al., 2021; Lee et al., 2021; Saunders et al., 2021; Sun YC et al., 2021; Mund et al., 2022].

Among those modalities, a series of single-cell sequencing methods for detecting heritable DNA methylation and altered chromatin configurations allow the description of epigenetic changes on a genome-wide scale. In particular, a single-cell sequencing assay for transposase-accessible chromatin (scATAC-seq) is the most commonly used method for studying epigenetic landscapes in single cells [Stuart and Satija, 2019; Armand et al., 2021]. Although scRNA-seq and scATAC-seq are different, both represent the activity of genes. Thus, it is also vital to analyze proteomics and metabolomics to understand actual cell function.

The data gained from scRNA-seq and other dissociated protocols lead to the loss of spatial information. By contrast, spatial biology and spatial transcriptomics include histological, cellular, and subcellular information of transcripts in spatial context and coordinates [Rao et al., 2021; Zhuang, 2021; Palla et al., 2022; Moffitt et al., 2022; Chen A et al., 2022; Moses and Pachter, 2022]. The type and function of cells are further designated by cell morphology and cell interactions, including neuronal connectivity. eCHIKIN (electroporation- and CRISPR-mediated Homology-Instructed Knock-IN) is a technique for CRISPR-mediated genome editing in somatic cells to insert GFP or Cre cDNA into genes identified as cell-type specific in scRNA-seq data [Yamagata and Sanes, 2021]. This technique will reveal cell morphology and connectivity and potentially describe the molecular and spatial networks that organize the proteome [Cho et al., 2022]. These post-transcriptional modalities, together with spatial and temporal information, facilitate the resolution of cell types and states and provide more critical information on cell function.

## Open Science and DeSci

In order to promote any scientific research, it is a prerequisite to have a supportive community, raise funds, and build facilities. Furthermore, all scientific data should be made openly accessible and maintained permanently. Nonetheless, most of the current “big science” projects have suffered from several drawbacks. For example, not all contributions and data submitted by individuals or institutions are satisfactorily credited. Instead, the management is often exceedingly centralized: only a handful of influential scientists are highly recognized as leaders of successful projects.

Decentralized science (DeSci) is an emerging movement that proposes to build a shared infrastructure for disseminating, assessing, funding, crediting, and storing data and knowledge using blockchain technologies
(https://ethereum.org/en/desci/) [Hamburg, 2021]. It aims to create an ecosystem where all the researchers are motivated to share their data and receive credit for their effort while allowing everyone open access to research materials using crypto-technologies such as non-fungible tokens (NFTs). More important, scientific organizations can be governed using tokenized incentive structures without powerful leaders by establishing decentralized autonomous organizations (DAOs). DAOs can provide more flexible and agile funds using retroactive or quadratic funding by working together with a consortium of academic, philanthropic, and corporate partners. The scientific data and achievements can be owned and credited using research NFTs (rNFTs) or intellectual property NFTs (ipNFTs). Peer-to-peer data storage such as Interplanetary File System (IPFS) warrants storing and distributing data enduringly. These cutting-edge crypto-technologies will be able to support a gamut of endeavors in various basic science projects, create an armamentarium, and organize research resources and incentives. Thus, I wish to propose an international collaboration among many researchers by establishing galusDAO to facilitate this chicken cell atlas project in a decentralized manner.

### Conflict of Interest Statement

There are no conflicts of interest.

## Standardized Chicken Gene Nomenclature to Support Functional and Comparative Studies

(Prepared by F.M. McCarthy, S. Davey, M.C. Young, K.C. Potter, A. Lanke, P. Barela, M.C. Casono, M. Chiodi, L. Cigan, N. Das, A. Goodell, S.M. Johnson, J.H. Keroack, and D. Webb)

Standardized gene nomenclature facilitates unambiguous communication about genes and allows accurate indexing of associated scientific literature. Genes are assigned a full-length name that succinctly describes what is known about its function and a short gene symbol that is unique, the latter being the nomenclature most often used in scientific communications. In addition, we annotate additional names used in published literature to support literature indexing (referred to as synonyms or aliases). Gene nomenclature efforts were developed for representative species of each of the vertebrate lineages [Burt et al., 2009; Kusumi et al., 2011; Fortriede et al., 2020; Tweedie et al., 2021; Bradford et al., 2022], and more recently these groups coordinate their efforts to ensure

that gene nomenclature is consistent across vertebrate species. This consistency expedites comparative studies amongst vertebrate species and enables discoveries about gene evolution in these animals. With support from public resources such as NCBI, standardized nomenclature can be propagated from a representative species such as chick to other species in the same taxonomical class (i.e., birds), further promoting comparative studies and revealing genetic differences between species.

The Chicken Gene Nomenclature Consortium (CGNC) was convened in 2009 to develop and promote standardized gene nomenclature for chicken genes [Burt et al., 2009]. The CGNC provides updates following each major annotation update of the chicken genome and has an ongoing manual annotation effort. Nomenclature provided by CGNC is displayed and continually updated in NCBI's Gene Database. CGNC cross-references both NCBI Gene IDs (formerly Entrez Gene IDs) and Ensembl Gene IDs. Nomenclature may be assigned automatically via multiple orthology and homology searches [Eyre et al., 2007] or by manually reviewing gene location, syn teny, and published papers. While we make every effort to assign standardized gene symbols and names that are based upon human nomenclature when there is a clear orthology between genes, there are exceptions to this rule. For example, exceptions are made in instances where the chicken gene names are well established in the literature (e.g., ovalbumin) or when human nomenclature refers to genomic features or physiology not common to birds (e.g., when gene names reference the X or Y chromosomes, the HLA region, or human blood groups). Close collaboration between vertebrate gene nomenclature groups allows curators to come to a consensus and ensure that gene nomenclature for orthologs conserved across a diverse range of vertebrate species can be practically applied to all vertebrate lineages.

The last major update of CGNC was completed in May 2022, and this update coincides with updated gene annotations for the GRCg7 assembly. This genome assembly represents a distinct change from other chicken genome assemblies, as it is based upon modern chicken lines rather than Red Jungle Fowl. Since both broiler and leghorn lines were sequenced, CGNC now provides information about genes from both the maternal broiler assembly and the paternal white leghorn layer. There are currently 22,315 genes automatically assigned nomenclature and 3,602 manually approved genes. While initial manual curation efforts focused on assigning nomenclature for 1:1 orthologs between chicken and human, more recently we are focused on assigning standardized names for genes
that are not found in mammals. This latter gene set includes gene family expansions found in chickens or birds, or genes found in other vertebrate lineages but lost or significantly altered in mammals. Several of these projects have only been made possible with assistance from community experts, and some examples of these are described below.

## MHC-Y Region Genes

The chicken major histocompatibility complex (MHC) contains two independent regions which segregate independently, MHC-B and MHC-Y [Briles et al., 1993; Miller et al., 1994]. The unique nature of MHC regions means that these genes need to be manually reviewed, both to identify the genes and to assign nomenclature in a standardized manner. With assistance from Dr. Miller, we reviewed and provided standardized gene nomenclature for genes associated with the MHC-Y locus. This includes a set of 107 genes and encompasses class I and II molecules, c-type lectin molecules, zinc finger proteins, and leukocyte receptor cluster members [Goto et al., 2022]. Class I genes are denoted with the “MHCY” prefix for gene symbols and the full gene name has the format “major histocompatibility complex Y, class I heavy chain”, with individual genes numbered sequentially and pseudogenes indicated with the designation “P” in the symbol and “pseudogene” for the full gene name. Likewise, class II genes are named using the “MHCY2B” gene symbol prefix and standardized gene names starting with “major histocompatibility complex Y, class II beta”. The MHC-Y region includes eight leukocyte receptor cluster genes that are most closely related to the human leukocyte receptor cluster member 9 family (LENG9; HGNC:16306). To clearly indicate this similarity while ensuring that these genes are not mistaken for direct homologs, the chicken leukocyte receptor genes are designated as “leukocyte receptor cluster member 9 like 3, MHCY region” and gene symbols have the prefix “LENG9L.” Likewise, c-type lectin gene names include the information that these genes are located in the MHC-Y region. We note that 40 of these genes are not annotated on current reference assemblies and may be breed-specific; CGNC has reserved these gene names for future use. This systematic review of MHC-Y genes provides a clear and unambiguous naming system that can be extended to genes found in other chicken lines. We expect to extend this work to provide standard gene nomenclature for genes in the MHC-B regions and the nucleolus organizer region (NOR), although this effort may require assistance from other community experts.

## Chicken Immunoglobulin Receptor Genes

Chicken microchromosome 31 contains more than 100 chicken Ig-like receptors, which function as activating (CHIRA), inhibitory (CHIRB), or bifunctional receptors (CHIRAB) [Nikolaidis et al., 2005; Laun et al., 2006]. This nomenclature represents a challenge because, although it is well established in scientific literature (including in closely related species) [Windau et al. 2013], it is desirable to remove references to a specific species to promote transfer across multiple species. To resolve this dichotomy, we propose to retain the “CHIR” symbol for these genes but adjust the gene name to “cluster homolog of immunoglobulin like receptor” and to add appropriate “chicken Ig-like receptor” synonyms used in publications to support searching and indexing. We are currently working with Ms Brandi Sparling in Dr. Drechsler’s laboratory to ensure that these chicken genes are correctly identified as activating (CHIRA), inhibitory (CHIRB), or bifunctional (CHIRAB) receptors and named systematically.

## Heat Shock Proteins

Heat shock proteins are an example of functionally related proteins that are expressed under stress conditions [Nakai and Ishikawa, 2001]. Multiple lines of independent research are actively studying stress responses in chickens [Olanrewaju et al., 2010; Zhao et al., 2011; Jastrebski et al. 2017; Sarsour and Persia, 2022], and the GRCg7 assemblies have 16 genes currently named with a temporary “LOC” designator that are assigned as “heat shock transcription factor, X-linked-like.” Using gene names and literature searches we identified a more comprehensive group of 55 chicken genes associated with heat shock responses. These genes were manually reviewed, and corrected gene nomenclature assigned where required. This process included assigning new names for the expanded gene set of heat shock transcription factors which had formerly been associated with mammalian X and Y chromosomes (HSFX1-HSFX4 and HSFY1). This project highlights the importance of manual review for gene nomenclature automatically assigned across taxonomic groups, and we are currently reviewing additional gene sets with names that reference X and Y chromosomes.

Manually assigning gene nomenclature requires that biocurators have a good working knowledge of how public databases store, integrate, and share information, as well as understanding key principles in genetics, molecular and cellular biology, evolutionary biology, and physiology. CGNC offers research experiences for undergraduate and high school students who wish to complete a
research project. Students are introduced to the concept of standardized gene nomenclature and why it is important. They work collaboratively to review automatically assigned gene nomenclature, including reviewing associated publications and completing their own phylogenetic analyses. Students review each other's work, and all nomenclature is quality checked before data entry to the CGNC database. This allows students the opportunity to see the results of their research reflected in databases which use CGNC information.

The current focus for manual curation efforts is to assign nomenclature for gene families and chicken genes that are currently designated with the “LOC” prefix in the NCBI gene set. Gene families frequently have contractions and expansions between species and homology searches often cannot identify orthologs within a family. Many of these LOC genes are similar to (or “-like”) genes found in other vertebrates but have no clear ortholog and are either new members of a gene family or represent lineage-specific genes. For example, avian gene families with notable changes compared to mammals include keratins, histones, ribosomal RNAs, cadherins, and olfactory receptors, and many of these gene families are included in the current LOC genes. Moreover, this process is complicated by the removal and merging of LOC genes during the process of assembly and annotation updates. However, the manual review of these genes necessarily includes closer inspection of annotation and assembly and can provide useful information about regions of the genome that can be improved.

### Acknowledgments

D.W. and F.M.M.: We wish to acknowledge Dr. Marcia Miller (City of Hope), Ms Brandi Sparling (Western University of Health Sciences), and Dr. Carl Schmidt (University of Delaware) for their invaluable expertise in MHC-Y region, chicken immunoglobulin receptors, and heat shock genes, respectively. We also acknowledge and sincerely thank gene nomenclature biocurators at HGNC, the Vertebrate Gene Nomenclature Committee (VGNC), Xenbase and Zfin for their ongoing feedback and support.

### Conflict of Interest Statement

The authors have no conflicts of interest to declare.

### Funding Sources

CGNC was developed and supported by National Institutes of Health (R24 GM079326) and US Department of Agriculture National Institute of Food and Agriculture (competitive grant no. 2011-67015-30332 to F.M.M.). F.M.M is partly supported by USDA NIFA Hatch funding from National Research Support Project-8 (NRARZT-1370570-R50-131) and Multistate Research Project NC1170 (ARZT-1370520-R50-130).

## Distinct Hypothalamus and Breast Muscle Transcriptomic Response to Heat Stress under Newcastle Disease Virus Infection

(Prepared by Y. Wang, P. Saelao, C. Kern, B. Zhao, R.A. Gallardo, T. Kelly, J.M. Dekkers, S.J. Lamont, and H. Zhou)

Global warming has had adverse effects on mammals, birds, and plants [Bellard et al., 2012]. Physiological and behavioral disorders resulting from heat stress raise concerns for poultry production, health, and welfare. Decreased feed efficiency, growth rate, and egg production, as well as increased susceptibility to disease, are major challenges associated with heat stress in the poultry industry. Total economic losses due to heat stress amount to approximately 1.69–2.36 billion USD annually for the US livestock production industry, of which 128–165 million USD is lost from the poultry industry [St-Pierre et al., 2003]. Currently, in genetic selection and breeding programs, high-performance broilers and layers are superior in productivity; however, they cannot maintain their production in the presence of heat stress compared to non-selected birds [Kumar et al., 2013]. In addition, heat stress suppresses the chicken’s immune response, which increases the bird’s susceptibility to infectious disease threats [Monson et al., 2018]. Therefore, it is essential to better understand the deleterious effects and underlying mechanisms of heat stress on immunity and growth performance.

The hypothalamus in the host plays a critical role as a central regulator of temperature in chickens. It links the nervous system to the endocrine system and primarily mediates thermoregulation, food intake, and stress response [Chen et al., 2015]. It functions as a central hub to interact and coordinate downstream metabolism [Bohler et al., 2021]. Heat stress activates the hypothalamic-pituitary-adrenal (HPA) axis by initiating an appropriate response in the preoptic area of the hypothalamus [Tzschentke and Basta, 2002], which simulates the secretion of gamma-aminobutyric acid and then inhibits the dorsomedial nucleus (DMN) of the hypothalamus. As a result of this inhibition, circuits to muscle, fat, and the cardiovascular system are deactivated reducing heat production and increasing heat loss [Molinas et al., 2019]. Activation of the HPA axis results in the release of both corticotrophin releasing hormone (CRH) and arginine vasopressin (AVP), and secretion of adrenocorticotropic hormone (ACTH) by the pituitary gland, which stimulates the release of glucocorticoids from the adrenal cortex. Circulating glucocorticoids interact with a variety of
cells to regulate both metabolism and immune function [Chen et al., 2015]. Disruption of the metabolic and immune system function in birds gives rise to lower productivity, such as low breast muscle yield, and an increased risk of disease and mortality.

Genetics plays an important role in the host response to heat in chickens [Wolc et al., 2019]. Therefore, genetic selection offers a feasible and sustainable option to improve heat stress resistance and immune response in chickens. Genetic selection of chickens with improved resilience to heat stress requires a better understanding of the genetic contribution and molecular mechanisms of the heat stress response. Previous studies conducted by our group have characterized transcriptome profiles of several chicken tissues under heat stress. Candidate genes differentially expressed in the comparisons of heat-stressed and control groups in chicken lungs, tracheas, Harderian gland, and livers were identified as a first step toward identifying targets for validation to serve as selection markers [Saelao et al., 2018, 2021; Wang Y et al., 2020]. A series of genome-wide association studies (GWAS) have identified several quantitative trait loci (QTL) associated with heat stress response such as body temperature, growth phenotypes such as body weight and breast muscle yield, and immune-related phenotypic traits such as mortality, viral replication, and antibody levels, highlighting the genetic contribution to heat stress tolerance in different chicken populations [Van Goor et al., 2015; Rowland et al., 2018; Saelao et al., 2019; Wolc et al., 2019]. However, underlying mechanisms of genetic resistance and susceptibility to heat stress are still not fully understood. Physiological parameters from two genetically distinct, highly inbred lines (Fayoumi and Leghorn) that were heat stressed under NDV infection, were reported [Wang Y et al., 2018a]. The relatively more heat tolerant Fayoumi birds showed distinct responses during acute and chronic heat stress compared to the relatively heat stress susceptible Leghorn birds [Wang Y et al., 2018a]. Better hemostasis was maintained in the Fayoumi birds under both abiotic and biotic stressors. Acid-base balance in the host was one of the key factors accounting for the genetic differences between these two lines [Wang Y et al., 2018a]. The liver transcriptome profile of these two inbred lines revealed that the relative heat tolerant and disease resistant Fayoumi line activated not only metabolic but also immune regulation with heat stress and viral infection, and the susceptible Leghorn birds were only able to maintain the basic metabolic responses [Wang Y et al., 2020].

Thermoregulation is mediated by both the nervous and endocrine systems with a key regulator of the nervous system being the hypothalamus. A major downstream organ, the breast muscle, deserves intensive study to help better understand different aspects of thermoregulation. Therefore, the current study was designed to survey the transcriptomic profiles of the hypothalamus and the breast muscle of these two inbred lines to characterize the global gene expression response to heat stress in distinct genetic backgrounds. Identifying candidate genes and pathways associated with heat tolerance in a variety of tissues will aid discovery of important molecular markers to target in genetic selection to improve heat resilience in poultry.

## Experimental Populations

Two genetically distinct, highly inbred lines, Fayoumi (M 15.2) and Leghorn (GHs 6) having inbreeding coefficient of 99.95% [Zhou and Lamont, 1999], from the Iowa State University Poultry Teaching and Research Farm (Ames, IA) were used in the current study. Fifty-five Leghorns and 56 Fayoumi birds were housed in two temperature and humidity-controlled isolators. Birds were provided with ad libitum access to food and water. On day (d) 1 of age, 30 Leghorn and 31 Fayoumi birds were randomly chosen as the treatment groups and housed in one isolator and the rest of the birds were used as the non-treatment groups in another isolator. The two genetic lines were mixed in each isolator. From 14 days of age to the end of the experiment (41 days of age), the heat-treated groups were exposed to continuous heat stress of 38°C for the first 4 h and then decreased to 35°C, while the non-treatment groups were maintained at 29.4°C for the first week and then 25°C throughout the whole experiment. On d21, birds in the heat-treatment groups were inoculated with $10^7$ EID$_{50}$ (one EID$_{50}$ unit is the amount of virus that will infect 50% of inoculated embryos) Newcastle disease virus (NDV) La Sota strain through both eyes and nares (50 µL per each eye and nostril). The non-treatment birds were given 200 µL phosphate-buffered saline (PBS) via the same routes. The animal experiment was performed according to the guidelines approved by the Institutional Animal Care and Use Committee, University of California, Davis (IACUC #17853).

## Blood Parameter Measurements

Physiological blood parameters were measured at three stages: acute heat (AH) at 4 h, chronic heat (CH) at 7 days, and chronic heat combined with NDV infection
(CH&NDV) at 10 days post-initiation of heat treatment. Thirteen parameters including four chemistry/electrolyte parameters (concentrations of sodium [Na⁺], potassium [K⁺], ionized calcium [iCa²⁺], and glucose [Glu]); seven blood gas parameters (blood pH, carbon dioxide partial pressure [PCO₂], oxygen partial pressure [PO₂], total carbon dioxide [TCO₂], bicarbonate [HCO₃], base excess [BE], and oxygen situation [sO₂]) were measured immediately by using an i-STAT Portable Blood Analyzer as we described in the previous study [Wang Y et al., 2018a].

## Tissue Sample Collection and RNA Isolation

A total of 32 chickens (4 birds per line per treatment at d14 [4-h post-heat stress treatment, acute phase (AH)] and d23 [9 days post-treatment, chronic phase (CH)]) were randomly selected from treatment and non-treatment Leghorn and Fayoumi birds. The birds were euthanized, and the hypothalamus and the breast muscle samples were collected for RNA isolation. Total RNA was isolated using Trizol (Invitrogen, Carlsbad, CA) according to the manufacturer's protocol. DNase I (Ambion, Austin, TX) digestion was carried out after RNA isolation, and the RNA concentration and purity were determined by measuring absorbance at 260 nm and A260/A280 ratio using a NanoDrop ND-1000 spectrophotometer (Nano-Drop Technologies, Wilmington, DE), and RNA quality was checked by Agilent Bioanalyzer (Agilent, Santa Clara, CA). The RNA samples were stored at -80°C until further use.

## RNA-Seq Library Preparation and Data Analysis

For each sample, 500 ng of total RNA was used to construct a cDNA library by using the NEBNext® Ultra™ RNA library prep Kit for Illumina® (New England Biolabs, Ipswich, MA). In total, there were 64 RNA-seq libraries. The cDNA libraries were quantified by Qubit (ThermoFisher Scientific, Waltham, MA) and validated by Agilent Bioanalyzer High Sensitivity DNA Assay (Agilent, Santa Clara, CA) and then sequenced on the HiSeq4000 platform (Illumina, San Diego, CA) for 100-bp paired-end reads (DNA Core Facility, University of California, Davis, CA). Sequencing data can be accessed at GEO (PRJNA896699).

Raw reads were trimmed to remove adapter sequences and low-quality bases were removed using the Trim Galore program (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). The RNA-seq analysis was carried out by using the same method as described in the liver transcriptome study [Wang Y et al., 2020]. The statistical model included the effects of line and condition

for each treatment in each tissue, along with the interactions between condition and line. DEGs were declared with a false discovery rate (FDR) <0.05. Gene expression profiles were compared between each pair of the 8 groups in each tissue: Leghorn non-treated (LC) and Leghorn treated (LT) with AH and CH&NDV; Fayoumi non-treated (FC) and Fayoumi treated (FT) with AH and CH&NDV. In each tissue, the comparison groups with the AH and CH&NDV treatment were FTFC (Fayoumi treatment vs. Fayoumi non-treatment), LTLC (Leghorn treatment vs. Leghorn non-treatment), FCLC (Fayoumi non-treatment vs. Leghorn non-treatment), and FTLT (Fayoumi treatment vs. Leghorn treatment).

## Gene Ontology

Statistics related to the overrepresentation of functional categories were performed using DAVID [Huang et al., 2009; Sherman et al., 2022]. A fold enrichment >2 and FDR <20% was considered significant.

## Pathway Analysis

Pathway analysis using the DEGs of within-line contrasts was performed by using the Ingenuity Pathway Analysis software (IPA; Qiagen, Redwood City, CA, USA) [Kramer et al., 2014]. Significant associations (*p* < 0.05) and a Z-score cutoff of |z| > 1.64 were used to identify significantly activated or inhibited pathways.

## Gene Co-Expression Network Analysis

The Weighted Gene Co-expression Network Analysis (WGCNA) package in R was used for gene co-expression network analysis [Langfelder and Horvath, 2008, 2012]. A soft threshold of 13 was set for generating an adjacency matrix based on co-expression, and the minimum module size was arbitrarily set at 30. To evaluate associations of co-expressed gene clusters with line and treatment, the Leghorn and Fayoumi lines were given nominal values of 1 and 2 and non-treatment and treatment the nominal values of 0 and 1. For the continuous traits collected from the previous physiological study [Wang Y et al., 2018a], association of co-expressed gene clusters with the continuous traits pH, PCO₂, HCO₃, BE, PO₂, sO₂, Glu, Na⁺, K⁺, and iCa²⁺ were also evaluated. The driver genes were identified by high absolute values of gene significance (GS > 0.5) and module membership (MM > 0.5) with a threshold of *p* < 0.05 [Horvath and Dong, 2008].

## Summary of RNA-Seq Analysis

Sixty-four chicken cDNA libraries were prepared from hypothalamus and breast muscle samples and sequenced
<table><caption>**Table 3.** Summary statistics of RNA-Seq output</caption><thead><tr><th></th><th>Line</th><th>Treatment</th><th>Time point</th><th>Raw reads</th><th>Aligned reads</th><th>Alignment rate, %</th></tr></thead><tbody><tr><td rowspan="8">Hypothalamus</td><td>Leghorn</td><td>Non-treated</td><td>Acute heat</td><td>90,131,605</td><td>82,596,745</td><td>91.91</td></tr><tr><td>Leghorn</td><td>Treated</td><td>Acute heat</td><td>99,594,536</td><td>92,251,909</td><td>92.91</td></tr><tr><td>Fayoumi</td><td>Non-treated</td><td>Acute heat</td><td>88,123,256</td><td>81,628,366</td><td>92.85</td></tr><tr><td>Fayoumi</td><td>Treated</td><td>Acute heat</td><td>88,021,678</td><td>81,552,876</td><td>92.88</td></tr><tr><td>Leghorn</td><td>Non-treated</td><td>CH&amp;NDV</td><td>85,586,252</td><td>78,499,724</td><td>93.38</td></tr><tr><td>Leghorn</td><td>Treated</td><td>CH&amp;NDV</td><td>82,230,451</td><td>75,662,675</td><td>93.28</td></tr><tr><td>Fayoumi</td><td>Non-treated</td><td>CH&amp;NDV</td><td>85,981,153</td><td>80,602,593</td><td>93.72</td></tr><tr><td>Fayoumi</td><td>Treated</td><td>CH&amp;NDV</td><td>84,433,550</td><td>78,218,756</td><td>93.76</td></tr><tr><td rowspan="8">Breast muscle</td><td>Leghorn</td><td>Non-treated</td><td>Acute heat</td><td>202,976,814</td><td>171,248,763</td><td>84.65</td></tr><tr><td>Leghorn</td><td>Treated</td><td>Acute heat</td><td>197,916,434</td><td>168,158,916</td><td>85.16</td></tr><tr><td>Fayoumi</td><td>Non-treated</td><td>Acute heat</td><td>199,002,995</td><td>172,003,441</td><td>86.63</td></tr><tr><td>Fayoumi</td><td>Treated</td><td>Acute heat</td><td>167,292,872</td><td>145,714,038</td><td>86.93</td></tr><tr><td>Leghorn</td><td>Non-treated</td><td>CH&amp;NDV</td><td>97,856,074</td><td>85,109,242</td><td>99.91</td></tr><tr><td>Leghorn</td><td>Treated</td><td>CH&amp;NDV</td><td>101,741,687</td><td>87,181,252</td><td>99.91</td></tr><tr><td>Fayoumi</td><td>Non-treated</td><td>CH&amp;NDV</td><td>100,416,625</td><td>86,497,655</td><td>99.91</td></tr><tr><td>Fayoumi</td><td>Treated</td><td>CH&amp;NDV</td><td>102,815,215</td><td>88,592,510</td><td>99.92</td></tr></tbody></table>

by the Illumina HiSeq 4000, which included four treatment groups: treated and non-treated Leghorn (LT and LC); treated and non-treated Fayoumi (FT and FC), at the acute heat stress (AH) and the chronic heat stress &amp; NDV inoculation (CH&amp;NDV), respectively. Of the 24,357 annotated chicken genes in the chicken Galgal 6.0 database, more than 80% of the annotated genes were identified in both lines. The detailed mapping statistics are listed in Table 3. The two tissues showed distinct transcriptome profiles. Many genes (2,242 at the AH stage and 3,052 at the CH&amp;NDV stage) were specifically identified in the hypothalamus. For the breast muscle, 627 genes at the AH stage and 384 genes at the CH&amp;NDV stage were specifically expressed compared to the hypothalamus, many of which are related to muscle biogenesis (online suppl. Material 6).

*Differential Gene Expression in Different Comparisons*

Gene expression profiles were compared between each pair of the four groups with two treatments in the two tissues, respectively. The comparison groups were the same as previously reported, which were FC versus LC (FCLC), FT versus LT (FTLT), LT versus LC (LTLC) and FT versus FC (FTFC) at AH and CH&amp;NDV for both hypothalamus and breast muscle [Wang Y et al., 2020]. The number of DEGs is shown in Figure 14 for the between-line comparisons and Figure 15 for the within-line comparisons of the two tissues. No genes were identified by the DEG analysis which demonstrated an interaction effect in each

tissue. Detailed gene information and fold changes in each comparison are listed in online supplementary Material 7.

For DEGs between genetic lines, upregulated genes with higher expression in Fayoumi birds were Fayoumi-favored genes, while down-regulated genes with higher expression in Leghorn birds were Leghorn-favored genes. In the hypothalamus, more Leghorn-favored genes were identified than Fayoumi-favored genes at the AH stage, however, the pattern was switched by having more Fayoumi-favored DEGs than Leghorn's at the CH stage (Fig. 14). Meanwhile, more DEGs were identified in the non-treated comparisons than in the treated comparisons. The numbers of DEGs identified between the treated Leghorn and Fayoumi line dramatically decreased at the CH&amp;NDV stage with 320 Fayoumi-favored genes and 224 Leghorn-favored genes. In the breast muscle, there were always more Fayoumi-favored genes than Leghorn-favored genes across all comparisons (Fig. 14).

There were fewer DEGs within genetic lines than between lines, the DEG numbers were significantly decreased especially in the hypothalamus compared to the between-line comparisons (Fig. 15). More DEGs were upregulated at the AH stage for both lines, as compared to the CH&amp;NDV stage in which more DEGs were downregulated. There were no DEGs shared by the four within-line comparisons in the hypothalamus (Fig. 16). Ten DEGs, including 3 heat shock protein family members (heat shock protein family B member 9 [HSPB9],
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Numbers of differentially expressed genes identified between genetic lines. A false discovery rate <0.05 was used to classify genes as differentially expressed. AH, acute heat stress; CH, chronic heat stress and NDV infection at 2 dpi; FCLC, Fayoumi non-treated versus Leghorn non-treated; FTLT, Fayoumi treated versus Leghorn treated.</figcaption></figure>

heat shock 70kDa protein 2 [HSPA2], and heat shock protein 90 beta family member 1 [HSP90B1]), were shared by the Leghorn and Fayoumi birds at the AH stage (Table 4). Four DEGs were identified in both lines with the CH&NDV treatment, three of which were hemoglobin genes (Table 4). The list of DEGs specifically identified in each comparison provided more information about the line-specific gene regulatory response to heat stress under NDV infection. Six heat shock related genes were specifically identified in the LTLCAH contrast, 3 immune related genes (T cell leukemia homeobox 3 [TLX3], interferon alpha inducible protein 6 [IFI6], and hemoglobin subunit epsilon 1 [HBE1]) were specifically identified in the LTLCCH group, 2 heat shock related genes (heat shock protein 30C-like [HSP30C] and heat shock protein 90 alpha family class B member 1 [HSP90AB1]), and 3 immune-related genes (Myxovirus resistance 1, interferon-inducible protein [MX1], radical S-adenosyl methionine domain con-

taining 2 [RSAD2], and interferon induced protein with tetratricopeptide repeats 5 [IFIT5]) were only identified in the FTFCAH contrast, and 7 metabolism-associated genes were only identified in the FTFCCH contrast (Table 5).

The DEG numbers were higher in the breast muscle than hypothalamus, particularly during the CH&NDV stage, in which there were more upregulated DEGs than downregulated DEGs for all comparisons (Fig. 15). The heat shock protein family H member 1 (HSPH1) gene was identified in all four groups (Fig. 17). Nine DEGs were identified in both lines at the AH stage and 200 DEGs were shared by the two lines at the CH&NDV stage. Of the line-specific DEGs, 53 DEGs for the LTLCAH, 33 DEGs for the FTFCAH, 316 for the LTLCCH&NDV, and 677 DEGs for the FTFCCH&NDV, the majority of DEGs were metabolism-related genes except a few immune-related genes such as interferon regulatory factor 2 (IRF2), interferon alpha inducible protein 6 (IFI6) and interferon regulatory
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Numbers of differentially expressed genes identified within genetic lines. A false discovery rate <0.05 was used to classify genes as differentially expressed. AH, acute heat stress; CH, chronic heat stress and NDV infection at 2 dpi; LTLC, Leghorn treated versus non-treated; FTFC, Fayoumi treated versus non-treated.</figcaption></figure>

factor 2 binding protein 2 (IRFBP2). Detailed gene information and fold changes for breast muscle line-specific DEGs are listed in online supplementary Material 8.

## Functional Categories of Differentially Expressed Genes

Gene ontology (GO) was used to evaluate the function of DEGs from different comparisons. Because characterizing the line-specific responses to heat stress under NDV infection was one of our major objectives, the gene function analysis focused on the within-line comparisons in the two lines and two tissues. All DEGs in the within-line comparisons were performed by functional enrichment analysis through the DAVID program (DAVID Bioinformatics Resources 6.8). The biological process and KEGG pathways were presented as functional clusters. The significant enriched GO terms and KEGG pathways were presented if *p* < 0.05 and FDR <20%.

## Hypothalamus

In the Leghorn line, even with fewer DEGs for the within-line comparisons identified in the hypothalamus, 7 GO terms including the protein folding, and one KEGG pathway (protein processing in endoplasmic reticulum) were enriched by the upregulated DEGs in the LTLCAH comparison. Most of them are metabolic-associated except 1 GO term, peptide antigen assembly with MHC class I protein complex, which was immune related (Fig. 18a). More GO terms (20) were significantly enriched by the downregulated DEGs in the Leghorn birds with the AH treatment, which included iron transport and cholesterol and cellular glucose homeostasis functions (Fig. 18b). No GO terms were enriched by the uregulated DEGs in the Leghorn line with the CH&NDV treatment. Three metabolic functional GO terms were enriched by the downregulated DEGs in the LTLCCH comparison (Fig. 18c).
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Venn diagram of differentially expressed genes within genetic lines in the hypothalamus. LTLCAH, Leghorn treated versus non-treated with acute heat stress; FTFCAH, Fayoumi treated versus non-treated with acute heat stress; LTLCCH&NDV, Leghorn treated versus non-treated with chronic heat stress and 2 dpi NDV infection; FTFCCH&NDV, Fayoumi treated versus non-treated with chronic heat stress and 2 dpi NDV infection.</figcaption></figure>

In the Fayoumi line, 2 metabolic functions were enriched by the upregulated DEGs with the AH treatment (Fig. 18d), same as the Leghorn line having the protein folding function enriched. Seven GO terms were enriched by the downregulated DEGs in the FTFCAH comparison, which included defense response to virus (Fig. 18e). With the CH&NDV treatment, 11 GO terms and 1 KEGG pathway (ABC transporters) were enriched by the upregulated DEGs, and 2 metabolic functional terms were enriched by the downregulated DEGs (Fig. 18f, g), and all of them were associated with metabolic and catabolic processes.

## Breast Muscle

At the AH stage, upregulated DEGs in Leghorns enriched 5 GO terms including skeletal muscle cell differentiation and negative regulation of apoptotic signaling pathway (Fig. 19a). Downregulated Leghorn DEGs enriched 1 GO term, positive regulation of vasodilation, and 1 KEGG pathway, alanine, aspartate, and glutamate metabolism (Fig. 19b). Four GO terms were enriched by the upregulated DEGs in Fayoumi birds, which included both metabolic-related functions, such as responding to cold, and immune-related functions, such

<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Venn diagram of differentially expressed genes within genetic lines in the breast muscle. LTLCAH, Leghorn treated versus non-treated with acute heat stress; FTFCAH, Fayoumi treated versus non-treated with acute heat stress; LTLCCH&NDV, Leghorn treated versus non-treated with chronic heat stress and 2 dpi NDV infection; FTFCCH&NDV, Fayoumi treated versus non-treated with chronic heat stress and 2 dpi NDV infection.</figcaption></figure>

as positive regulation of T-cell activation (Fig. 19c). No GO terms were enriched by the downregulated Fayoumi DEGs.

At the CH&NDV stage, many GO terms and pathways were enriched by DEGs identified in the within-line comparisons due to the larger number of DEGs. Seventeen GO terms and 2 KEGG pathways were enriched by the upregulated Leghorn DEG genes (Fig. 19d). The innate immune response is one of the enriched functions. Downregulated Leghorn DEGs enriched 7 GO terms and 8 KEGG pathways, all of which are metabolic-associated functions, except the biosynthesis of antibiotics enriched by the downregulated Leghorn DEGs (Fig. 19e). In the Fayoumi line, both up-and downregulated DEGs are involved in metabolism functions such as skeletal muscle cell differentiation by upregulated DEGs and fatty acid beta-oxidation using acyl-CoA by the downregulated DEGs (Fig. 19f, g). Two immune related KEGG pathways were enriched by Fayoumi DEGs at this stage. Adipocytokine signaling pathway was enriched by upregulated Fayoumi DEGs and biosynthesis of antibiotics was enriched by downregulated Fayoumi DEGs (Fig. 19g).
<table><caption>**Table 4.** DEGs identified in multiple contrasts in the hypothalamus</caption><thead><tr><th>DEG</th><th>Gene description</th><th>LTLC AH</th><th>FTFC AH</th><th>LTLC CH&amp;NDV</th><th>FTFC CH&amp;NDV</th></tr></thead><tbody><tr><td>HSPB9</td><td>Heat shock protein family B (small) member 9</td><td><strong>2.58<sup>a</sup></strong></td><td><strong>3.06</strong></td><td></td><td></td></tr><tr><td>HSPA2</td><td>Heat shock 70kDa protein 2</td><td><strong>1.23</strong></td><td><strong>1.07</strong></td><td></td><td></td></tr><tr><td>CRELD2</td><td>Cysteine rich with EGF like domains 2</td><td><strong>1.07</strong></td><td><strong>1.04</strong></td><td></td><td></td></tr><tr><td>HSP90B1</td><td>Heat shock protein 90 beta family member 1</td><td><strong>0.96</strong></td><td><strong>0.80</strong></td><td></td><td></td></tr><tr><td>GPX1</td><td>Glutathione peroxidase</td><td><strong>0.77</strong></td><td><strong>0.94</strong></td><td></td><td></td></tr><tr><td>CALR</td><td>Calreticulin</td><td><strong>0.64</strong></td><td><strong>0.90</strong></td><td></td><td></td></tr><tr><td>SDF2L1</td><td>Stromal cell derived factor 2 like 1</td><td><strong>0.58</strong></td><td><strong>0.59</strong></td><td></td><td></td></tr><tr><td>MANF</td><td>Mesencephalic astrocyte derived neurotrophic factor</td><td><strong>0.57</strong></td><td><strong>0.89</strong></td><td></td><td></td></tr><tr><td>HYOU1</td><td>Hypoxia up regulated 1</td><td><strong>0.37</strong></td><td><strong>0.43</strong></td><td></td><td></td></tr><tr><td>PDIA3</td><td>Protein disulfide isomerase family A member 3</td><td><strong>0.28</strong></td><td><strong>0.38</strong></td><td></td><td></td></tr><tr><td>HBBA</td><td>Hemoglobin beta, subunit A</td><td></td><td></td><td>-1.18<sup>b</sup></td><td>-1.54</td></tr><tr><td>HBA1</td><td>Hemoglobin subunit alpha 1</td><td></td><td></td><td>-1.26</td><td>-1.47</td></tr><tr><td>HBAD</td><td>Hemoglobin alpha, subunit D</td><td></td><td></td><td>-1.27</td><td>-1.59</td></tr><tr><td>EPB42</td><td>Erythrocyte membrane protein band 4.2</td><td></td><td></td><td>-2.73</td><td>-1.77</td></tr></tbody></table>

<table><caption>**Table 5.** Contrast specific DEGs identified in hypothalamus</caption><thead><tr><th>DEG</th><th>Gene description</th><th>LTLC AH</th><th>FTFC AH</th><th>LTLC CH&amp;NDV</th><th>FTFC CH&amp;NDV</th></tr></thead><tbody><tr><td>CLIC5</td><td>Chloride intracellular channel 5</td><td><strong>1.03<sup>a</sup></strong></td><td></td><td></td><td></td></tr><tr><td>HSPB8</td><td>Heat shock protein family B member 8</td><td><strong>0.89</strong></td><td></td><td></td><td></td></tr><tr><td>BAG3</td><td>BCL2 associated athanogene 3</td><td><strong>0.67</strong></td><td></td><td></td><td></td></tr><tr><td>HSPH1</td><td>Heat shock protein family H member 1</td><td><strong>0.64</strong></td><td></td><td></td><td></td></tr><tr><td>HSPA5</td><td>Heat shock 70kDa protein 5</td><td><strong>0.43</strong></td><td></td><td></td><td></td></tr><tr><td>DNAJB6</td><td>DnaJ heat shock protein family (Hsp40) member B6</td><td><strong>0.24</strong></td><td></td><td></td><td></td></tr><tr><td>DDX3X</td><td>DEAD-box helicase 3, X-linked</td><td><strong>0.22</strong></td><td></td><td></td><td></td></tr><tr><td>DNAJC3</td><td>DnaJ heat shock protein family member C3</td><td><strong>0.21</strong></td><td></td><td></td><td></td></tr><tr><td>SLC38A2</td><td>Solute carrier family 38 member 2</td><td>-0.30<sup>b</sup></td><td></td><td></td><td></td></tr><tr><td>NR1D1</td><td>Nuclear receptor subfamily 1 group D member 1</td><td>-0.61</td><td></td><td></td><td></td></tr><tr><td>HSP30C</td><td>Heat shock protein 30C-like</td><td></td><td><strong>2.23</strong></td><td></td><td></td></tr><tr><td>HSP90AB1</td><td>Heat shock protein 90 alpha family class B member 1</td><td></td><td>-0.47</td><td></td><td></td></tr><tr><td>MX1</td><td>Myxovirus resistance 1, interferon-inducible protein</td><td></td><td>-0.85</td><td></td><td></td></tr><tr><td>RSAD2</td><td>Radical S-adenosyl methionine domain containing 2</td><td></td><td>-1.50</td><td></td><td></td></tr><tr><td>IFIT5</td><td>Interferon induced protein with tetratricopeptide repeats 5</td><td></td><td>-1.60</td><td></td><td></td></tr><tr><td>TLX3</td><td>T cell leukemia homeobox 3</td><td></td><td></td><td><strong>1.46</strong></td><td></td></tr><tr><td>IFI6</td><td>Interferon alpha inducible protein 6</td><td></td><td></td><td><strong>1.17</strong></td><td></td></tr><tr><td>HBE1</td><td>Hemoglobin subunit epsilon 1</td><td></td><td></td><td>-2.08</td><td></td></tr><tr><td>SLC2A12</td><td>Solute carrier family 2 member 12</td><td></td><td></td><td></td><td><strong>1.15</strong></td></tr><tr><td>MRC1</td><td>macrophage mannose receptor 1-like</td><td></td><td></td><td></td><td><strong>1.05</strong></td></tr><tr><td>SLC13A5</td><td>Solute carrier family 13 member 5</td><td></td><td></td><td></td><td><strong>0.91</strong></td></tr><tr><td>TGM2</td><td>Transglutaminase 2</td><td></td><td></td><td></td><td><strong>0.71</strong></td></tr><tr><td>CRCSB</td><td>Cytokine receptor common subunit beta-like</td><td></td><td></td><td></td><td><strong>0.64</strong></td></tr><tr><td>HS3ST1</td><td>Heparan sulfate-glucosamine 3-sulfotransferase 1</td><td></td><td></td><td></td><td><strong>0.43</strong></td></tr><tr><td>NKAIN4</td><td>Sodium/potassium transporting ATPase interacting 4</td><td></td><td></td><td></td><td>-0.22</td></tr></tbody></table>

<footnote>a Log2 fold change; bold values indicate upregulation. b Negative values (–) indicate downregulation.</footnote>

<footnote>a Log2 fold change; bold values indicate upregulation. b Negative values (–) indicate downregulation.</footnote>
<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Gene ontology (GO) biological processes and KEGG pathway overrepresentation (*p* < 0.05 and FDR <20%) for within-line comparisons in the hypothalamus. **a** GO terms and KEGG pathways significantly enriched by upregulated genes in the LTLCAH comparison. **b** GO terms and KEGG pathways significantly enriched by downregulated genes in the LTLCAH comparison. **c** GO terms and KEGG pathways significantly enriched by downregulated genes in</figcaption></figure>

the LTLCCH comparison. **d** GO terms and KEGG pathways significantly enriched by upregulated genes in the FTFCAH comparison. **e** GO terms and KEGG pathways significantly enriched by downregulated genes in the FTFCAH comparison. **f** GO terms and KEGG pathways significantly enriched by upregulated genes in the FT-FCCCH comparison. **g** GO terms and KEGG pathways significantly enriched by downregulated genes in the FT-FCCCH comparison.
<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Gene ontology (GO) biological processes and KEGG pathway overrepresentation (*p* < 0.05 and FDR <20%) for within-line comparisons in the breast muscle. **a** GO terms and KEGG pathways significantly enriched by upregulated genes in the LTLCAH comparison. **b** GO terms and KEGG pathways significantly enriched by downregulated genes in the LTLCAH comparison. **c** GO terms and KEGG pathways significantly enriched by upregulated genes in the FTFCAH comparison. **d** GO terms and KEGG pathways significantly enriched by downregulated genes in the LTLCAH comparison.</figcaption></figure>

ways significantly enriched by upregulated genes in the LTLCCH comparison. **e** GO terms and KEGG pathways significantly enriched by downregulated genes in the LTLCCH comparison. **f** GO terms and KEGG pathways significantly enriched by upregulated genes in the FTFCCH comparison. **g** GO terms and KEGG pathways significantly enriched by downregulated genes in the FT-FCCCH comparison.
<figure><img src="image_22.png" /><figcaption>**Fig. 20.** Comparative analysis of significantly enriched canonical pathways through Ingenuity Pathway Analysis among differentially expressed genes by genetic line and treatment in the breast muscle (*p* < 0.05 and *z* >1.64], where orange (positive *z*-score) refers to predicted activation and blue (negative *z*-score) to predicted inhibition.</figcaption></figure>

### Gene Network and Signaling Pathway Analysis

In addition to the GO analysis, the Ingenuity Pathway Analysis was used to identify gene networks and signaling pathways by using significant DEGs. Because of the limited DEG numbers in the hypothalamus, there was only one canonical pathway, unfolded protein response, identified in the Leghorn line with the AH treatment. In the breast muscle, no canonical pathways were identified in the Leghorn and Fayoumi lines at the AH stage. Twenty-five canonical pathways were identified in the two genet-

ic lines with the CH&NDV treatment (Fig. 20). Thirteen out of 25 had similar activation and inhibition patterns. Five pathways were only activated in Leghorn birds, and 3 pathways were only identified in Fayoumi birds. Four pathways, phospholipase C signaling, calcium signaling, androgen signaling and synaptic long-term potential, were activated in Leghorns and inhibited in the Fayoumi line with the CH&NDV treatment.

### Weighted Gene Co-Expression Network Analysis

To overcome the limitation of DEG numbers specifically identified in the hypothalamus, a correlation-based WGCNA analysis was conducted. Seven gene modules were identified in total by using all the hypothalamus and breast muscle data (Fig. 21). Four out of 7 gene modules, black, turquoise, blue, and yellow modules, showed significant correlations with tissue, genetic lines, the CH&NDV treatment, and phenotypic traits. Four trait-correlated modules were further analyzed and the results are described in Table 6. The black module was the only one positively correlated with tissue, which means genes from the black module had higher expression levels in the hypothalamus. The turquoise module was only negatively correlated with tissue. Genes in the turquoise module had higher expression levels in the breast muscle. The blue module was positively correlated with line and Na⁺ and negatively correlated with pH, HCO₃, TCO₂, BE, PO₂, sO₂, and glucose levels in chicken blood. Genes in the blue module are highly expressed in Leghorn birds and correlated with higher sodium and lower blood gas levels. The yellow module was negatively correlated with the CH&NDV treatment and the blood PCO₂ level.

Furthermore, genes driving the biology within these trait-correlated modules were identified to understand the potential biological processes these two genetic lines underwent during treatments. Driver genes were not discovered in all correlated gene modules. Driver genes were only identified from tissue-correlated black and turquoise modules, the CH&NDV treatment-correlated yellow module, and the line, pH, PO₂, sO₂, and Glu correlated blue module. The top-5 driver genes with the highest absolute gene significance (GS) and module membership (MM) are listed in Table 7. All the detailed driver gene information is presented in online supplementary Material 9 and 10.

For example, driver genes, glutamate ionotropic receptor NMDA type subunit1 (GRIN1), tropomyosin 1 (TMP1), solute carrier family 17 member 6 (SLC17A6), and calcium voltage-gated channel auxiliary subunit