Genomic status of yellow-breasted bunting following recent rapid population decline.

Global biodiversity is facing serious threats. However, knowledge of the genomic consequences of recent rapid population declines of wild organisms is limited. Do populations experiencing recent rapid population decline have the same genomic status as wild populations that experience long-term declines? Yellow-breasted Bunting (Emberiza aureola) is a critically endangered species that has been experiencing a recent rapid population decline. To answer the question, we assembled and annotated the whole genome of Yellow-breasted Bunting. Furthermore, we found high genetic diversity, low linkage disequilibrium, and low proportion of long runs of homozygosity in Yellow-breasted Bunting, suggesting that the populations following recent rapid declines have different genomic statuses from the population that experienced long-term population decline.

Introduction

Wild species are suffering population declines and some species have gone extinct (IUCN, 2021). 2.9 billion breeding adult birds have been lost since 1970 in North America, with particularly strong losses among grassland birds (53%) (Rosenberg et al., 2019). It is, therefore, urgent to understand the genomic consequences of the rapid population decline in wild organisms. Genetic status, including genetic diversity and inbreeding level, provides important information for planning conservation strategy (DeWoody et al., 2021). Endangered species with different genetic characteristics may have different survival fates. However, few at-risk species have had their genomic status evaluated. What is worse, it is still unclear what the genomic status of recent rapid population declines of wild species is.

With the development of genomic sequencing, the availability of population genomic data allows the exploration of genetic diversity, inbreeding level, mutation load, and historical effective population size fluctuations. Previous studies have found endangered species have low genetic diversity (e.g., Crossoptilon mantchuricum), high inbreeding levels (e.g., Pavo muticus), and damaging mutation loads (e.g., Nipponia nippon), but these species have experienced long-term population decline (Dong et al., 2021; Feng et al., 2019; Wang et al., 2021). Limited studies have explored the genomic status of species that have been subjected to recent rapid population declines.

Yellow-breasted Bunting (Emberiza aureola; Emberizidae; Passeriformes) is an important endangered species because it experienced a recent rapidly decline. It was one of the most abundant songbirds in northern Eurasia as recently as 1980s (Kamp et al., 2015). They breed from Finland in the west to Kamchatka in the east, and winter in South and Southeast Asia (International, 2021). Yellow-breasted Bunting was even considered a serious agricultural pest in the 20th century, as it often formed large flocks thousands-strong and fed on wheat and rice during the non-breeding season in China (Shaw, 1936; Zheng, 1956). However, this species is now considered to be critically endangered (International, 2021), with a global population decline of 84.3–94.7% from 1980 to 2013 (Kamp et al., 2015), which was mainly attributed to the excessive trapping at migration and wintering sites (IUCN, 2021; Kamp et al., 2015). As the speed and magnitude of the population decline is unprecedented, there was concern about its extinction. It is unclear whether the bunting retains the population’s genetic capacity to recover.

To explore the genomic status of the Yellow-breasted bunting following its recent rapid population decline, and the species' genetic capacity to recover, we de novo assembled a whole genome of Yellow-breasted Bunting and sequenced whole genomes of ten individuals to estimate the genetic diversity, inbreeding level, potential mutational load, and demographic history. The current study will provide genomic information for the conservation of Yellow-breasted Bunting.

Results

Genome assembly and annotation

