https://orcid.org/0000-0003-2051-2516
Abstract
The crested ibis (Nipponia nippon) is endangered worldwide. Although a series of conservation measures have markedly increased the population size and distribution area of these birds, the high mortality of embryos and nestlings considerably decreases the survival potential of this bird species. High-throughput sequencing technology was utilized to compare whole genomes between ten samples from dead crested ibises (including six dead embryos and four dead nestlings aged 0–45 days) and 32 samples from living birds. The results indicated that the dead samples all shared the genetic background of a specific ancestral subpopulation. Furthermore, the dead individuals were less genetically diverse and suffered higher degrees of inbreeding compared with these measures in live birds. Several candidate genes (KLHL3, SETDB2, TNNT2, PKP1, AK1, and EXOSC3) associated with detrimental diseases were identified in the genomic regions that differed between the alive and dead samples, which are likely responsible for the death of embryos and nestlings. In addition, in these regions, we also found several genes involved in the protein catabolic process (UBE4A and LONP1), lipid metabolism (ACOT1), glycan biosynthesis and metabolism (HYAL1 and HYAL4), and the immune system (JAM2) that are likely to promote the normal development of embryos and nestlings. The aberrant conditions of these genes and biological processes may contribute to the death of embryos and nestlings. Our data identify congenital factors underlying the death of embryos and nestlings at the whole genome level, which may be useful toward informing more effective conservation efforts for this bird species.
Introduction
The crested ibis (Nipponia nippon), a species of bird, endangered worldwide. It is classified as order Pelecaniformes, family Threskiornithidae (Birdlife International 2018). This bird was historically distributed widely across East Asia with numerous flocks and individuals (Ding 2004). During the 19th to mid-20th century, the number of crested ibises suffered a rapid decrease because of climate change, habitat loss, overhunting, and agrochemical use (Ding 2004; Li et al. 2014). In 1981, after the remaining five wild crested ibises in Japan were captured and subjected to captive breeding, this bird was considered to be nearly extinct in the wild (Ding 2010). Although these birds had not been observed for over 17 years in China, a 3-year intensive investigation eventually identified two additional breeding pairs and three chicks in Yangxian County (Shaanxi Province, China) in 1981 (Ding 2010; Li et al. 2014). Subsequently, a series of conservation measures were immediately taken by the Chinese government, such as habitat and foraging site conservation, legislation for a ban on hunting, and the establishment of a special protection agency (Shi and Cao 2001; Wu et al. 2017). These measures effectively increased the number of individuals and distribution areas of the created ibis. Over the past 30-plus years, the IUCN Red List categorization of the crested ibis has been upgraded from critically endangered (CE) to endangered (EN), with the population and distribution areas of the crested ibis improving from only one remnant population with seven birds in one area to more than ten areas with nearly 2,500 birds in three countries (supplementary fig. S1, Supplementary Material online) (Ding 2010; Wajiki et al. 2015; Li et al. 2016; Birdlife International 2018).
Even though the population size and distribution of the crested ibis have been increasing steadily, this species still faces many survival challenges. The limited number of founder individuals renders it susceptible to the loss of genetic variation, fixation of deleterious mutations, and inbreeding depression, which can negatively impact viability, disease resistance, fertility, embryo survival, and number of litters (Soulé 1980; Lehmkuhl 1984; Falconer 1989; Lande 1995, 1998; Fischer et al. 2003). Accordingly, the mean death rates of the embryos and nestlings aged 0–42 days of the crested ibis were ∼37% in the Shannxi (China) Crested Ibis Feeding and Breeding Center during 1995–2000 and the mean death probability of the yearlings (0–1 years old) of a reintroduced crested ibis population in Ningshan (China) during 2008–2015 was ∼61.6% (Xi et al. 2001; Li et al. 2018). In comparison, the mean death rate of nestlings aged 0–30 days of 12 species of birds in captive populations was 22.73%, ranging from 5.48% to 46.46% (Boakes et al. 2007). The high mortality of yearlings has significantly diminished the survival probability of the crested ibis population; hence, effective protection of yearlings is considered critical to population survival (Li et al. 2018). Conversely, previous research related to the crested ibis focused on external factors such as viruses, bacteria, and parasites that influence survival, and incorporated microsatellite markers, major histocompatibility complex loci, and defensins that only provide limited genomic information for the development of particular measures for species conservation (Zhang et al. 2004, 2006, 2015; He et al. 2006; Xi et al. 2007; Chen et al. 2013, 2015; Urano et al. 2013; Hao et al. 2014; Lan 2014; Taniguchi et al. 2014). The rapid developments in high-throughput sequencing technology increased the research in genetic disorders and multifactorial etiologies (Hennekam and Biesecker 2012). Although numerous human genetic disorders have been studied via whole genome resequencing, animal hereditary diseases are rarely investigated using this technology (Dey et al. 2018; Estrada et al. 2018; Fahey et al. 2018; Tang et al. 2018; Vardarajan et al. 2018; Wang et al. 2018). Notably, as the whole crested ibis genome was assembled by Li et al. (2014), it is now feasible to study the intrinsic factors that affect embryonic and nestling death of the crested ibis (Li et al. 2014).
The captive crested ibis population of Deqing County (Zhejiang, China), established in 2008, constitutes the first southern population of the crested ibis in China and has five breeding couples introduced from the captive crested ibis population of Zhouzhi County (Shaanxi, China) as founders. In order to improve the genetic diversity of this population, ten and eight novel individuals were introduced from the captive population of Japan to this population in 2009 and 2011, respectively, and then four novel individuals from the captive population of Beijing were introduced into this population in 2014. By the end of the breeding season in 2017, after 9 years of assisted breeding with artificial incubation, the size of this captive population has gradually increased to 209 individuals. However, the death rate of embryos and nestlings aged 0–45 days of this captive population is high, ranging from 0% to 61.8% from 2009 to 2017 with a mean death rate of 38.8%.
To better understand potential genetic factors underlying the high embryo and nestling mortality in the Deqing population, in this study, we sequenced the whole genome of ten samples from six dead embryos and four dead nestlings aged 0–45 days, together with 32 samples from surviving crested ibises that were all incubated and fed under the same experimental conditions. These whole-genomes were then compared with explore genomic differences between the two groups. The findings are expected to substantially contribute to increasing the size and viability of the crested ibis population and will provide a better understanding at the whole genome scale of factors influencing the conservation of this endangered species.
Materials and Methods
Sample Collection
All samples used in this research were collected from the Xiazhu Lake Crested Ibis Breeding Base in Deqing County (Zhejiang, China) from 2009 to 2017. During the crested ibis breeding season, we collected eggs every morning for artificial incubation and recorded the exact time that the eggs started incubation. The hatching conditions and brooding conditions used were the same as those reported in a previous study (Xi et al. 2001). After 28 days, if the chick had still not pecked through the shell and a faint cheeping could not be detected, the shell was carefully peeled away and the viability of the chick was determined. If the embryo was dead, a tissue sample was collected and stored at −20°C. We also collected a tissue sample if the nestling had died from unknown reasons within 45 days posthatching, and stored it at −20°C. The peripheral venous blood samples from living birds were drawn via wing venipuncture when birds reached an age of 45 days, and stored at −20°C. Samples were collected under the supervision of the Animal Ethics Committee of Zhejiang University and permission was acquired from the Xiazhu Lake Crested Ibis Breeding Base in Deqing County (Zhejiang, China). We tried to collect as many samples of individuals that died at an early age as possible; however, owing to the problem of sample collection and genomic degradation, ten samples including six dead embryos and four dead nestlings aged 0–45 days (DEN: Dead embryos and nestlings) were amenable for whole genome resequencing. Of those, five dead embryos and two dead nestlings were collected in 2017, one dead nestling was collected in 2016, one dead embryo was collected in 2014, and one dead nestling was collected in 2013. About 32 samples from living crested ibises (AL: Alive individuals) were accessible for whole genome resequencing, including 12 birds that were born in 2017, 6 birds that were born in 2016, 2 birds that were born in 2009, and 5 founder breeding couples and 2 second spouses (the first spouses of 2 founder breeding couples died). Pedigree information on this captive population is missing (supplementary table S1, Supplementary Material online).
DNA Extraction and Sequencing
Genomic DNA was extracted using the DNeasy Tissue & Blood Kit (Generay, China). DNA degradation and contamination were detected by 1% agarose gel electrophoresis and its concentration was measured using the Qbit DNA Assay Kit with a Qbit 2.0 Fluorometer (Life Technologies, Carlsbad, CA). Sequencing libraries were constructed using the TruSeq Nano DNA HT Sample preparation Kit (Illumina, San Diego, CA) following manufacturer’s instructions. The generated libraries, with insert sizes of ∼350 bp, were sequenced using the Illumina HiSeq X Ten system to generate 150-bp paired-end reads at the Beijing Novogene Co., Ltd (Beijing, China). Approximately 1,386.27 Gb of raw data were generated in total. The Q30 ranged from 88.50% to 94.47% and GC content comprised 43.35–47.49% for all samples (supplementary table S1, Supplementary Material online).
Sequence Quality Checking and Read Mapping
We removed low quality paired reads based on the following criteria: 1) ≥10% unidentified nucleotides, 2) >50% bases having phred quality <5, 3) >10 nt aligned to the adapter, allowing ≤10% mismatches, and 4) putative PCR duplicates (i.e., read1 and read2 of two paired-end reads were completely identical). These low-quality reads were mainly caused by base-calling duplicates and adapter contamination. The cleaned data accounted for ∼99.26% of the raw data and the amount of high-quality sequence data of all samples ranged from 30.32 to 47.41 Gb (supplementary table S1, Supplementary Material online). We mapped the remaining high-quality reads to the N. nippon reference genome (GCF_000708225.1) using BWA-MEM (0.7.10-r789) with default parameters (Li and Durbin 2009; Li et al. 2014). We used SAMtools1.8 (Li et al. 2009), Picard2.18.7 (http://broadinstitute.github.io/picard/), and GATK4.0 (DePristo et al. 2011) (https://software.broadinstitute.org/gatk/) successively for sorting and merging the preliminary alignment results, assigning read group information, and performing local realignment, respectively. The cleaned reads of all samples appropriately mapped to the crested ibis reference genome. The effective sequencing depth exceeded 22.70× (22.70–33.98×) and the genome coverage surpassed 95.24% (95.24–96.92%) for all samples (supplementary table S1, Supplementary Material online).
Variant Calling
We used Picard2.18.7 and GATK4.0 (DePristo et al. 2011) with the HaplotypeCaller tool and followed the standard pipeline to perform SNP calling. GATK4.0 was used to identify insertions and deletions (InDels) with a fragment size <50 bp and it was also used to perform variant quality score recalibration, which utilized the top10 of the best-quality InDels as the training set to model InDels. We used the linear regression function in R3.5.1 to test for assessing differences in the number of SNPs and InDels between the AL and DEN samples, because regression has been shown to be robust to violations of the normality assumption (Knief and Forstmeier 2018).
Runs of Homozygosity
Population Genetic Structure
The program ADMIXTURE1.3.0 (Alexander et al. 2009) was implemented to infer the population genetic structure. The number of clusters was varied from K = 2 to K = 5 and the lowest cross-validation error was used to determine the optimal K value. The neighbor-joining method was performed in VCF-kit (http://vcf-kit.readthedocs.io/en/latest/) to reconstruct a phylogenetic tree of the 42 samples, which was visualized using FigTree v1.4.3. EIGENSOFT v6.0.1 was utilized to perform the principal component analysis and the Tracey–Widom test was conducted to ascertain the significance of eigenvectors (Patterson et al. 2006).
Identification of Genomic Regions with Differences between AL and DEN
In order to identify genomic regions with differences between the AL and DEN samples, we calculated the difference in sequence diversity (as the ratio ln(πAL/πDEN)), the population differentiation statistic (FST), and Tajima’s D values between AL and DEN samples in a genome-wide sliding window approach with a window size of 30 kb and a step size of 15 kb. The windows with simultaneously the lowest or highest 5% ln(πAL/πDEN) ratios (i.e., ln(πAL/πDEN) < −0.4645 or ln(πAL/πDEN) > 0.60822) and the highest 5% FST values (i.e., FST > 0.134475) were selected as representative regions (part A and part B) with differences between AL and DEN. These regions were further confirmed by testing whether there was a significant difference in values of Tajima’s D between AL and DEN (Li et al. 2013; Qi et al. 2013; Qiu et al. 2015). We used the linear regression function in R3.5.1 to test for differences in values of Tajima’s D of the representative regions with difference between AL and DEN. KOBAS 3.0 was used to perform functional enrichment analysis with the candidate genes located in the genomic outlier regions using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations (Wu et al. 2006; Xie et al. 2011).
Results
SNPs and InDels
Using SAMtools, 1.29 million SNPs were identified, with the majority (811,567; 62.70%) being located in intergenic regions, and 20,691 located in exonic regions (supplementary table S2, Supplementary Material online). The average number of SNPs across all individuals was 766,748. The number of SNPs per individual was significantly lower in the DEN samples (estimate ± SE: 748,737 ± 8,991) than in the AL samples (estimate ± SE: 772,377 ± 5,026) (P = 0.03, fig. 1A and supplementary table S3, Supplementary Material online). The sequence diversity (π) of all samples was 4.65 × 10−4, and the sequence diversity of DEN (4.60 × 10−4) was also lower than that of AL (4.64 × 10−4). We obtained a total of 114,794 insertions and 102,935 deletions among all samples. The average number of InDels across all individuals was 123,085, and the number of InDels per individual in the DEN samples (estimate ± SE: 121,046 ± 1,103) was significantly lower than that in the AL samples (estimate ± SE: 123,722 ± 617) (P = 0.04, fig. 1B and supplementary table S3, Supplementary Material online).
—Genetic variation in the crested ibis genome. (A) The number of SNPs between the alive (AL) and dead nestling and embryo (DEN) samples. (B) The number of InDels between the AL and DEN samples. Black dots indicate outliers.
Runs of Homozygosity
The mean number of ROHs across all individuals was 2,979 and the number of ROHs per individual in the DEN samples (estimate ± SE: 3,080 ± 94) was significantly higher than that in the AL samples (estimate ± SE: 2,711 ± 52) (P = 0.0014, fig. 2A and supplementary table S4, Supplementary Material online). The mean value of the total length of ROHs across all individuals was 219,219,330 bp, and the total length of ROHs per individual in the DEN samples (estimate ± SE: 258,397,906 ± 10,482,167 bp) was significantly higher than that in the AL samples (estimate ± SE: 206,976,024 ± 5,859,710 bp) (P = 0.0001, fig. 2B and supplementary table S4, Supplementary Material online). The average Froh was 0.20 for all samples. The Froh of the DEN group (estimate ± SE: 0.231 ± 0.009) was significantly higher than that of the AL group (estimate ± SE: 0.184 ± 0.005) (P = 0.0001, fig. 3 and supplementary table S4, Supplementary Material online).
—Runs of homozygosity (ROHs) in the crested ibis genome. (A) The number of ROHs per individual between the alive (AL) and dead nestling and embryo (DEN) samples. (B) The total length of ROHs per individual between the AL and DEN samples. Black dots indicate outliers.
—Genomic inbreeding coefficient (Froh) of the alive (AL) and dead embryo and nestling (DEN) samples. Only autosomal sequences were analyzed; samples were ordered by their values of Froh.
Population Genetic Structure of the Captive Crested Ibis Population of Deqing
We inferred that the captive crested ibis population of Deqing was derived from three ancestral subpopulations by running the ADMIXTURE1.3 program from K = 2 to K = 5 and calculating the cross-validation error (fig. 4A and B). We also found that all DEN samples possessed the genetic background from ancestral subpopulation 1 (fig. 4C). We grouped the 42 samples on the basis of the contribution percentage of the three ancestral subpopulations (fig. 5A and supplementary table S5, Supplementary Material online). The genomic components of three samples, all of which were DEN, were completely derived from ancestral subpopulation 1, whereas eight and two samples (all AL) were entirely derived from ancestral subpopulation 2 and 3, respectively. The results of neighbor joining tree and principal component analysis were consistent with the results of population genetic structure analysis with K = 3 (fig. 5B and C).
—Population genetic structure of the captive crested ibis population of Deqing. (A) Population structure plots with K = 2 to K = 5. (B) Broken line graph of cross-validation error. (C) Samples grouped according to dead and alive state while K = 3.
—Phylogenetic analyses of the captive crested ibis population of Deqing. (A) Statistics of samples based on composition of the three ancestral subpopulations (e.g., samples labeled with ancestral subpopulation 1_3_2 indicates that the genomic components of these samples were derived from ancestral subpopulations 1, 3, and 2, with the order of the genomic components being ancestral subpopulation 1 > ancestral subpopulation 3 > ancestral subpopulation 2). (B) Neighbor joining tree based on the whole genome single nucleotide polymorphisms of the crested ibis. Scale bar represents the level of similarity. The tree is colored according to the ancestor subpopulation composition of each sample from the result of ADMIXTURE; branches with red dots represent samples belonging to the dead embryo and nestling group. (C) Results of principal components analysis of 42 crested ibises. The percentage of variance explained by PC1, PC2, and PC3 is 24.23%, 23.85%, and 23.29%, respectively.
Identification of Genomic Regions That Differ between the AL and DEN Samples
The windows with simultaneously the lowest 5% ln(πAL/πDEN) ratios and the highest 5% FST values, and the windows with simultaneously the highest 5% ln(πAL/πDEN) ratios and the highest 5% FST values were selected as part A and part B of representative regions with differences between AL and DEN, respectively (fig. 6 and supplementary fig. S2A and B, Supplementary Material online). These regions were further confirmed by significantly different values in Tajima’s D between the AL and DEN (P < 2.2e − 16 and P < 2.2e − 16 of part A and part B, respectively). We identified 858 (15.78 Mb in total) and 269 windows (4.65 Mb in total) in the part A and part B of the representative regions with difference between AL and DEN, respectively, with these windows containing 107 and 31 nonsynonymous SNPs in 117 and 57 protein-coding genes, respectively. Of those, 56 and 21 genes contained the nonsynonymous SNPs, respectively. To obtain further information regarding the functional roles of these genes in the crested ibis genome, we used KOBAS to perform GO and KEGG pathway analysis.
—Genomic regions with difference between the alive (AL) and dead embryo and nestling (DEN) samples of crested ibis. Distribution of the population statistic (FST) values and the ln(πAL/πDEN) ratios of 30-kb windows with 15-kb steps. Blue and red dots highlight part A and part B of the representative regions with difference between AL and DEN, respectively.
A total of 110 GO terms and 4 KEGG pathways were significantly enriched in part A of the representative regions at P < 0.05(supplementary tables S6 and S7, Supplementary Material online). Among these, five significantly enriched GO terms (P < 0.05, GO: 0030163; GO: 0006511; GO: 0019941; GO: 0051603; GO: 0044257) included four genes (UBE4A, ANAPC2, USP40, and LONP1) that were related to the protein catabolic process. In addition, ADAM11 (with one nonsynonymous SNP) and JAM2, which fall within the same GO term of binding (GO: 0005488), may play a role as tumor suppressors. The four significantly enriched pathways are associated with lipid metabolism, glycan biosynthesis and metabolism, and the immune system, respectively. Specifically, ACOT1 (with three nonsynonymous SNPs) pertains to the pathway of lipid metabolism; HYAL1 (with one nonsynonymous SNP) and HYAL4, which can catabolize hyaluronic acid, belong to the pathway of glycan biosynthesis and metabolism, and JAM2, in the immune system pathway, plays a crucial role in the migration of natural killer cells.
Moreover, 58 GO terms and 2 KEGG pathways were significantly enriched in part B of the representative regions at P < 0.05 (supplementary tables S8 and S9, Supplementary Material online). Among the GO terms, four genes were associated with diseases that likely result in lethality: KLHL3, SETDB2 (with one nonsynonymous SNP), TNNT2, and PKP1, all within the GO term “intracellular” (GO: 0005622). In addition, AK1, which pertains to the pathway of thiamine metabolism and EXOSC3, which belongs to the pathway of RNA degradation, are also associated with severe, potential lethal diseases.
Discussion
Population Genetic Structure of the Crested Ibis
Our findings suggest that the captive crested ibis population of Deqing derived from three ancestral subpopulations. This population was founded by five breeding couples that were introduced from Zhouzhi County in 2008, merged with ten and eight novel individuals from the captive population of Japan in 2009 and 2011, respectively, then acquired four individuals from a captive population of Beijing in 2014 to improve its genetic diversity. This population thus represents the entire global population of the crested ibis to some extent. Thus, all living individuals of the global crested ibis population may be descendants of the three ancestral subpopulations. However, individuals of the Sado captive crested ibis population were reported to derive from only two ancestral subpopulations (Urano et al. 2013). This outcome may have been reached because the Sado captive population was founded by only five individuals from China and the results were obtained using microsatellite markers with limited polymorphism. In addition, we also found that all DEN individuals possessed the genetic background of a specific subpopulation, with the genomic components of three DEN samples being completely inherited from this ancestral subpopulation. This result suggested that the genome of this ancestral subpopulation may contain some recessive deleterious mutations, with homozygosity of these mutations potentially constituting one of the reasons causing the death of affected embryos and nestlings. In subsequent breeding projects, we should try to introduce more novel individuals from other populations into the captive crested ibis population of Deqing to increase its genetic diversity. In addition, we should try to collect accurate pedigree information by careful observation, and perform paternity testing on newborn birds. Based on the complete pedigree data, we would thus be able to minimize inbreeding with appropriate human intervention and reduce the risks of species extinction.
Inbreeding Depression as a Cause of Embryo and Nestling Death
Mating among individuals or organisms that have close affinities with each other through common ancestry is termed inbreeding, which often leads to inbreeding depression, that is, a reduction in fitness (Falconer 1989; Charlesworth and Willis 2009). Research from birds and mammals in the wild revealed that inbreeding depression has a marked impact on birth weight, survival, proliferation, and disease resistance (Keller and Waller 2002). Furthermore, by evaluating inbred wild species under natural conditions, it was discovered that wild animals commonly displayed moderate to high degrees of inbreeding depression in fitness traits (Crnokrak and Roff 1999). Moreover, the disease prevalence of human gout, depression, peptic ulcer, and cancer increased with higher inbreeding coefficients, and 23–48% of disease morbidity can be ascribed to inbreeding (Rudan et al. 2003).
Researchers have found that the crested ibis shows a high inbreeding coefficient and low genetic diversity (Zhang et al. 2004; Urano et al. 2013; Li et al. 2014; Wajiki et al. 2015). The crested ibis experienced a severe bottleneck and all the birds of this species in the world at present are descendants of the two breeding couples that were found in Yangxian County (Shaanxi Province, China) in 1981 (Ding 2010). The limited number of founders and small population size thus account for the present state of high inbreeding coefficient and low genetic diversity of this species, which render it more susceptible to inbreeding depression. Even though the size of the captive crested ibis population of Deqing has increased gradually to reach 209 individuals, it is still far lower than the number suggested by Franklin (1980) that is necessary to maintain sufficient genetic variation within the population over a long period (Franklin 1980). Small populations are especially sensitive to the effects of deleterious mutations and inbreeding depression (Fischer et al. 2003). Our study identified that the number of SNPs and InDels in the DEN samples were significantly lower than in the AL samples (P < 0.05), and the sequence diversity of DEN was also lower than that of AL. In addition, the number and total length of ROHs, along with Froh values in DEN samples were significantly higher than those in AL. These results all suggest that in comparison with AL, the DEN samples contained less genetic diversity and suffered higher degrees of inbreeding. Therefore, inbreeding depression may be the reason for the death of crested ibis embryos and nestlings.
Even though the effect of inbreeding depression appears to vary among species, meta-analyses of 119 zoo populations indicated that inbreeding depression has a significant impact on neonatal survival (Boakes et al. 2007). With inbreeding, the neonatal mortality of the Speke’s gazelle (Gazella spekei) increased (Boakes et al. 2007). Inbreeding populations of the domestic Pacific blue shrimp Penaeus (Litopenaeus) stylirostris also exhibit higher mortality than interbreeding populations (Goyard et al. 2008). In addition, inbreeding has adverse effects on fitness-related traits such as body mass and horn growth of the Alpine ibex (Capra ibex) (Brambilla et al. 2014). Moreover, inbred individuals may have genetic diseases, such as the chondrodystrophy observed in the California condor (Gymnogyps californianus), and may be more susceptible to infectious diseases, as avian malaria for the crow (Corvus brachyrhynchos) (Ralls and Ballou 2004; Townsend et al. 2018). These results all suggest that it is important to minimize inbreeding and protect genetic diversity within populations for preserving the population’s fitness, especially in small populations with limited founder individuals.
Congenital Factors That Cause the Death of Embryos and Nestling
Among the significantly enriched GO terms and KEGG pathways in part A of the representative regions with difference, five significantly enriched GO terms were related to protein catabolic process. This process is extremely important for embryonic survival and development as the incubated eggs lack carbohydrate, and most of the glucose that sustains embryo vigor and survival is derived from protein metabolism by glycogenesis and gluconeogenesis (Christensen et al. 2001; De Oliveira et al. 2008). The two genes UBE4A and LONP1 play particular roles in the protein catabolic process. UBE4A participates in diverse biochemical processes that are necessary for growth and differentiation (Contino et al. 2004). LONP1 can degrade proteins that are damaged by oxidation for metabolic purposes (Bezawork-Geleta et al. 2015). In addition, ADAM11 and JAM2, which belong to the GO term of binding, may act as tumor suppressors of breast cancer and colorectal cancer, respectively (Katagiri et al. 1995; Zhao et al. 2016).
The significantly enriched pathways are associated with lipid metabolism, glycan biosynthesis and metabolism, and the immune system in part A of the representative regions with difference. In the midincubation stage until 2 or 3 days prior to piping, lipid metabolism comprises the most active metabolism in embryos. In particular, the biosynthesis of fatty acids and β-oxidation are very active during this embryonic developmental stage (De Oliveira et al. 2008). Specifically, ACOT1 which pertains to the pathway of lipid metabolism can sustain free fatty acid homeostasis (Hunt et al. 2006).
During the last stage of incubation, the majority of threats that bird embryos encounter are related to energy metabolism, and the processes of glycogen synthesis and degradation are crucial for their survival (De Oliveira et al. 2008). Embryos metabolize glycogen to maintain the activity of muscle and generate heat. If the stored glycogen becomes depleted, this will affect embryo vigor and survival (Christensen et al. 1999; De Oliveira et al. 2008). HYAL1 and HYAL4 can catabolize hyaluronic acid, a kind of glycosaminoglycan (Laurent and Fraser 1992; Gakunga et al. 1997; Csoka et al. 1999).
In addition, in the last-incubation stage, embryos are immunocompetent and can generate innate and adaptive immune responses to causative agents (Ribatti 2008; Schilling et al. 2018). JAM2 belongs to the immunoglobulin superfamily and plays a crucial role in the migration of natural killer cells and cytotoxic T cells (Johnson-Léger et al. 2002; Ueki et al. 2008). These genes, GO terms, and KEGG pathways are likely to promote normal embryo development and healthy nestling growth. Conversely, the aberrant conditions of these genes, GO terms, and pathways in DEN samples may contribute to the eventual death of embryos and nestlings.
Among the genes associated with the significantly enriched GO terms in part B of the representative regions with difference, four are related to detrimental diseases. Mutations in KLHL3 can cause pseudohypoaldosteronism II (OMIM: 145260) (Louis-Dit-Picard et al. 2012; Wu and Peng 2013). SETDB2 plays oncogenic roles in gastric cancer (Nishikawaji et al. 2016). TNNT2 encodes a part of the troponin complex; this protein has been reported as an effective prognosticator for all-cause mortality of patients with acute coronary syndrome (Arnadottir et al. 2017). Mutations in PKP1 can cause ectodermal dysplasia-skin fragility syndrome, an autosomal recessive genodermatosis (Hernandez-Martin et al. 2013).
The thiamine metabolism pathway represents one of the significantly enriched KEGG pathways in part B of the representative regions with difference. Thiamine is a vital compound for multicellular organisms that participates in several metabolic processes including immune system activation, energy metabolism, and interorganellar signal transduction (Manzetti et al. 2014). Families with “near-miss” sudden infant death syndrome (SIDS) exhibit higher mortality from thiamine deficiency, with thiamine deficient infants being more likely to die because of SIDS. Such “near-miss” SIDS families might carry genetic defects of thiamine metabolism (Jeffrey et al. 1985). In addition, rats deficient in thiamine exhibit neuropathological damage of brainstem structures (Butterworth et al. 1985). Moreover, researchers have shown that as a result of thiamine deficiency, hatchlings of chickens may die from polyneuritis (Polin et al. 1962; Charles et al. 1972). AK1 pertains to the pathway of thiamine metabolism, and its mutation can cause chronic hemolytic anemia (Bianchi 1999). In addition, EXOSC3 is included in the significantly enriched pathway of RNA degradation. Mutation in this gene can cause pontocerebellar hypoplasia and spastic paraplegia, which is a serious autosomal recessive neurodegenerative disease with prenatal onset (Wan et al. 2012; Eggens et al. 2014; Halevy et al. 2014). This disease usually causes patient death within 1 year after birth (Schwabova et al. 2013). Thus, these genes, GO terms, and KEGG pathways may comprise causal factors in embryos and nestlings’ death.
Human genetic disorders like the Hirschsprung disease, spontaneous coronary artery dissection, late-onset Alzheimer’s disease, and neuromyelitis optica have been studied through whole genome resequencing (Estrada et al. 2018; Fahey et al. 2018; Tang et al. 2018; Vardarajan et al. 2018). The animal hereditary disease ascites syndrome for the chicken was also studied through this technology (Dey et al. 2018). Whole genome resequencing technology might be utilized to study other threatened bird species with high mortality of embryos and nestlings, such as the Congo peafowl (Afropavo congensis) and Guam Micronesian kingfisher (Halcyon cinnamomina cinnamomina), as the mortality of nestlings aged 0–30 days of these two birds exceeds 40% (Boakes et al. 2007). Knowing the congenital factors that cause the nestling death of these species can help acquire more knowledge of detrimental diseases and better protect these species.
In summary, high-throughput sequencing technology was utilized to resequence the whole crested ibis genome of AL and DEN individuals in this study. The DEN group exhibited a genetic background that was probably enriched for deleterious genes (KLHL3, SETDB2, TNNT2, PKP1, AK1, and EXOSC3) and it also showed higher extent of inbreeding, ultimately resulting in the death of embryos and nestlings. Moreover, in the genomic regions that differ between the alive and dead samples, we also identified several genes involved in the protein catabolic process (UBE4A and LONP1), lipid metabolism (ACOT1), glycan biosynthesis and metabolism (HYAL1 and HYAL4), and the immune system (JAM2) that are likely to promote the normal development of embryos and nestlings. The aberrant conditions of these genes and biological processes may also contribute to the eventual death of embryos and nestlings. Even though additional experiments are required to verify our results, our study provides analyses at the whole genome level of the congenital factors underlying the death of crested ibis embryos and nestlings. These data establish a foundation for more effectively managing the captive breeding populations of the crested ibis, with the goal toward ensuring the survival of this endangered species.
Supplementary Material
Supplementary data are available at Genome Biology and Evolution online.
Data deposition: Sequencing data of the crested ibis genomes of this study have been deposited at NCBI BioProject (http://www.ncbi.nlm.nih.gov/bioproject) under the accession number PRJNA508426.
Acknowledgments
This work was supported by the National Key Program (2016YFC0503200) from the Ministry of Science and Technology of China, a special grant from the State Forestry Administration, and the Fundamental Research Funds for the Central Universities of China.