To estimate the whole genomic status and the amount of deleterious mutations of Yellow-breasted Bunting, we de novo assembled the genome and annotated its protein-coding regions. The length of the assembled Yellow-breasted Bunting genome was 1,149,406,387 bp (1.15 Gb) with the longest scaffold 37,240,703 bp long. The coverage of the Pacbio subreads was almost 162 ×. The GC content in the Yellow-breasted Bunting genome was 43.16%. The scaffold N50 of the assembly was 8,222,108 bp (8.22 Mbp). The assembly covered 95.7% of BUSCOs and 91.94% of core gene, suggesting that the completeness of the assembly was good. The mapping rate and coverage of the Illumina reads were 98.84% and 99.87%, respectively, suggesting that the assembly has good consistency, integrity, and accuracy. Based on the alignment information, 153 scaffolds were identified as sex chromosome-linked scaffolds. The total length of these scaffolds was 106,560,505 bp, approximately 9.27% of the whole genome. The new assembly has been deposited in the Genome Warehouse in National Genomics Data Center (https://ngdc.cncb.ac.cn/gwh), Beijing Institute of Genomics, Chinese Academy of Sciences/China National Center for Bioinformation, NGDC: GWHBJBD00000000.

In the genome of Yellow-breasted Bunting, 18.86% genome was identified as repeat sequences (Table S1). Based on the combined evidence, the Yellow-breasted Bunting genome was predicted to have 18,348 genes, among which 84.30% of genes assigned functions (Figures S1 and S2). The average transcript length and CDS lengths were 23,435.04 bp and 1,531.24 bp, respectively. There are 0.02% of the genome was predicted to be non-coding RNA (Table S2).

Yellow-breasted bunting has high genetic diversity

To perform population genetic analysis, we re-sequenced 10 Yellow-breasted Buntings (Table S3). The mean coverage and mapping rates of the 10 Yellow-breasted Buntings were 20.59 (standard deviation (SD) = 1.31) and 99.37% (SD = 0.10%), respectively (Table S3). After filtering, 18,938,565 high-quality SNPs were used in the population genetic analysis. There were no paired samples with relationship closer than 3rd-degree (Figure S3), so the 10 samples were all used in the population genetic analysis.

To understand the genomic status of Yellow-breasted Bunting among Passeriformes, we compared its genomic heterozygosity with 17 Passeriformes birds. The 17 Passeriformes birds included Florida scrub jay (Aphelocoma coerulescens), Noisy scrubbird (Atrichornis clamosus), Chestnut-collared longspur (Calcarius ornatus), Seychelles magpie robin (Copsychus sechellarum), Ribbon-tailed drongo (Dicrurus megarhynchus), Painted honeyeater (Grantiella picta), Akiapola’au (Hemignathus wilsoni), Philippine fairy-bluebird (Irena cyanogastra), Kirtland’s warbler (Setophaga kirtlandii), Loggerhead shrike (Lanius ludovicianu), Bali myna (Leucopsar rothschildi), Yellowhead (Mohoua ochrocephala), Velvet myiagra (Myiagra hebetior), Inaccessible island finch (Nesospiza acunhae), Stitchbird (Notiomystis cincta), White-necked rockfowl (Picathartes gymnocephalus), and Golden-crowned babbler (Sterrhoptilus dennistouni). The result tells us the Yellow-breasted Bunting has high heterozygosity. The mean heterozygosity of the 10 samples is 0.46 × 10−2 (0.44 × 10−2–0.47 × 10−2) (Table S4). The vulnerable chestnut-collared longspur has higher heterozygosity than Yellow-breasted Bunting (Figure 1, Table S4). Excepting the chestnut-collared longspur, all other Passeriformes assessed had lower heterozygosity than the Yellow-breasted Bunting, no matter their IUCN threatened status (Figure 1, Table S4).

Figure 1

Heterozygosity of Yellow-breasted Bunting (indicated by yellow bird icon) compared to that of 17 other Passeriformes bird species

Heterozygosity equals the number of heterozygous SNPs divided by the length of the genome without gaps. The red, brown, sepia, and black colors correspond to critically endangered, endangered, vulnerable, and near-threatened status, respectively.

Yellow-breasted bunting is not severely inbred

To understand the inbreeding pattern of Yellow-breasted Bunting, we constructed the linkage disequilibrium (LD) decay curve and calculated the amount of runs of homozygosity (ROH) in the species. The LD decay curve becomes smooth at a pairwise distance of less than 100 kb (Figure 2A). The maximum correlation coefficient between any two pairs of SNPs in a 100 bp windows is 0.18. The average total length of long and short ROH in the 10 samples is 1143.68 kb (SD = 2458.15) and 6153.12 kb (SD = 2165.16), respectively. The average percentage of long and short ROH on autosomes is 0.11% (SD = 0.25%) and 0.59% (SD = 0.22%), respectively (Figure 2B). The total length of short ROH is significantly larger than that of long ROH (KS test, p < 0.01). Both the LD decay curve and the percentage of ROH on autosomes suggest inbreeding is not serious in Yellow-breasted Bunting (Figure 2).

Figure 2

Linkage disequilibrium (LD) decay and runs of homozygosity (ROH) in Yellow-breasted Bunting

(A) Linkage disequilibrium decay of Yellow-breasted Bunting. The r2 is a correlation coefficient between any pair of SNPs.

(B) Box plot of runs of homozygosity (ROH) in Yellow-breasted Bunting. The percentage is the total length of ROH divided by the length of autosomes x 100%. Long ROH are defined as those longer than 1.5 Mb and short ROH as shorter than 1 Mb. A Kolmogorov-Smirnov test suggests the total length of short ROH is significantly larger than that of long ROH. ∗∗ means the difference is significant.

Potential mutational load

To recover the genetic load in Yellow-breasted Bunting, we identified the synonymous, missense, and lose of function (LOF) mutations in each sample. There were 130,622 synonymous, 61,977 missenses, 813 LOF mutations in Yellow-breasted Bunting. The frequency of homozygous sites is less than that of heterozygous sites for all three kinds of mutations (Figure 3A). Among 1000 random SNPs in coding sites, there was no significant difference in the relative frequency of homozygous sites between synonymous and missense mutations (Figure 3A). However, the relative frequency of homozygous sites of LOF mutations was significantly less than that of synonymous mutations (Figure 3A), suggesting that most LOF mutations are potential homozygous lethal mutations. There were 743 genes harboring LOF mutations, and they were significantly enriched in six KEGG pathways (Figure 3B). The six pathways are oocyte meiosis (gga04114), metabolic pathways (gga01100), NOD-like receptor signaling pathway (gga04621), apoptosis (gga04210), ECM-receptor interaction (gga04512), and C-type lectin receptor signaling pathways (gga04625).

Figure 3

Mutation load and enriched KEGG pathways of the genes that harbor loss-of-function (LOF) mutations in Yellow-breasted Bunting

(A) Average difference between homozygous site frequency and heterozygous site frequency. The y-axis is the relative frequency of homozygous sites and negative values mean the frequency of heterozygous sites is larger than that of homozygous sites. S, synonymous mutations; M, missense mutations; L, loss of function mutations. Kolmogorov-Smirnov comparison was used to calculate the p value.

(B) The significant enriched KEGG pathways. Fisher's exact test was used to calculate the p value.

Yellow-breasted bunting began to decline 100–200 generations ago

To reconstruct the demographic history of the Yellow-breasted bunting, we employed three methods: PSMC, stairway plot, and fastsimcoal. The complementary methods provide a detailed demographic history of Yellow-breasted Bunting, suggesting that Yellow-breasted Bunting experienced a bottleneck, and that the most recent population decline started 100–200 generations ago (Figures 4 and 5C). PSMC and the stairway plot both suggest Yellow-breasted Bunting began to expand about 2 × 106 years ago (Figure 4). After a period of population stability, a decline began about 4 × 105 years ago (Figures 4A and 4B). Recovery likely began about 2 × 105 years ago (Figure 4A). The best fit model inferred by fastsimcoal is consistent with the results of PSMC and stairway plot (Figure 5C). The best fit model supports Yellow-breasted Bunting having experienced a bottleneck and the recently decline of the population started at about 147 (95% confidence interval (95% CI) = 146–149) generations ago. The average current population size of Yellow-breasted Bunting is 5119488 (95% CI = 5117928–5121048). Before the recent population decline, the population size of Yellow-breasted Bunting was 7320650 (95% CI = 7319004–7322296). Nearly 30.07% of the population disappeared. The Akaike Information Criterion (AIC) value of the best model is 49761.58, and its average likelihood from 100 replications was a significant difference from other models (p < 0.01).

Figure 4

The fluctuations of historical effective population size (Ne)

(A) Historical Ne inferred by stairway plot. The red line represents the median inferred Ne. The grey line represents the 75% confidence interval.

(B) Historical Ne inferred by pairwise sequentially Markovian coalescent (PSMC). The numbers before the colored lines are the numbers in the sample IDs in Table S3. The generation time was assumed to be two years. The mutation rate was set as 3.3 × 10−9 per site per generation.

Figure 5

The demographic history models estimated by fastsimcoal

(A) Bottleneck and stable model (BS).

(B) Bottleneck and decline model (BD).

(C) Bottleneck, ancestry stable, and early decline model (BSED).

(D) Bottleneck, ancestry stable, and recently decline model (BSRD). NPOA, NPOP, NPOS, NPOG, and NPOB represent the effective population size (Ne) of the different historical stages. TG, TB, TA, and TC represent the time that Ne began to change; Current, 2019 A. D.; GRO, population decline rate; g, generations ago; AIC, Akaike Information Criterion.

Discussion

The current project assembled the whole genome of Yellow-breasted Bunting and combined it with population genomic data to infer the genomic status of the recent rapid population decline of Yellow-breasted Bunting, and to determine whether the species could recover its abundant population size. The results suggest Yellow-breasted Bunting has rich genetic diversity, and low levels of inbreeding, but potentially homozygous lethal mutations. Demographic history shows Yellow-breasted Bunting has experienced a bottleneck and the last population decline began about 147 generations ago.

Genomic status following recent rapid population decline differs from long-term decline

Wild populations that experienced long-term declines may have low genetic diversity, high levels of inbreeding, and genetic load. In small populations, genetic drift would eliminate genetic diversity (Ellstrand and Elam, 1993; Feng et al., 2019). Individuals in small populations have much more opportunity to be inbred compared to that in large populations. In turn, inbreeding leads small populations to accumulate deleterious mutations (Charlesworth and Willis, 2009). Empirical studies have found that wild populations that experienced long-term populations have low genetic diversity, high frequency of inbreeding, and serious mutation loads. Such species include Brown Eared Pheasant (Crossoptilon mantchuricum) (Wang et al., 2021) and Crested Ibis (Nipponia nippon) (Feng et al., 2019). Long-term declines make populations have small population size, which may facilitate genetic drift toward low genetic diversity, high levels of inbreeding, and genetic load.

Compared with populations that experienced long-term declines, those experiencing recent rapid declines may have high genetic diversity, low levels of inbreeding, and low mutation load. The populations declining recently and rapidly may have relatively large population sizes, in which case genetic drift has less power to erase the genetic variations. Over short time frames, the population may accumulate most genetic characteristics from its parental population facing genetic drift.

The Yellow-breasted Bunting population began to decline about 100–200 generations ago (Figure 5). Even though genetic drift could act on the Yellow-breasted Bunting, the species has high genetic diversity (Figure 1). During most of the recent decline, the Yellow-breasted Bunting may have had a large population size, which makes genetic drift have less power to eliminate the genetic diversity. The short bottleneck, nearly 147 generations, may be the other reason that explains Yellow-breasted Bunting has high genetic diversity. The declining population may need a long time to lose genetic diversity. The Yellow-breasted Bunting has low levels of inbreeding (Figure 2), which may be indirect evidence to support a large recent population size because individuals have less chance of being inbred with each other in larger populations. The low level of inbreeding will help the population have only a slight mutation load. The result suggests Yellow-breasted Bunting has potential homozygous lethal mutations (Figure 3), which may be from a historical bottleneck (Figure 4).

The yellow-breasted bunting has the capability to recover

Genetic diversity is an important metric in conservation and biodiversity, because it is intimately tied to fitness and adaptive potential (DeWoody et al., 2021). In taxonomically diverse organisms, such as fruit flies (Drosophila melanogaster) (Frankham, 1995), field mice (Peromyscus leucopus) (Lacy et al., 2013), and monkeyflower (Mimulus guttatus) (Willis, 1993), genetic diversity is positively related to fitness (Chapman et al., 2009). In threatened desert tortoises (Gopherus agassizii), individual heterozygosity predicts translocation success (Scott et al., 2020). Besides fitness, adaptive potential is correlated with genetic diversity. For example, the genetic diversity in the coding region and conserved non-coding region is low in Brown Eared Pheasant, a species with low genome-wide diversity (Wang et al., 2021). Moreover, neutral diversity in the genome provides a genetic template for the future evolution of adaptive diversity (Harrisson et al., 2014). The rich genetic diversity in the Yellow-breasted Bunting population (Figure 1) suggests it has high fitness and adaptive potential relative to other endangered species, such as Painted Honeyeater (Grantiella picta) and Velvet Myiagra (Myiagra hebetior) (Figure 1), which could contribute to its recovery. For example, the genome-wide heterozygosity of crested ibis is 0.04 × 10−2 (Li et al., 2014), much less than the heterozygosity of Yellow-breasted Bunting (0.46 × 10−2). Only seven crested ibis remained in the world in 1981, but the species has recovered to more than 5,000 individuals over the last 40 years (Valchuk, 2001). The rich genetic diversity in Yellow-breasted Bunting and the successful conservation efforts of crested ibis give us confidence that the Yellow-breasted Bunting has the genetic ability to recover.

Inbreeding can accelerate the extinction of threatened species because it can accelerate the spread of potential homozygous lethal mutations through inheritance (Frankham, 1995; Willis, 1993). In wild populations of the Glanville fritillary butterfly (Melitaea cinxia) and in laboratory cultures of fruit flies (Drosophila simulans), extinction risk increased significantly with increasing inbreeding (Saccheri et al., 1998; Wright et al., 2007). We did not find genetic signals of recent inbreeding in Yellow-breasted Buntings because the species has short LD extension and ROH (Figure 2) (Ceballos et al., 2018).

In Yellow-breasted Bunting, the average difference between homozygous and heterozygous frequency in LOF sites is significantly lower than that in synonymous sites (Figure 3), suggesting that the homozygous LOF sites are under purifying selection. The homozygous LOF mutations may reduce the survival rate. Although the Yellow-breasted Bunting may have potential homozygous lethal mutations, low levels of inbreeding may restrain the negative influence of these mutations. If we protect the remaining wild population, and conserve its habitat, the Yellow-breasted Bunting could recover, despite having potential homozygous lethal mutations. Low levels of inbreeding could reduce the negative influence of potential homozygous lethal mutations on survival. Certainly, we need to test the effects of the potential homozygous lethal mutations on fitness in the future.

The factors linked with rapid population decline

Natural selection, climate change, and human activities are three factors that affect the abundance of wild species (Feng et al., 2019; Murray et al., 2017; Zhao et al., 2012). Natural selection reduces genetic diversity through linkage selection (Buffalo, 2021); however, the LD extension in Yellow-breasted Bunting is short (Figure 2A). The short LD extension suggests natural selection has a lower probability of reducing genetic diversity through linkage selection in Yellow-breasted Bunting. The high genetic diversity of Yellow-breasted Bunting (Figure 1) also suggests natural selection may be not the main factor in prior population declines.

Climate change could act species’ abundance and cause extinction, and is in fact inferred as the one driver of long-term fluctuations in giant panda population size (Zhao et al., 2012). Demographic history suggests the Yellow-breasted Bunting experienced a bottleneck between 105 and 4 × 105 years ago (Figure 4) in the Middle Pleistocene. This period was characterized by climate oscillations, suggesting that climate change may have caused the bunting bottleneck. The Yellow-breasted Bunting began to decline about 100–200 generations ago (Figure 5). Anthropogenic global climate change has only been progressing since the late 19th century (Zanna et al., 2019), so it is unlikely that climate change directly caused the recent population decline of Yellow-breasted Bunting.

Interestingly, our study suggests that the time of decline in Yellow-breasted Bunting was much earlier than the prevailing view that the decline has mainly occurred since the 1980s (Zheng, 1956). Human activities have dramatically changed the global environment in many aspects as the industrial revolution, including climate warming, pollution, and fragmentation or loss of habitat (Steffen et al., 2016). All these processes could lead to population declines or extinction. As the industrial revolution, the global extinction rate is estimated to be at least 100 times greater than that in the fossil record (Hassan et al., 2005). Because the time that the Yellow-breasted Bunting began to decline is similar to the time of the industrial revolution, it is reasonable to infer that the decline may be directly or indirectly caused by human activities.

Human activities may be one driver of the population decline. The Yellow-breasted Bunting mainly breeds in wet meadows with scrub nearby. This includes wetlands in Russia. The species winters in South-east Asia and South China (11). In its place, cropland and pastures have increased by 358.1 million ha during 1850–1990 (Goldewijk, 2001). Hence, habitat loss caused by human activity could be the driver of the recent population decline. In addition, unsustainable trapping of the species is one driver of the recent decline. The Yellow-breasted Bunting has historically been trapped for food by the indigenous people in South China for at least 2,100 years (Yuren, 1996). China’s population has increased by about one billion people in the past 200 years (Poston and Yaukey, 2013), which could exacerbate anthropogenic pressure causing the Yellow-breasted Bunting decline. Furthermore, environmental pollution and climate change as the industrial revolution also could play an important role during the population decline (Gils et al., 2016; Kubelka et al., 2018; Zheng, 1956).

However, we still know little about the specific mechanisms of population decline in Yellow-breasted Bunting. We need to conduct many more studies, e.g., to determine the effects of heavy metal pollution on fitness, and the effects of habitat change coupled with climate warming, to infer the specific reason for population decline in the future.

Conservation implications

With a de novo genome assembly and population genomic data from 10 individuals, we show that Yellow-breasted Buntings have high genetic diversity and low levels of inbreeding, giving confidence in the potential to recover the species. Increasing population size over the long term may be the most important management goal, especially when the population has the mutation load and inbreeding is associated with inbreeding depression. It will be critical to increase its habitat size in future conservation efforts. Helping the local community to increase the paddy and wheat fields without the use of pesticides may be an effective method to increase the habitat size and population size of Yellow-breasted Bunting. It is critical to understand the fluctuations of the local population size and genetic variation with non-invasive sampling during conservation management.

On the other hand, protecting Yellow-breasted Bunting from illegal trapping is key to conservation. The Chinese government has placed Yellow-breasted Bunting in the highest protection classes (Ministry of Agriculture and Rural affairs of the People’s Republic of China, 2021), and imposed a ban on wildlife trade and consumption Zhang (2020). As a result, some positive signs have emerged. The number of criminal cases involving Yellow-breasted Bunting in the last two years has been significantly lower than before. A number of small breeding populations also recovered considerably in Russia (Heim et al., 2021). Ongoing and increased collaboration will help gain a better understanding of genomic status across its whole population and better advice for conservation efforts.

Limitations of the study

The current project used 10 samples to investigate the genetic status of Yellow-breasted Bunting, which did not cover all its populations. In future, evaluating the genetic status of all the breeding populations is necessary. What is more, we need to monitor the changes in genetic status for lasting years.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Peng Chen (capricorncp@163.com).

Materials availability

All the samples in the study could be shared with researchers who aim to conserve Yellow-breasted Buntings. Further requests for the samples will be fulfilled by the first author, Pengcheng Wang (wpc@njnu.edu.cn).

Method details

Sampling

On 15 September 2019, local policemen found Yellow-breasted Buntings were illegally trapped by mist net at Tangshan, East China, and brought them to Daqinghe Bird Rescue Center. Two buntings died soon after they arrived at the rescue center. These two Yellow-breasted Buntings were refrigerated immediately, and their tissues were sampled for de novo assembly and to annotate the genome. To analyze the population genomic status of Yellow-breasted Bunting, trace brachial vein blood was sampled from the remaining 10 rescued Yellow-breasted Bunting. After sampling and physical exam, the 10 healthy Yellow-breasted Buntings were released to the wild. All samples were kept in a −80 °C freezer in College of life sciences, Nanjing Normal University. Because the Yellow-breasted Bunting is a critically endangered species, the current project did not collect samples from all populations. All sample collection was done by an experienced veterinarian and sample information can be found in Table S3.

DNA/RNA extraction and sequencing

DNA and RNA extraction was completed using the phenol-chloroform method and trizol method, respectively. After passing a quality inspection, the DNA from fresh muscle was used to construct Illumina (Insert size was 350 bp, target data volume was 120 Gb) and Pacbio (target data volume was 120 Gb) genome libraries for genome assembly. The qualified RNA was used to construct Illumina (Insert size was 250–300 bp, target data volume was 10 Gb for each tissue) and Pacbio (target data volume was 30 Gb for each tissue) transcriptome libraries for genome annotation. The 10 qualified DNA samples from blood were used to construct Illumina DNA libraries for whole genome sequencing (Insert size was 350 bp, target data volume was 30 Gb for each sample). The Illumina DNA and RNA libraries were sequenced on Illumina HiSeq and Illumina NovaSeq 6000 platforms with paired-end 150 bp sequencing strategies, respectively. The Pacbio libraries were sequenced on Pacbio Sequel sequencing platform. The DNA/RNA extraction, genome/transcriptome library construction and sequencing were completed by the Novogene sequencing company, Beijing.

Genome assembly and evaluation

Nearly 186 Gb Pacbio subreads were used to assemble the genome in FALCON software (Kronenberg et al., 2019). After removing the reads whose N proportion was larger than 10%, and any low-quality reads (quality score less than 5), we retained reads with more than 30 bp, nearly 128 Gb of clean Illumina reads to polish the genome in the Nextpolish software (Hu et al., 2019). After polishing, we calculated the Scaffold N50 of the assembly. The Benchmarking Universal Single-Copy Orthologs (BUSCO) Aves dataset (Simão et al., 2015) and Core Eukaryotic Genes Mapping Approach (CEGMA) core gene dataset (Parra et al., 2008) were used to evaluate the completeness of the genome. To further evaluate the completeness of the genome and the sequencing uniformity, we mapped the clean Illumina reads to the assembled genome with BWA software (Li and Durbin, 2009) and calculated the mapping rate and coverage with Samtools (Danecek et al., 2021).

Genome annotation

For the Yellow-breasted Bunting genome, we performed repeat annotation, gene coding region annotation, and non-coding RNA annotation. To identify genome repeats, we used homology alignment and de novo predict methods. The Repbase database and RepeatMasker software was used for homology prediction (Bao et al., 2015; Chen, 2004). The LTR_FINDER, RepeatScout, and RepeatModeler with default parameters were used to de novo predict the transposable elements (Chen, 2004; Xu and Wang, 2007). Tandem repeats were predicted by the Tandem Repeats Finder (Benson, 1999). The DNA-level repeat identification was performed in RepeatMasker (Chen, 2004). We used the program tRNAscan-SE to predict the tRNAs (Lowe and Eddy, 1997). The software with Rfam database was used to identify ncRNAs (Griffiths-Jones et al., 2005; Nawrocki and Eddy, 2013).

After identifying genome repeats and non-coding RNA, we combined de novo prediction, homology sequences-based prediction, and RNA sequencing-assisted prediction to annotate the gene structure in MAKER pipeline (Cantarel et al., 2008). The protein sequences of turkey (Meleagris gallopavo), chicken (Gallus gallus), and zebra finch (Taeniopygia guttata) were downloaded from the Ensembl database (http://ftp.ensembl.org/pub/release-100/fasta/), and the protein sequences of Song Sparrow (Melospiza melodia) were downloaded from NCBI database (https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/013/398/205/GCA_013398205.1_ASM1339820v1/). These protein sequences were aligned to the Yellow-breasted Bunting genome using TblastN v. 2.2.26 with the E-value less than 1 × 10−5 (Camacho et al., 2009). We then used usedGeneWise v. 2.4.1 (Birney et al., 2004) to predict gene structure based on the matching proteins and homologous genome sequences. The Augustus v. 3.2.3 with Geneid v. 1.4, Genescan v. 1.0, GlimmerHmm v. 3.04 and SNAP (2013-11-29) were used to de novo predict gene structure (Blanco et al., 2007; Korf, 2004; Majoros et al., 2004; Stanke et al., 2006; Tiwari et al., 1997).

The raw Illumina transcriptome sequences were filtered in the Trimmomatic software with default parameters (Bolger et al., 2014). There were 13.46 Gb, 14.19 Gb, 14.22 Gb, and 11.85 Gb of clean Illumina transcriptome sequences of muscle, brain, heart and liver, respectively, used to correct the error of the corresponding Pacbio subreads. After error correction, the Pacbio transcriptome sequences from muscle, brain, heart and liver were each used to produce the consensus sequences in the isoseq pipeline (Gonzalez-Garay, 2016). TopHat v. 2.0.11 was used to align the RNA sequences from different tissues to the genome (Trapnell et al., 2009). The transcripts were then used as input for Cufflinks v. 2.2.1 to predict the gene structure (Trapnell et al., 2010). The consensus gene structures from the three methods were produced by EvidenceModeler v 1.1.1 with PASA (Program to Assemble Spliced Alignment) (Haas et al., 2008).

To assign the gene functions, we aligned the protein sequences to the Swiss-Prot using Blastp (E-value ≤ 1e-5) (Bairoch and Apweiler, 2000; Camacho et al., 2009). The InterProScan70 v. 5.31 was used to annotate the motifs and domains using the public databases, including Pfam, ProDom, PRINTS, SMRT, PANTHER and PROSITE (Mulder and Apweiler, 2007). The Gene Ontology (GO) IDs of each gene was assigned based on the corresponding InterPro entry (Mulder and Apweiler, 2007). We also used the NR database to predict the gene function in the BLAST software (Camacho et al., 2009) and clustered the genes to KEGG pathway (Kanehisa and Goto, 2000).

Identifying sex chromosome-linked scaffolds

Because sex chromosomes have different recombination rates compared to that of autosomes, we needed to identify the sex chromosome-linked scaffolds. We used the chicken genome as reference because it has relatively continuous sex chromosomes among published avian genomes. We aligned the draft Yellow-breasted Bunting genome to the chicken genome (ftp://ftp.ensembl.org/pub/release-98/fasta/gallus_gallus/dna/) with MUMmer v. 3.23 (Kurtz et al., 2004). After alignment, the longest consistent alignment was retained. The alignment results were filtered such that the minimum alignment length retained was 500, and minimum alignment identity was 50. After filtering, the position of the longest alignment delineated the location of each scaffold.

Read mapping and SNP calling

The raw Illumina sequences were trimmed using Trim Galore v. 0.6.6 (https://github.com/FelixKrueger/TrimGalore/releases/tag/0.6.6). We trimmed the adapter contents and removed five base pairs from both ends of each read. The FastQC v. 0.11.9 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to evaluate the quality of the clean reads. After trimming, the clean reads were mapped to the draft Yellow-breasted Bunting genome with a Burrows–Wheeler alignment-maximal exact matches (BWA-MEM) algorithm in the Burrows-Wheeler Alignment Tool (BWA) v. 0.7.17 (Li and Durbin, 2009). After alignment, we used Samtools v. 1.13 (Li et al., 2009) to calculate the coverage for each site of the genome, and then calculated the mean coverage of each sample. Samtools v. 1.13 was also used to estimate the mapping rate of each sample.

After mapping, the GATK v. 4.1.4.0 package (McKenna et al., 2010) was used to sort the reads and mark duplicates in the ‘bam’ file. The HaplotypeCaller, GenomicsDBImport, and GenotypeGVCFs tools in GATK v. 4.1.4.0 package were used to identify SNPs and indels according to the GATK best practices recommendations. vcftools v. 0.1.13 (Danecek et al., 2021) was used to extract autosome-linked SNPs and the quality metrics of each raw SNP. We used a density curve of each parameter as a reference to set filtering criteria as follows: QUAL <30.0, AC < 2, FS > 60.0, QD < 2.0, MQ < 40.0, SOR >3.0, MQRankSum < −4.0, ReadPosRankSum < −2.0, DP < 50, and DP > 361. Raw SNPs that fit any one of these filtering criteria were removed. The explanation of these parameters can be found from the GATK guide book (McKenna et al., 2010). The bi-allelic sites that on autosome were used for following analyses.

Family relationship inference

It is necessary to subsample non-closely related samples from the population genetic analysis. However, it is difficult to determine relationships from the pedigree of wild samples. Estimating kinship relationships among the samples based on genotype data can help confirm their relationships. We used KING v. 2.2.7 (Manichaikul et al., 2010) to estimate the relationship among samples. Following the KING manual, we set the threshold of the close relationship as kinship coefficient larger than 0.0442.

Estimating genetic diversity

To quantify genetic diversity of Yellow-breasted Buntings relative to other Passeriformes, we estimated the heterozygosity of the 10 Yellow-breasted Bunting and of 17 other Passeriformes birds. Heterozygosity refers to the number of heterozygous SNPs divided by the length of the genome without gaps. The heterozygosity of SNPs from each sample was extracted using vcftools v. 0.1.13 (Danecek et al., 2021). The 17 Passeriformes birds included Florida scrub jay (Aphelocoma coerulescens), Noisy scrubbird (Atrichornis clamosus), Chestnut-collared longspur (Calcarius ornatus), Seychelles magpie robin (Copsychus sechellarum), Ribbon-tailed drongo (Dicrurus megarhynchus), Painted honeyeater (Grantiella picta), Akiapola'au (Hemignathus wilsoni), Philippine fairy-bluebird (Irena cyanogastra), Kirtland's warbler (Setophaga kirtlandii), Loggerhead shrike (Lanius ludovicianu), Bali myna (Leucopsar rothschildi), Yellowhead (Mohoua ochrocephala), Velvet myiagra (Myiagra hebetior), Inaccessible island finch (Nesospiza acunhae), Stitchbird (Notiomystis cincta), White-necked rockfowl (Picathartes gymnocephalus), and Golden-crowned babbler (Sterrhoptilus dennistouni). The reference genomes and the fragment genome libraries of these 17 birds were downloaded from the National Center for Biotechnology Information search database (NCBI) (Feng et al., 2020). We used same pipeline as for SNPs of Yellow-breasted Bunting to locate and filter SNPs of these 17 birds.

Estimating inbreeding pattern

To infer the inbreeding status of Yellow-breasted Bunting, we calculated the linkage disequilibrium (LD) pattern and runs of homozygosity (ROH). Haploview v. 4.2 was used to calculate the correlation coefficient (r2) of any pair of SNPs on autosomes (Barrett et al., 2004). The greater the R2 value, the greater the degree of linkage disequilibrium. We set the maximum distance between any pair of SNPs to 500 kb and the minor allele frequency to 0.01. To plot the LD decay curve, we calculated the average R2 in each 100 bp window. ROH always form due to inbreeding; the offspring inherit the identical-by-decent haplotype from parents (Ceballos et al., 2018). Thus, the number and length of ROH in the populations can help us infer the inbreeding pattern (Ceballos et al., 2018). In inferring the number and length of ROH, we first used vcftools v. 0.1.13 (Danecek et al., 2021) to convert vcf file to PLINK ped format for each scaffold and then used plink v. 1.90 (Howrigan et al., 2011) to infer ROH with the following parameters: ROH of at least 50 SNPs, at least 3 heterozygous SNPs, at least one SNP per 50 kb, a minimum of 300 kb length for an ROH, maximum distance between adjacent SNPs of 1000 kb per ROH, each scanning window having length of 50 SNPs, and 10 missing SNPs. Recombination can break the ROH. Thus, long and short ROH combined can reflect ancestral and recent inbreeding, respectively (Ceballos et al., 2018). We calculated the number and the total length of short ROH (<1Mb) and long ROH (>1.5 Mb) of each sample. We also calculated the percentage of the two kinds of ROH within full autosomes. Kolmogorov-Smirnov comparison was used to test whether the percentage has significant difference between two ROH groups.

Estimating potential inbreeding depression

To infer whether the Yellow-breasted Bunting genome contains signs of historical inbreeding depression, we used SnpEff v. 5.0c to identify the synonymous, missense, and lose of function (LOF) mutations in the 10 samples (Cingolani et al., 2012). After that, we calculated the average difference between the frequency of homozygous sites and heterozygous sites in 1000 random SNPs within the coding region (relative frequency of homozygous sites), according to the following formula:

fhom is the frequency of homozygous mutations in each site. fhet is the frequency of heterozygous mutations in each site. C is the set of synonymous, missense, and LOF mutations. n is the number of different sets of mutations among the 1000 random SNPs. To test whether the relative frequency of homozygous sites is significantly different between missense and synonymous mutations, and between LOF and synonymous mutations, we performed 100 bootstrap replications of the coding sites. If the difference was not significant, this indicated the missense or LOF mutations were similar to neutral mutations; if the relative frequency of homozygous sites of missense or LOF is significantly larger than that of the synonymous mutations, then the Yellow-breasted Bunting likely has homozygous non-lethal mutations; if the relative frequency of homozygous sites of missense or LOF is significantly smaller than that of synonymous mutations, the Yellow-breasted Bunting has potential homozygous lethal mutations. We used a Kolmogorov-Smirnov comparison to test whether the difference was significant. To explore the potential influence of deleterious mutations on biological function, we performed KEGG pathway enrichment analysis for the genes that harbor LOF mutations in KOBAS (Bu et al., 2021).

Reconstructing yellow-breasted bunting demographic history

We aimed to detect the historical factors that caused Yellow-breasted Bunting’s rapid decline. To do so, we used three methods to infer the species’ effective population size (Ne) over time: Pairwise Sequentially Markovian Coalescent (PSMC) modeling (Li and Durbin, 2011), Stairway plots (Liu and Fu, 2015), and fastsimcoal (Excoffier et al., 2013; Excoffier and Foll, 2011). PSMC models use a hidden Markov model and local density of heterozygotes to infer the time of the most recent common ancestor between two alleles in a diploid genome sequence (Li and Durbin, 2011). We followed the instruction of PSMC to analysis the historical effective population size fluctuations, but we adjusted some parameters as follows: minimum accepted mapping quality for an alignment was 1, minimum read depth was 6, maximum read depth was 46, maximum number of iterations was 30, and the initial theta/rho ratio was 5. While PSMC models provide detailed information on the fluctuation of Ne, stairway plots can produce more accurate inference for recent Ne changes (Liu and Fu, 2015). Hence, we also used stairway plots to infer the historical Ne fluctuations using the SNP frequency spectrum (SFS) (Liu and Fu, 2015). The ANGSD v. 0.910 was used to produce the folded SFS from the bam files (Korneliussen et al., 2014). The folded SFS file was used as input file for stairway plot v. 2.1.1 (Liu and Fu, 2015).

To simulate the demographic history of Yellow-breasted Bunting in recent centuries, we used fastsimcoal (Excoffier and Foll, 2011), a complementary approach to PSMC and stairway plots. Fastsimcoal is an ancestral recombination graph (ARG) model based simulation program that can simulate complex demographic scenarios with SFS and then select the best fit model based on the Akaike Information Criterion (AIC) and the distribution of the likelihood of the models (Excoffier and Foll, 2011). Because linkage between any two sites will reduce the recombination rates, we filtered the linkage SNPs before running fastsimcoal. The plink software was used to filter the SNPs that their linkage disequilibrium coefficient larger than 0.2 in 1000 kb windows (Chang et al., 2015). After filtering, the 14,808 no-missing bi-allelic sites were used to produce folded SFS file in easySFS (https://github.com/isaacovercast/easySFS). The new SFS file was used for fastsimcoal26 (Excoffier and Foll, 2011).

We devised four demographic models based on the results from PSMC and the stairway plot. For each model, we performed 100 replicates of the optimum procedure, each with 100 expectations conditional maximization (ECM) cycles and 1,000,000 simulations per cycle. The replicate with the highest likelihood was considered the best model of each demographic model. To identify the best fit model, we re-ran each demographic model 100 times with the best parameter from the initial 100 replications. If the average likelihood of one model is significant different with other models, and its AIC value is less than other models, the model was considered as the best fit demographic model. The Welch two sample t-test was used to determine whether the likelihood between two models is significant. The AIC value was calculated based on the likelihood value and the number of demographic parameters.

The four models included (the schematic diagram is shown in Figure S4): (a) Bottleneck and stable model (BS), in which Yellow-breasted Bunting has a constant population size after a bottleneck; (b) Bottleneck and decline model (BD), in which after experiencing a bottleneck, Yellow-breasted Bunting had a constant population and then declined continuously; (c) Bottleneck, which had stable ancestry and early decline model (BSED), the model is similar to model BD, but with time of population decline between 100 and 200 generations ago; (d) Bottleneck, stable ancestry and recent decline model (BSRD), similar to the BD model, but with population decline at 6–30 generations ago. The 200 generations ago of the Yellow-breasted Bunting was probably the beginning of the industrial revolution. In inferring the historical Ne fluctuations, we set the mutation rate as 3.3 × 10−9 per site per generation (Zhang et al., 2014) and assumed the generation time to be two years. After we identified the best fit model, we did 1000 bootstraps to estimate the demographic parameters and their confidence interval.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
The blood of Yellow-breasted Buntings (Emberiza aureola).	Daqinghe Bird Rescue Center, Tangshan, Hebei, China.	N/A
Deposited data
The reference genomes and the fragment genome libraries of Florida scrub jay (Aphelocoma coerulescens)	National Center for Biotechnology Information search database (NCBI)	SAMN12253992
The reference genomes and the fragment genome libraries of noisy scrubbird (Atrichornis clamosus)	National Center for Biotechnology Information search database (NCBI)	SAMN12253955
The reference genomes and the fragment genome libraries of chestnut-collared longspur (Calcarius ornatus)	National Center for Biotechnology Information search database (NCBI)	SAMN12253923
The reference genomes and the fragment genome libraries of Seychelles magpie robin (Copsychus sechellarum)	National Center for Biotechnology Information search database (NCBI)	SAMN12254010
The reference genomes and the fragment genome libraries of ribbon-tailed drongo (Dicrurus megarhynchus)	National Center for Biotechnology Information search database (NCBI)	SAMN12253804
The reference genomes and the fragment genome libraries of painted honeyeater (Grantiella picta)	National Center for Biotechnology Information search database (NCBI)	SAMN12253950
The reference genomes and the fragment genome libraries of Akiapola’au (Hemignathus wilsoni)	National Center for Biotechnology Information search database (NCBI)	SAMN10867508
The reference genomes and the fragment genome libraries of Philippine fairy-bluebird (Irena cyanogastra)	National Center for Biotechnology Information search database (NCBI)	SAMN12253784
The reference genomes and the fragment genome libraries of Kirtland’s warbler (Setophaga kirtlandii)	National Center for Biotechnology Information search database (NCBI)	SAMN12253801
The reference genomes and the fragment genome libraries of loggerhead shrike (Lanius ludovicianu)	National Center for Biotechnology Information search database (NCBI)	SAMN12253815
The reference genomes and the fragment genome libraries of Bali myna (Leucopsar rothschildi)	National Center for Biotechnology Information search database (NCBI)	SAMN12253829
The reference genomes and the fragment genome libraries of yellowhead (Mohoua ochrocephala)	National Center for Biotechnology Information search database (NCBI)	SAMN12253963
The reference genomes and the fragment genome libraries of velvet myiagra (Myiagra hebetior)	National Center for Biotechnology Information search database (NCBI)	SAMN12253791
The reference genomes and the fragment genome libraries of Inaccessible Island finch (Nesospiza acunhae)	National Center for Biotechnology Information search database (NCBI)	SAMN12254004
The reference genomes and the fragment genome libraries of stitchbird (Notiomystis cincta)	National Center for Biotechnology Information search database (NCBI)	SAMN12253957
The reference genomes and the fragment genome libraries of white-necked rockfowl (Picathartes gymnocephalus)	National Center for Biotechnology Information search database (NCBI)	SAMN12253907
The reference genomes and the fragment genome libraries of golden-crowned babbler (Sterrhoptilus dennistouni)	National Center for Biotechnology Information search database (NCBI)	SAMN12253785
chicken genome	Ensembl	ftp://ftp.ensembl.org/pub/release-98/fasta/gallus_gallus/dna/
Raw sequencing data	This paper	NGDC: PRJCA009403
Software and algorithms
FALCON	Kronenberg et al. (2019)	https://github.com/PacificBiosciences/FALCON
Nextpolish	Hu et al. (2019)	https://github.com/Nextomics/NextPolish
RepeatMasker v. 4.1.2	N/A	http://www.repeatmasker.org/RepeatModeler/
RepeatModeler v. 2.0.3	N/A	http://www.repeatmasker.org/RepeatModeler/
SNAP v. 20131129	Korf (2004)	http://korflab.ucdavis.edu/software.html
Augustus v. 3.4.0	Stanke et al., (2006)	http://bioinf.uni-greifswald.de/augustus/
Maker v. 3.1.4	Cantarel et al., (2008)	http://weatherby.genetics.utah.edu/MAKER/wiki/index.php/Main_Page
MUMmer v. 3.23	Kurtz et al. (2004)	http://mummer.sourceforge.net/
Burrows-Wheeler Alignment Tool (BWA) v. 0.7.17	Li and Durbin (2009)	http://bio-bwa.sourceforge.net/
Samtools v. 1.13	Li et al. (2009)	http://samtools.sourceforge.net/
GATK v. 4.1.4.0	McKenna et al. (2010)	https://gatk.broadinstitute.org/hc/en-us/sections/360007279452-4-1-4-0
vcftools v. 0.1.13	Danecek et al., 2021	http://vcftools.sourceforge.net/
KING v. 2.2.7	Manichaikul et al. (2010)	https://www.kingrelatedness.com/
Haploview v. 4.2	Barrett et al. (2004)	https://www.broadinstitute.org/haploview/haploview
PLINK	Howrigan et al. (2011)	https://www.cog-genomics.org/plink/
SnpEff v. 5.0	Cingolani et al. (2012)	http://pcingola.github.io/SnpEff/
Pairwise Sequentially Markovian Coalescent (PSMC)	Li and Durbin (2011)	https://github.com/lh3/psmc
Stairway plots	Liu and Fu (2015)	https://github.com/xiaoming-liu/stairway-plot-v22
fastsimcoal26	Excoffier and Foll (2011)	http://cmpg.unibe.ch/software/fastsimcoal26/
easySFS	https://github.com/isaacovercast/easySFS	https://github.com/isaacovercast/easySFS
ANGSD v. 0.910	Korneliussen et al. (2014)	http://www.popgen.dk/angsd/index.php/ANGSD