Report The Genomic Footprints of the Fall and Recovery of the Crested Ibis Highlights Authors d Recent population decline of crested ibis in 10 kya was likely Shaohong Feng, Qi Fang, shaped by human activity Ross Barnett, ..., Tomas Marques-Bonet, d Modern group has lost almost half of genetic variations in pre- M. Thomas P. Gilbert, Guojie Zhang bottleneck ancestors d Modern group suffers from a high inbreeding coefficient and Correspondence deleterious mutation load guojie.zhang@bio.ku.dk d Ancestral balancing selection is outweighed by genetic drift In Brief in the modern group Feng et al. use whole-genome sequencing of contemporary and historic crested ibis, an iconic endangered bird species, to explore how their genetic diversity has changed through time. Their analyses reveal the roles of genetic drift and intensive inbreeding on the loss of genetic diversity in today’s population. Feng et al., 2019, Current Biology 29, 340–349 January 21, 2019 ª 2018 The Author(s). Published by Elsevier Ltd. https://doi.org/10.1016/j.cub.2018.12.008 Current Biology Report The Genomic Footprints of the Fall and Recovery of the Crested Ibis Shaohong Feng,1,2,3 Qi Fang,3,4 Ross Barnett,5 Cai Li,3,6 Sojung Han,7 Martin Kuhlwilm,7 Long Zhou,3 Hailin Pan,3,4 Yuan Deng,3 Guangji Chen,1,3 Anita Gamauf,8 Friederike Woog,9 Robert Prys-Jones,10 Tomas Marques-Bonet,7,11,12,13 M. Thomas P. Gilbert,5,14 and Guojie Zhang2,3,4,15,16,* 1University of Chinese Academy of Sciences, Beijing 100049, China 2State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of Sciences, Kunming 650223, China 3China National GeneBank, BGI–Shenzhen, Shenzhen 518083, China 4Section for Ecology and Evolution, Department of Biology, University of Copenhagen, DK-2100 Copenhagen, Denmark 5Centre for GeoGenetics, Natural History Museum of Denmark, University of Copenhagen, Øster Voldgade 5-7, 1350 Copenhagen, Denmark 6The Francis Crick Institute, London NW1 1AT, UK 7Institute of Evolutionary Biology (UPF-CSIC), PRBB, Dr. Aiguader 88, 08003 Barcelona, Spain 8Museum of Natural History Vienna, 1st Zoological Department - Ornithology, Burgring 7, A-1010 Vienna, Austria 9Staatliches Museum für Naturkunde, Rosenstein 1, 70191 Stuttgart, Germany 10Bird Group, Department of Life Sciences, Natural History Museum, Akeman St, Tring, Herts HP23 6AP, UK 11Catalan Institution of Research and Advanced Studies (ICREA), Passeig de Lluı́s Companys, 23, 08010, Barcelona, Spain 12CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Baldiri i Reixac 4, 08028 Barcelona, Spain 13Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, c/ Columnes s/n, 08193 Cerdanyola del Vallès, Barcelona, Spain 14Norwegian University of Science and Technology, University Museum, 7491 Trondheim, Norway 15Center for Excellence in Animal Evolution and Genetics, Chinese Academy of Sciences, Kunming 650223, China 16Lead Contact *Correspondence: guojie.zhang@bio.ku.dk https://doi.org/10.1016/j.cub.2018.12.008 SUMMARY tion. Though its once tiny population has today recovered to >2,000 individuals [9], this process Human-induced environmental change and habitat was accompanied by almost half of ancestral loss fragmentation pose major threats to biodiversity of genetic variation and high deleterious mutation and require active conservation efforts to mitigate load. We furthermore show how genetic drift their consequences. Genetic rescue through coupled to inbreeding following the population translocation and the introduction of variation into bottleneck has largely purged the ancient polymor- imperiled populations has been argued as a phisms from the current population. In conclusion, powerful means to preserve, or even increase, the we demonstrate the unique promise of exploiting genetic diversity and evolutionary potential of en- genomic information held within museum samples dangered species [1–4]. However, factors such as for conservation and ecological research. outbreeding depression [5, 6] and a reduction in available genetic diversity render the success of RESULTS such approaches uncertain. An improved evaluation of the consequence of genetic restoration requires Historical and Contemporary Population Structure of knowledge of temporal changes to genetic diversity the Crested Ibis before and after the advent of management pro- We collected 57 historic crested ibis samples dating to be- grams. To provide such information, a growing tween 1841 and 1922. The geographic origin of these samples number of studies have included small numbers of broadly covers the historical distribution of the species, genomic loci extracted from historic and even including east China, northwest China, northeast China, Ko- rea, Japan, and Russia (Figure 1A). We initially generated an ancient specimens [7, 8]. We extend this approach average of 179 million (M) paired-end short sequencing to its natural conclusion, by characterizing the com- reads for historic samples (49 base pairs [bp] for most sam- plete genomic sequences of modern and historic ples and 100 bp for two; Data S1). After filtering to remove population samples of the crested ibis (Nipponia contaminants and other sequencing artifacts, we obtained nippon), an endangered bird that is perhaps the final genomic coverage at ca. 7-fold for each sample. Analysis most successful example of how conservation effort with mapDamage2.0 revealed signs of cytosine deamination has brought a species back from the brink of extinc- at the ends of the sequences, which are normally expected 340 Current Biology 29, 340–349, January 21, 2019 ª 2018 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). A Figure 1. Sampling Map and Population Russia Structure of the Crested Ibis (A) Map of sampling locations. Mongolia (B) NJ phylogeny of all 65 individuals (57 historic and 8 contemporary) correlates with geographic distri- Northeast Asia bution in (A), except for the one (marked in red) Northwest China North Korea collected from the NW, which is closely related to Gansu the contemporary group. Three NE subgroups South Korea Japan (NE-1, NE-2, and NE-3), EA, NW, and contemporary Shannxi groups (MC) are labeled by different background Anhui colors. Shanghai Chekiang (C) Principal component analysis for all 65 in- East China dividuals. The observed result is consistent with that of the NJ phylogeny. (D) Population structure of all 65 individuals (K = 2, 3, Northwest China and 4). The population origin of each individual is East China Northeast Asia indicated on x axis. Each individual is represented Contemporary China by a bar that is segmented into colors based on the ancestry proportions given the assumption of K B C populations. Contemporary Group 0.1 NE SZAIPI003460 61 See also Data S1 and S2. PI0034 EA CHKP 59 SZ SZAIPI0034 62 SZAI PI0034 EI14 NE-2 58 NW 03465 3 34 0802 CHK KPEI1 64 346 NE-1 34 I00 SZAI I00 I00 AIP 65 SZAIPI0 CHKS CH AIP 1 IP 6022 MC 0 SZ CH KP SZA CHKSEI14060217 022 CHKPEI14080262 CHKPEI1 CHKP CHKPE CH KS CHKP EI EI1 I140 EI14 406 CH SEI1 4080 EI1 1406 0802 408 KP 0.0 EI1 4 69 EI1 KSE 06 4060 281 EI1408 0 25 408 I14080282 0213 02 KS 4080258 025 14 CH 9 CH 1 219 0260 CHKPE I14080250 NW further divided into three subgroups 201 12 CHKPEI14080268 02 CH HKS 4060 EI1 406 4 19 4 EI1 22 KS EI1 60 CHKSEI14060215 60 40 PC2 KS 40 99 (NE-1, NE-2, and NE-3) according to their EI1 601 0200 CH C KS I140 06 −0.1 CH KSE SEI14 8 0196 CH CHK 406019 EI1406 EI1 CHKS S CHKPEI14080253 genetic distances. The remaining sam- CHK CHKSEI1 406022 2 60223 CHKPEI14080256 CHKSEI140 CHKPEI1 408025 5 ples fall into two groups: east China (EA) CHKPEI14080279 CHKSEI14060209 202 63 0802 8 408027 NE-3 CHKSEI14060195 EI14 CHKPEI1 4060 CHKP 2 −0.2 8025 18 0197 PEI1 40 2 and northwest China (NW). All contempo- EI1406 280 SEI1 0 CHKS 406 4080 206 CHK EI1 71 PEI1 60 CHK KS CHK EI140 204 CH 0802 EI14 3 KS 60 CHKP 08027 275 CH 40 5 60 93 rary samples (MC) are from northwest EI1 20 8 40 0 1 4 0 KS 060 08 20 EI1 EI1 06 KP I14 6 CH I14 KS I14 CH PE 27 SE K CH KSE 80 408 74 CH K 64 406 7 6 40 EI1 802 CH 027 021 CH CHKPEI14080266 0 I1 0802 0267 China and form a monophyletic group. 27 E 40 KP EI1 EI14080 H EI14 C CHKPEI1408 EI1 KP −0.3 KP CH KS CHKP CH CH CHKP 0.0050 EA These clustering results were also sup- −0.2 −0.1 0.0 0.1 PC1 ported by principal component analysis D 1.0 0.8 (PCA) [13] (Figure 1C). K=2 0.6 We ran STRUCTURE v2.3.4 [14] to 0.4 explore the genetic composition of each 0.2 0.0 group and subgroup after initially removing NE EA NW MC 1.0 potential bias caused by the missing loci 0.8 (Figure 1D; Data S2). Sample clusters 0.6 K=3 0.4 were evaluated using the ad hoc statistic 0.2 (DK) [15]. The contemporary group sepa- 0.0 rated from the historical populations NE EA NW MC 1.0 when setting the clusters K as 2. The DK 0.8 0.6 value reaches its maximum at K = 3, indi- K=4 0.4 cating the uppermost level of structure. 0.2 At K = 4, the clusters reflect the geographic 0.0 NE EA NW MC distribution of historical samples. Overall the population structuring indicates that while the contemporary group is geneti- from the historic DNA sequences (Data S1). We combined cally distinct to all historical groups, it shares a moderate level these data with previously published sequencing data for of genetic variation with the EA and NW groups. eight contemporary individuals [10] and used ANGSD v0.615 [11] to identify 5,268,206 SNPs in total. On average, each his- Population Size Dynamics of the Crested Ibis torical sample contained 1,088,820 (±112,638) heterozygous We used the historic genomic data to confirm the decline of the sites, while each contemporary sample contained 788,162 effective population size (Ne) reported in a previous study [10] (±35,747) heterozygous sites. (Figure S1). Because the power of PSMC (pairwise sequentially After excluding linked SNP loci that could potentially bias Markovian coalscent) to reconstruct demography is poor over clustering results (see STAR Methods), we built a neighbor- recent history due to the limited number of recombination events joining (NJ) tree [12] using all 65 samples (Figure 1B). The NJ [16], we used an alternate method based on approximate tree assigned the historical samples to three major groups, Bayesian computation, PopSizeABC v2.1 [17], to better charac- consistent with their original geographic provenance. A north- terize the Ne dynamics over the last 10,000 years, during which east Asian group (NE) exhibits the most diversity and can be the climate has been relatively stable [18]. We found that all Current Biology 29, 340–349, January 21, 2019 341 A 10 10 10 MC NW EA 1 1 1 Effective population size (x10 ) 0.1 0.1 0.1 4 0.01 0.01 0.01 0.001 0.001 0.001 0 10 100 1,000 10,000 100 110 200 1,100 10,100 100 110 200 1,100 10,100 10 10 10 NE-1 NE-2 NE-3 1 1 1 0.1 0.1 0.1 0.01 0.01 0.01 0.001 0.001 0.001 100 110 200 1,100 10,100 100 110 200 1,100 10,100 100 110 200 1,100 10,100 Years before present B Contemporary period: 19th-20th C LGM: ca. 22 kya D LIG: ca. 120 - 140kya Figure 2. Inference of Recent Population Size Tendency and Prediction of Breeding Ranges for the Crested Ibis at the Contemporary, the Last Glacial Maximum, and the Last Interglacial Periods (A) The recent effective population size (Ne) for each group is inferred using PopSizeABC. A 90% confidence interval is indicated by dotted lines. (B–D) Reconstruction of ecological niche models for the contemporary (19th and 20th centuries, B), last glacial maximum (LGM, approximately 22 kya, C) and the last interglacial (LIG, approximately 120–140 kya, D) periods. Pictures show the point-wise mean of 10 replicates using 19 environmental layers. The colors indicate the predicted probability that conditions are suitable, with red indicating high and blue indicating low probabilities. The white square points in (B) indicate the geographic coordinates of the crested ibis samples. See also Figure S1 and Table S1. groups declined at different rates over the last 10,000 years, and breeding range of the crested ibis covered broad areas across such decline was followed by a short period of relative popula- south China, with extensive hotspots of habitat (Figure 2D). Dur- tion stability before the near extinction of the crested ibis in the ing the LGM, their potential breeding range showed extensive 20th century (Figure 2A). The recent demographic pattern implies contraction, corresponding to the climate change during this that the impact of human activity on the population decline may period [20] (Figure 2C). Today, our simulations predict that the have already commenced up to 600 years ago (the starting point suitable breeding range would largely resemble that for the of decline in the NW group). LGM, with multiple distribution hotspots with suitable probability Environmental suitability is an important factor that influences higher than 50% (Figure 2B). Nevertheless, our data shows that population fluctuation. Therefore, we used ecological niche Ne continued to decline in the contemporary period (Table S1), models (ENMs) [19] to reconstruct the suitable niche for crested consistent with the fact that the historical groups went extinct ibis during the contemporary period (19th and 20th centuries), the across most of the breeding range during the 19th and 20th cen- last glacial maximum (LGM, approximately 22 thousand years turies [9]. The reconstructed suitable niche pattern implies that ago [kya]) and the last interglacial period (LIG, approximately the population collapse during the late-19th and early-20th cen- 120–140 kya) based on 35 occurrences of crested ibis samples turies was more likely a result of human activity rather than associated with a latitude and longitude record during the 19th climatic factors. This is consistent with previous analyses on and 20th centuries in the VertNet database (Figures 2B–2D) the availability of habitat suitable for crested ibis, which demon- (see STAR Methods). During the LIG period, the potential strated how radical land-cover changes (e.g., through urban 342 Current Biology 29, 340–349, January 21, 2019 development, increased cultivation of land, and increased crop tional excess of homozygous deleterious-derived mutations diversification) between the 1970s and 2000s reduced the could be likely due to the varying efficacy of purifying selection, swamp environments that used to be their feeding areas [21, 22]. resulting from differences in the demographic histories of the populations [27, 28]. We found that the ratio of homozygous sites Severe Genetic Diversity Loss in the Contemporary to homo- and heterozygous sites in the contemporary population Population Caused by Inbreeding is significantly higher than that in the historic populations (0.328 ± Although the size of the contemporary population has been 0.041 in contemporary group, 0.207 ± 0.039 in historical group, steadily increasing owing to recent conservation efforts, its Welch two-sample t test, 2.524e 05) (Figure S2E). A similar dif- nucleotide diversity (p) is considerably lower than that of the his- ference could also be seen when using loss of function (LOF) as toric population (mean value is 0.112, while mean value of the an index (Figure S2F; see STAR Methods). Despite only 100 historic samples is 0.208; Figure 3A). We also found accelerated years difference between the collection dates of historic and fixation of polymorphisms in the contemporary population, with contemporary samples, the number of homozygous deleterious ca. 70% of the SNPs that were polymorphic in the historical mutations has doubled in the contemporary population, show- groups having been fixed in the contemporary group (see casing the severity of the bottleneck. STAR Methods). We further identified the identical by descent (IBD) regions, which means the alleles derive from the ancestor Role of Natural Selection in Shaping the Genetic and maintain in descendent populations [23]. We obtained a total Diversity of the Crested Ibis of 444.5 megabases (Mb) merged excess IBD regions in the The effect of selection is considered to be limited in contrast to contemporary group (Figure 3B), 31.8 Mb for the NW, 31.0 Mb that of genetic drift in populations with a small Ne [29]. Immune for the NE, and 0 Mb for the EA groups. genes are believed to be under pathogen-mediated balancing Such dramatic loss of genetic diversity could be caused by the selection, thus maintaining high polymorphism levels in wild unavoidable effects of inbreeding. To assess the extent of populations [30, 31]. However, we found that the genetic diver- inbreeding, we estimated FUNI based on the correlation between sity of the major histocompatibility complex (MHC) has experi- uniting gametes following Wright [24]. A lower (i.e., negative) enced drastic haplotype loss in the contemporary group. The value indicates an excess of heterozygosity and little inbreeding overall allelic richness (AR) of the entire MHC region in the effect, while a positive value indicates an excess of homozygos- contemporary group is significantly lower than that shown in ity and strong inbreeding effect. We subsequently calculated the historical groups (meanAR of the contemporary group is FUNI and Ne for 14 animal species spanning different conserva- 2.44, meanAR of the historical group is 2.92, Welch two-sample tion statuses. A regression correlation including the crested t test, p < 2.2e 16; see STAR Methods). Furthermore, haplotype ibis and other species clearly shows that the contemporary pop- blocks in the contemporary group were much longer and ex- ulation is a significant outlier, highlighting the severity of its hibited lower diversity (Figures 4 and S3). This dramatically inbreeding (Figure 3C; Table S2; see STAR Methods). This reduced diversity of the MHC in the contemporary population finding is consistent with the fact that the gene pool of all may confer a high risk of increasing the burden of transmissible contemporary individuals was inherited from only two breeding pathogens when captive animals are released into the wild [32]. pairs. We subsequently performed a broader scan across the coding region to identify genomic loci that were previously under Patterns of Deleterious Mutations across Populations balancing selection but have been fixed in the current population Support Ne Collapse, Bottlenecks, and Inbreeding in the (see STAR Methods). In total we found 53 genes with high Taji- Contemporary Population ma’s D values in at least two of the historic groups, implying Recent conservation genetic studies have demonstrated that that they were once under balancing selection. Of these, nine small populations are susceptible to allelic drift, leading to allele genes exhibit a negative Tajima’s D value in the contemporary loss or fixation, which can be evaluated using the site-frequency- group that might have been substantially affected by bottleneck spectrum (SFS) [25]. In our analysis, all historic populations effects during the breeding program. Several of these genes are showed a general pattern in which most mutations occurred at involved in processes relating to immunity, the nervous system, low frequencies, while few mutations segregated at high fre- and reproduction. For example, polymorphism has been quencies (Figures S2A–S2C). However, in the contemporary completely lost in half of the SNP loci in gene KPNA1, a member population, we observed a reduced number of low-frequency of the importin or karyopherin alpha family that has been re- mutations and a relatively higher fraction of mutations that segre- ported to be involved in the development and/or repair pro- gate at medium or higher frequencies (Figure S2D). We hypoth- cesses of the nervous system [33]. esize that these mutations occur at low rates and that relaxed purifying selection has probably led to random fixation, likely DISCUSSION due to the high level of inbreeding in the contemporary popula- tion. An alternative possibility might be that the contemporary Thanks to the development of paleogenomic sequencing tech- individuals sampled were from a sub-structured population; niques and the recovery of genetic data from museum samples however, we find no such evidence in the results of the admixture collected up to 100 years before the introduction of the modern analysis to support this hypothesis. management program, our study has reconstructed a much- We further used the Grantham score [26] as an index to identify improved demographic history of the crested ibis over the last potential deleteriousness of missense mutations in the coding 10,000 years. The continued population decline of this species regions of historic and contemporary populations. The propor- in this period therefore not only leads to the extreme bottleneck Current Biology 29, 340–349, January 21, 2019 343 A 0.30 0.25 0.20 0.15 π 0.10 Historical groups 0.05 Contemporary group 0.00 100 80 60 40 20 0 1 2 3 4 5 6 7 8 9 10 11 12 13 1415 Z Frequency Chromosome B Chr 1 Chr 2 Chr 9 Chr 3 Chr 10 Chr 4 Chr 11 Chr 5 Chr 12 Chr 6 Chr 13 Chr 7 Chr 14 Chr 8 Chr 15 Chr Z C 0.2 Conservation Status Least concern (LC) M MC EN Vulnerable (VU) Inbreeding coefficient (FUNI ) Endangered (EN) Critically endangered (CR) EN N Extinct (EX) .1 0.1 LC C LC VU U CR C LC EN E EX LC EX 0.0 EX CR R CR C LC EA NW VU LC NE EX −0.1 R2 : 0.7282 P value : 6.35e−05 0 100 200 300 Effective population size (Ne) Figure 3. Loss of Genetic Diversity and Elevated Inbreeding in the Crested Ibis (A) Left: A histogram describing mean p for 5 Mb windows across the crested ibis genome in the historical and contemporary groups. Right: Genomic distribution of individual pairwise estimates of mean p in 5 Mb windows across the crested ibis genome in historical and contemporary groups. Chromosome boundaries are indicated as vertical dashed lines. (B) The chromosomal distribution of IBD regions in the contemporary population. IBD regions are marked in orange, with additional red highlighting if supported by more than half samples. (C) Linear model of inbreeding coefficient (FUNI) and effective population size (Ne) based on 14 animal species using the published population datasets. R-squared value and the p value of F test indicated a significant positive correlation between FUNI and Ne, while the contemporary population (labeled as MC) is a significant outlier. The historical crested ibis populations are indicated in gray. See also Figure S2 and Table S2. that started in the late-19th century but also provides new in- climate of the past 10,000 years should have stopped the prior sights into the negative effect of human activities. Given that decline and incurred a recovery—something we see no evidence crested ibis did not appear to favor colder climates, the warmer of. One explanation for this discrepancy might be that the 344 Current Biology 29, 340–349, January 21, 2019 A B Figure 4. Haplotype Structure and Length Distribution of the Class I Region in Historical and Contemporary Groups (A) Haplotype blocks are detected based on the SNP datasets for historical and contemporary groups, respectively. Colored circles show the haplotype structure (each box indicates the start and end SNPs of one haplotype) and the diversity along the class I region (color indicates the allelic richness). The allelic richness of (legend continued on next page) Current Biology 29, 340–349, January 21, 2019 345 crested ibis had entered into an evolutionary dead end from environment, given its lack of competitors and predators, thus which they were unable recover, because of their small Ne and to some extent reducing the negative effect of genetic load. the low genetic diversity. However, given the success of the cur- Ultimately, recent successful re-introduction of the crested rent management program in increasing their population size to ibis into the wild has provided an important step forward in this ca. 2,600 individuals from only seven founding individuals, it is regard [41]. Long-term monitoring of their genetic diversity clear that the species still maintains a robust reproductive rate. in the wild will provide a crucial index for their long-term An alternate explanation for their population decline might there- persistence and management. Recently, efforts have been fore relate to human interference. Human activities, like over- commenced to establish a DNA identification profiling (DIP) plat- hunting and development of farm land, might have led to the form using short tandem repeat markers to help avoid close fragmentation of niches for the crested ibis and ultimately crosses [42]. Given that our analyses identified several candidate increased the inbreeding coefficient within the isolated subpop- genes that no longer exhibit high polymorphism, we suggest it ulations. Consistently, recent research on crested ibis notes how may be useful to incorporate them into the existing DIP for moni- high human-induced mortality risks play a significantly negative toring and preserving genomic diversity. Lastly, our study high- effect on the population growth and range expansion of the sur- lights how integrating genomic information from the modern viving population [34]. and historic samples can help to reveal the evolutionary history One factor that contributed to the crested ibis’s recent demog- of endangered species and in doing so provide guiding informa- raphy is the effect of inbreeding. There is little doubt that tion of benefit to future conservation programs. inbreeding reduces the fitness of a population [35–37], particu- larly for critically endangered species with small population STAR+METHODS sizes. However, previous studies that have argued that inbreeding depression can elevate extinction risk are based on Detailed methods are provided in the online version of this paper simulations [38], and thus there is a need for actual evidence and include the following: with which to verify whether the consequences of inbreeding, such as the accumulation of deleterious alleles, will lead to and d KEY RESOURCES TABLE even accelerate the loss of diversity. Our findings indicate that d CONTACT FOR RESOURCE SHARING the contemporary crested ibis population has dramatically d EXPERIMENTAL MODEL AND SUBJECT DETAILS reduced the ancestral genetic diversity, undergone higher allelic d METHOD DETAILS fixation, and carried an elevated load of homozygous deleterious B Whole-genome sequencing of museum samples mutations that associated with a raised inbreeding coefficient. B Read alignment and SNP calling Importantly, the inclusion of samples spanning 1841–2011 CE d QUANTIFICATION AND STATISTICAL ANALYSIS enabled us to accurately quantify the extent to which the most B Data pre-processing recent demographic decline has led to a loss in their genetic vari- B Population structure ation as well as increased inbreeding levels and genetic load, all B Population demographic analysis of which are regarded as useful genomic parameters for endan- B Ecological niche models analysis gered species [39]. Our analyses showed how the overall nucle- B Temporal population polymorphism otide diversity of the contemporary group is only 53.85% of that B Correlation between inbreeding coefficient and effec- which was present ca. 100 years ago. This has been accompa- tive population size (Ne) nied by an almost 4-fold elevation in their inbreeding coefficient B Patterns of deleterious mutations across populations and a doubling of their homozygous deleterious mutation load. B MHC diversity analysis The conservation program has had remarkable success in B Natural selection analysis increasing the size of the crested ibis population, from the seven d DATA AND SOFTWARE AVAILABILITY individuals during their near extinction to currently ca. 2,600 indi- SUPPLEMENTAL INFORMATION viduals. Nevertheless, the contemporary population suffers from dramatic diversity loss across the genome. This overall reduced Supplemental Information includes three figures, two tables, and two data files genetic diversity represents a risk when they are introduced to and can be found with this article online at https://doi.org/10.1016/j.cub.2018. the natural environment as part of the conservation program. 12.008. We acknowledge, however that it may be, that the protected environment of the management program will help them to over- ACKNOWLEDGMENTS come the potential harmful effects from the reduced diversity, We are grateful to the National Museum of Natural History (USNM) (Gary thus enabling them to survive into the future with low genetic Graves, Chris Milensky); Senckenberg Museum of Natural History (Gerald diversity, as for example has been the case for the endangered Mayr); Natural History Museum, Tring (Hein van Grouw); Department of island fox over the past millennia [40]. In that example, it is spec- Zoology, National Museums Liverpool (Tony Parker, Clemency Fisher); Center ulated that their isolated island habitat might provide a friendly for Natural History/Zoological Museum Hamburg (Alexander Has); Naturalis historical groups is estimated in the random way (see STAR Methods). Gaps in the colored circles are the intervals between two neighboring haplotype blocks. Genomic structures of the class I region are drawn as gray bars, with different sizes showing different gene loci of varied sizes. (B) Length distribution of haplotype blocks. Haplotype blocks with lengths larger than 2,000 bp are grouped into the last bar. The plotting area for the bin of 250– 2,000 bp is enlarged in the top right corner. See also Figure S3. 346 Current Biology 29, 340–349, January 21, 2019 Biodiversity Center Leiden, the Netherlands (Pepijn Kamminga); NRM Stock- and rebirth of the crested ibis and other endangered bird species. holm (Ulf Johansson); and Natural History Museum of Denmark (Jon Fjeldså) Genome Biol. 15, 557. for providing access to museum specimens and Nigel Collar for providing 11. Korneliussen, T.S., Albrechtsen, A., and Nielsen, R. (2014). ANGSD: background information critical to locating the samples. This work was sup- Analysis of Next Generation Sequencing Data. BMC Bioinformatics ported by the Strategic Priority Research Program of Chinese Academy of 15, 356. Sciences (XDB31020000, XDB13000000), Carlsberg Foundation grant to 12. Saitou, N., and Nei, M. (1987). The neighbor-joining method: a new G.Z. (CF16-0663), ERC Consolidator Grant 681396 ‘Extinction Genomics’ method for reconstructing phylogenetic trees. Mol. Biol. Evol. 4, (M.T.P.G.), BFU2017-86471-P (MINECO/FEDER, UE) (T.M.-B.), U01 406–425. MH106874 grant (T.M.-B.), Howard Hughes International Early Career (T.M.-B.), Obra Social ‘‘La Caixa’’ and Secretaria d’Universitats i Recerca 13. Luu, K., Bazin, E., and Blum, M.G. (2017). pcadapt: an R package to (T.M.-B.), CERCA Programme del Departament d’Economia i Coneixement perform genome scans for selection based on principal component anal- de la Generalitat de Catalunya (T.M.-B.), and Deutsche Forschungsgemein- ysis. Mol. Ecol. Resour. 17, 67–77. schaft (DFG) fellowship (KU 3467/1-1) (M.K.). 14. Falush, D., Stephens, M., and Pritchard, J.K. (2003). Inference of popula- tion structure using multilocus genotype data: linked loci and correlated AUTHOR CONTRIBUTIONS allele frequencies. Genetics 164, 1567–1587. 15. Evanno, G., Regnaut, S., and Goudet, J. (2005). Detecting the number of C.L., M.T.P.G. and G.Z. conceived the study; M.T.P.G. collected museum clusters of individuals using the software STRUCTURE: a simulation study. samples, which were prepared into sequencing libraries by R.B.; C.L. Mol. Ecol. 14, 2611–2620. managed sample sequencing; C.L., Y.D., L.Z., and H.P. developed, tested, 16. Li, H., and Durbin, R. (2011). Inference of human population history from and implemented the detailed SNPs calling pipeline; S.F., Q.F., C.L., and individual whole-genome sequences. Nature 475, 493–496. G.C. carried out genetic analyses described in the manuscript; S.H., 17. Boitard, S., Rodrı́guez, W., Jay, F., Mona, S., and Austerlitz, F. (2016). M.K. and T.M.-B. assisted with the mutation load analysis; S.F., Q.F., C.L., Inferring Population Size History from Large Samples of Genome-Wide T.M.-B., M.T.P.G. and G.Z. wrote and edited the manuscript; All authors Molecular Data - An Approximate Bayesian Computation Approach. read and approved the final manuscript. PLoS Genet. 12, e1005877. 18. Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., DECLARATION OF INTERESTS and Leuenberger, M. (2014). Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core. Clim. Past 10, 887–902. The authors declare no competing interests. 19. Hung, C.M., Shaner, P.J., Zink, R.M., Liu, W.C., Chu, T.C., Huang, W.S., and Li, S.H. (2014). Drastic population fluctuations explain the rapid Received: July 31, 2018 extinction of the passenger pigeon. Proc. Natl. Acad. Sci. USA 111, Revised: October 9, 2018 10636–10641. Accepted: December 6, 2018 Published: January 10, 2019 20. Tian, Z., and Jiang, D. (2016). Revisiting last glacial maximum climate over China and East Asian monsoon using PMIP3 simulations. Palaeogeogr. Palaeoclimatol. Palaeoecol. 453, 115–126. REFERENCES 21. Mochizuki, S., and Murakami, T. (2010). Temporal change of crested ibis 1. Weeks, A.R., Sgro, C.M., Young, A.G., Frankham, R., Mitchell, N.J., Miller, habitat in Shaanxi Province, China. Proceedings of the 31st Asian K.A., Byrne, M., Coates, D.J., Eldridge, M.D., Sunnucks, P., et al. (2011). Conference on Remote Sensing. Assessing the benefits and risks of translocations in changing environ- 22. Li, X.H., Tian, H.D., and Li, D.M. (2009). Why the crested ibis declined in the ments: a genetic perspective. Evol. Appl. 4, 709–725. middle twentieth century. Biol. Conserv. 18, 2165–2172. 2. Armstrong, D.P., and John, L.C. (1995). Effects of familiarity on the 23. Albrechtsen, A., Moltke, I., and Nielsen, R. (2010). Natural selection and outcome of translocations, II. A test using New Zealand Robins. Biol. the distribution of identity-by-descent in the human genome. Genetics Conserv. 71, 281–288. 186, 295–308. 3. Whiteley, A.R., Fitzpatrick, S.W., Funk, W.C., and Tallmon, D.A. (2015). 24. Wright, S. (1922). Coefficients of inbreeding and relationship. Am. Nat. 56, Genetic rescue to the rescue. Trends Ecol. Evol. 30, 42–49. 330–338. 4. Weeks, A.R., Heinze, D., Perrin, L., Stoklosa, J., Hoffmann, A.A., van 25. Bouzat, J.L. (2010). Conservation genetics of population bottlenecks: the Rooyen, A., Kelly, T., and Mansergh, I. (2017). Genetic rescue increases role of chance, selection, and history. Conserv. Genet. 11, 463–478. fitness and aids rapid recovery of an endangered marsupial population. 26. Grantham, R. (1974). Amino acid difference formula to help explain protein Nat. Commun. 8, 1071. evolution. Science 185, 862–864. 5. Olden, J.D., Leroy Poff, N., Douglas, M.R., Douglas, M.E., and Fausch, 27. Lohmueller, K.E. (2014). The distribution of deleterious genetic variation in K.D. (2004). Ecological and evolutionary consequences of biotic homoge- human populations. Curr. Opin. Genet. Dev. 29, 139–146. nization. Trends Ecol. Evol. 19, 18–24. 28. Lohmueller, K.E., Indap, A.R., Schmidt, S., Boyko, A.R., Hernandez, R.D., 6. Leimu, R., and Fischer, M. (2010). Between-population outbreeding af- Hubisz, M.J., Sninsky, J.J., White, T.J., Sunyaev, S.R., Nielsen, R., et al. fects plant defence. PLoS ONE 5, e12614. (2008). Proportionally more deleterious genetic variation in European 7. Draheim, H.M., Baird, P., and Haig, S.M. (2012). Temporal analysis of than in African populations. Nature 451, 994–997. mtDNA variation reveals decreased genetic diversity in least terns. 29. Woolfit, M. (2009). Effective population size and the rate and pattern of Condor 114, 145–154. nucleotide substitutions. Biol. Lett. 5, 417–420. 8. Wandeler, P., Hoeck, P.E., and Keller, L.F. (2007). Back to the future: 30. Reed, D.H., and Frankham, R. (2003). Correlation between fitness and ge- museum specimens in population genetics. Trends Ecol. Evol. 22, netic diversity. Conserv. Biol. 17, 230–237. 634–642. 31. Messaoudi, I., Guevara Patiño, J.A., Dyall, R., LeMaoult, J., and Nikolich- 9. Collar, N.J., Andreev, A., Chan, S., Crosby, M., Subramanya, S., and Zugich, J. (2002). Direct link between mhc polymorphism, T cell avidity, Tobias, J. (2001). Threatened birds of Asia: the BirdLife International red and diversity in immune defense. Science 298, 1797–1800. data book, Fifth Edition (Cambridge: BirdLife International). 32. Grenfell, B.T., Dobson, A.P., and Moffatt, H.K.E. (1995). Ecology of 10. Li, S., Li, B., Cheng, C., Xiong, Z., Liu, Q., Lai, J., Carey, H.V., Zhang, Q., Infectious Diseases in Natural Populations (Cambridge University Zheng, H., Wei, S., et al. (2014). Genomic signatures of near-extinction Press). Current Biology 29, 340–349, January 21, 2019 347 33. Wang, B., Li, Z., Xu, L., Goggi, J., Yu, Y., and Zhou, J. (2004). Molecular 51. Murray, G.G.R., Soares, A.E.R., Novak, B.J., Schaefer, N.K., Cahill, J.A., cloning and characterization of rat karyopherin alpha 1 gene: structure Baker, A.J., Demboski, J.R., Doll, A., Da Fonseca, R.R., Fulton, T.L., and expression. Gene 331, 149–157. et al. (2017). Natural selection shaped the rise and fall of passenger pigeon 34. Sun, Y., Wang, T., Skidmore, A.K., Palmer, S.C.F., Ye, X., Ding, C., and genomic diversity. Science 358, 951–954. Wang, Q. (2016). Predicting and understanding spatio-temporal dynamics 52. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., and Jarvis, A. of species recovery: implications for Asian crested ibis Nipponia nippon (2005). Very high resolution interpolated climate surfaces for global land conservation in China. Divers. Distrib. 22, 893–904. areas. Int. J. Climatol. 25, 1965–1978. 35. DeRose, M.A., and Roff, D.A. (1999). A Comparison of Inbreeding 53. Schubert, M., Ermini, L., Der Sarkissian, C., Jónsson, H., Ginolhac, A., Depression in Life-History and Morphological Traits in Animals. Schaefer, R., Martin, M.D., Fernández, R., Kircher, M., McCue, M., et al. Evolution 53, 1288–1292. (2014). Characterization of ancient and modern genomes by SNP detec- 36. Keller, L.F., and Waller, D.M. (2002). Inbreeding effects in wild populations. tion and phylogenomic and metagenomic analysis using PALEOMIX. Trends Ecol. Evol. 17, 230–241. Nat. Protoc. 9, 1056–1082. 37. Charlesworth, D., and Charlesworth, B. (1987). Inbreeding depression 54. Simonsen, M., Mailund, T., and Pedersen, C.N. (2008). Rapid neighbour- and its evolutionary consequences. Annu. Rev. Ecol. Syst. 18, joining (Springer). 237–268. 55. Phillips, S.J., Anderson, R.P., and Schapire, R.E. (2006). Maximum en- 38. O’Grady, J.J., Brook, B.W., Reed, D.H., Ballou, J.D., Tonkyn, D.W., tropy modeling of species geographic distributions. Ecol. Modell. 190, and Frankham, R. (2006). Realistic levels of inbreeding depression 231–259. strongly affect extinction risk in wild populations. Biol. Conserv. 133, 56. R Core Development Team (2016). R: a language and environment for sta- 42–51. tistical computing (R Foundation for Statistical Computing). 39. Frankham, R. (2005). Genetics and extinction. Biol. Conserv. 126, 57. Do, C., Waples, R.S., Peel, D., Macbeth, G.M., Tillett, B.J., and Ovenden, 131–140. J.R. (2014). NeEstimator v2: re-implementation of software for the estima- 40. Robinson, J.A., Ortega-Del Vecchyo, D., Fan, Z., Kim, B.Y., vonHoldt, tion of contemporary effective population size (Ne ) from genetic data. Mol. B.M., Marsden, C.D., Lohmueller, K.E., and Wayne, R.K. (2016). Ecol. Resour. 14, 209–214. Genomic flatlining in the endangered island fox. Curr. Biol. 26, 1183– 58. Yang, J., Lee, S.H., Goddard, M.E., and Visscher, P.M. (2011). GCTA: a 1189. tool for genome-wide complex trait analysis. Am. J. Hum. Genet. 88, 41. Yu, X.P., Chang, X.Y., Li, X., Chen, W.G., and Shi, L. (2009). Return of the 76–82. Crested Ibis Nipponia nippon: a reintroduction programme in Shaanxi 59. Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M.A., Bender, province, China. BirdingASIA 11, 80–82. D., Maller, J., Sklar, P., de Bakker, P.I., Daly, M.J., and Sham, P.C. (2007). 42. Oldoni, F., Castella, V., and Hall, D. (2015). A novel set of DIP-STR markers PLINK: a tool set for whole-genome association and population-based for improved analysis of challenging DNA mixtures. Forensic Sci. Int. linkage analyses. Am. J. Hum. Genet. 81, 559–575. Genet. 19, 156–164. 60. Gusev, A., Lowe, J.K., Stoffel, M., Daly, M.J., Altshuler, D., Breslow, J.L., 43. Prado-Martinez, J., Sudmant, P.H., Kidd, J.M., Li, H., Kelley, J.L., Lorente- Friedman, J.M., and Pe’er, I. (2009). Whole population, genome-wide Galdos, B., Veeramah, K.R., Woerner, A.E., O’Connor, T.D., Santpere, G., mapping of hidden relatedness. Genome Res. 19, 318–326. et al. (2013). Great ape genetic diversity and population history. Nature 61. Paradis, E. (2010). pegas: an R package for population genetics with an in- 499, 471–475. tegrated-modular approach. Bioinformatics 26, 419–420. 44. Kukekova, A.V., Johnson, J.L., Xiang, X., Feng, S., Liu, S., Rando, 62. Cingolani, P., Platts, A., Wang, L., Coon, M., Nguyen, T., Wang, L., Land, H.M., Kharlamova, A.V., Herbeck, Y., Serdyukova, N.A., Xiong, Z., S.J., Lu, X., and Ruden, D.M. (2012). A program for annotating and predict- et al. (2018). Red fox genome assembly identifies genomic regions ing the effects of single nucleotide polymorphisms, SnpEff: SNPs in the associated with tame and aggressive behaviours. Nat Ecol Evol 2, genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 1479–1491. (Austin) 6, 80–92. 45. Liu, S., Lorenzen, E.D., Fumagalli, M., Li, B., Harris, K., Xiong, Z., Zhou, L., 63. Kimmel, G., and Shamir, R. (2005). GERBIL: Genotype resolution and Korneliussen, T.S., Somel, M., Babbitt, C., et al. (2014). Population geno- block identification using likelihood. Proc. Natl. Acad. Sci. USA 102, mics reveal recent speciation and rapid evolutionary adaptation in polar 158–162. bears. Cell 157, 785–794. 64. Pfeifer, B., Wittelsbürger, U., Ramos-Onsins, S.E., and Lercher, M.J. 46. Zhang, Z., Jia, Y., Almeida, P., Mank, J.E., van Tuinen, M., Wang, Q., (2014). PopGenome: an efficient Swiss army knife for population genomic Jiang, Z., Chen, Y., Zhan, K., Hou, S., et al. (2018). Whole-genome rese- analyses in R. Mol. Biol. Evol. 31, 1929–1936. quencing reveals signatures of selection and timing of duck domestica- tion. Gigascience 7, giy027. 65. Rohland, N., and Hofreiter, M. (2007). Ancient DNA extraction from bones and teeth. Nat. Protoc. 2, 1756–1762. 47. Li, D., Che, T., Chen, B., Tian, S., Zhou, X., Zhang, G., Li, M., Gaur, U., Li, Y., Luo, M., et al. (2017). Genomic data for 78 chickens from 14 popula- 66. Vilstrup, J.T., Seguin-Orlando, A., Stiller, M., Ginolhac, A., Raghavan, M., tions. Gigascience 6, 1–5. Nielsen, S.C., Weinstock, J., Froese, D., Vasiliev, S.K., Ovodov, N.D., et al. n, M.S., Maqbool, K., Grabherr, M., (2013). Mitochondrial phylogenomics of modern and ancient equids. PLoS 48. Lamichhaney, S., Berglund, J., Alme ONE 8, e55950. Martinez-Barrio, A., Promerová, M., Rubin, C.J., Wang, C., Zamani, N., et al. (2015). Evolution of Darwin’s finches and their beaks revealed by 67. Paijmans, J.L.A., Barnett, R., Gilbert, M.T.P., Zepeda-Mendoza, M.L., genome sequencing. Nature 518, 371–375. Reumer, J.W.F., de Vos, J., Zazula, G., Nagel, D., Baryshnikov, G.F., Leonard, J.A., et al. (2017). Evolutionary History of Saber-Toothed Cats 49. Nadachowska-Brzyska, K., Burri, R., Smeds, L., and Ellegren, H. (2016). Based on Ancient Mitogenomics. Curr. Biol. 27, 3330–3336.e5. PSMC analysis of effective population sizes in molecular ecology and its application to black-and-white Ficedula flycatchers. Mol. Ecol. 25, ale, D., Coltman, D.W., and Pontier, D. (2007). Detecting 68. Kaeuffer, R., Re 1058–1072. population structure using STRUCTURE software: effect of background 50. Zhao, S., Zheng, P., Dong, S., Zhan, X., Wu, Q., Guo, X., Hu, Y., He, W., linkage disequilibrium. Heredity (Edinb) 99, 374–380. Zhang, S., Fan, W., et al. (2013). Whole-genome sequencing of giant 69. Nadachowska-Brzyska, K., Li, C., Smeds, L., Zhang, G., and Ellegren, H. pandas provides insights into demographic history and local adaptation. (2015). Temporal Dynamics of Avian Populations during Pleistocene Nat. Genet. 45, 67–71. Revealed by Whole-Genome Sequences. Curr. Biol. 25, 1375–1380. 348 Current Biology 29, 340–349, January 21, 2019 70. Li, W.H., Wu, C.I., and Luo, C.C. (1984). Nonrandomness of point mutation regions of extensive and minimal linkage disequilibrium in human Xq25 as reflected in nucleotide substitutions in pseudogenes and its evolu- and Xq28. Nat. Genet. 25, 324–328. tionary implications. J. Mol. Evol. 21, 58–71. 73. Daly, M.J., Rioux, J.D., Schaffner, S.F., Hudson, T.J., and Lander, E.S. 71. Chen, L.C., Lan, H., Sun, L., Deng, Y.L., Tang, K.Y., and Wan, Q.H. (2015). (2001). High-resolution haplotype structure in the human genome. Nat. Genomic organization of the crested ibis MHC provides new insight into Genet. 29, 229–232. ancestral avian MHC structure. Sci. Rep. 5, 7963. 72. Taillon-Miller, P., Bauer-Sardiña, I., Saccone, N.L., Putzel, J., Laitinen, T., 74. Tajima, F. (1989). Statistical method for testing the neutral mutation hy- Cao, A., Kere, J., Pilia, G., Rice, J.P., and Kwok, P.Y. (2000). Juxtaposed pothesis by DNA polymorphism. Genetics 123, 585–595. Current Biology 29, 340–349, January 21, 2019 349 STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Biological Samples 27 Nipponia nippon tissue samples American Museum of Data S1 Natural History 19 Nipponia nippon tissue samples Natural History Museum Data S1 (London) 2 Nipponia nippon tissue samples World Museum, Liverpool Data S1 1 Nipponia nippon tissue sample Naturalis Biodiversity Center Data S1 Leiden, the Netherlands 1 Nipponia nippon tissue sample Senckenberg Natural History Data S1 Museum 2 Nipponia nippon tissue samples Stuttgart State Museum of Data S1 Natural History 2 Nipponia nippon tissue samples Smithsonian Institute National Data S1 Museum of Natural History 2 Nipponia nippon tissue samples Natural History Museum Data S1 Vienna 1 Nipponia nippon tissue sample Natural History Museum of Data S1 Denmark Deposited Data Whole-genome sequencing data of 57 historical This paper NCBI Project ID: PRJNA481495; Nipponia nippon samples CNSA: CNP0000084 Whole-genome sequencing data of 8 contemporary [10] NCBI Project ID: PRJNA232572 Nipponia nippon samples Nipponia nippon reference genome [10] NCBI Project ID: PRJNA232572 Whole-genome sequencing data of 13 Pan paniscus [43] SRA: SRP018689 Whole-genome sequencing data of 5 Pongo abelii [43] SRA: SRP018689 Whole-genome sequencing data of 5 Pongo pygmaeus [43] SRA: SRP018689 Whole-genome sequencing data of 25 Pan troglodytes [43] SRA: SRP018689 Whole-genome sequencing data of 31 Gorilla gorilla [43] SRA: SRP018689 Whole-genome sequencing data of 10 Vulpes vulpes [44] SRA: SRP100625 Whole-genome sequencing data of 18 Ursus maritimus [45] SRA: SRA092289 Whole-genome sequencing data of 10 Ursus arctos [45] SRA: SRA092289 Whole-genome sequencing data of 14 Anas [46] SRA: SRP125660 platyrhynchos Whole-genome sequencing data of 6 Tibetan fowl [47] SRA: SRP040477 from Qinghai Whole-genome sequencing data of 5 Platyspiza [48] SRA: SRP048643 crassirostris Whole-genome sequencing data of 20 Ficedula [49] SRA: ERA362744, albicollis ERA362838, ERA391986, ERA391991 Whole-genome sequencing data of 34 Ailuropoda [50] SRA: SRP013618 melanoleuca Whole-genome sequencing data of 6 Ectopistes [51] SRA: SRP042357, migratorius SRP103042 MHC class 1 region sequences of Nipponia nippon NCBI KP182409.1 MHC class 2 region sequences of Nipponia nippon NCBI KP182408.1 (Continued on next page) e1 Current Biology 29, 340–349.e1–e7, January 21, 2019 Continued REAGENT or RESOURCE SOURCE IDENTIFIER MHC extended region sequences of Nipponia nippon NCBI KP182407.1 Locality data of crested ibis fromm American Museum VertNet database http://ipt.vertnet.org:8080/ipt/resource.do? of Natural History. AMNH Bird Collection (published on r=amnh_birds 2013-10-30) Locality data of crested ibis from California Academy VertNet database http://ipt.calacademy.org:8080/ipt/resource. of Sciences. CAS Ornithology (ORN) (published on do?r=orn 2016-10-31) Locality data of crested ibis from Museum of Comparative VertNet database http://digir.mcz.harvard.edu/ipt/resource.do? Zoology, Harvard University (published on 2017-07-11) r=mcz_subset_for_vertnet Locality data of crested ibis from Museum of Vertebrate VertNet database http://ipt.vertnet.org:8080/ipt/resource.do? Zoology, UC Berkeley. MVZ Bird Collection (Arctos) r=mvz_bird (published on 2015-10-27) Locality data of crested ibis from National Museum of VertNet database https://collections.nmnh.si.edu/ipt/resource? Natural History, Smithsonian Institution. NMNH Birds r=nmnh_extant_dwc-a (published on 2016-07-21) Locality data of crested ibis (1842.1.19.90) from Natural VertNet database http://data.nhm.ac.uk/object/7e36ff64- History Museum. Natural History Museum (London) 1cf0-44f6-b30d-ae0715a85a53 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1852.2.5.12) from Natural VertNet database http://data.nhm.ac.uk/object/5a4a8a93- History Museum. Natural History Museum (London) ead3-42ce-9229-a8d34a8d0a0e Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1891.10.19.20) from Natural VertNet database http://data.nhm.ac.uk/object/b1a0a925- History Museum. Natural History Museum (London) 5132-40c7-a9c3-2a999d336660 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1892.4.2.492) from Natural VertNet database http://data.nhm.ac.uk/object/1a589acb- History Museum. Natural History Museum (London) 26c9-48c0-9365-6a7620d12e5d Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1892.4.2.493) from Natural VertNet database http://data.nhm.ac.uk/object/7be93e6c- History Museum. Natural History Museum (London) f56d-4922-95ce-31635a6e3c30 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1897.10.30.2) from Natural VertNet database http://data.nhm.ac.uk/object/f7ca4de8- History Museum. Natural History Museum (London) 0772-4743-9c82-59656231504d Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1897.10.30.3) from Natural VertNet database http://data.nhm.ac.uk/object/f9b28058- History Museum. Natural History Museum (London) f974-40c1-9300-0248271cf128 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1897.10.30.4) from Natural VertNet database http://data.nhm.ac.uk/object/927fb79d- History Museum. Natural History Museum (London) 4100-49a0-894a-3a8cfa644e1f Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1897.10.30.5) from Natural VertNet database http://data.nhm.ac.uk/object/218dc4cd- History Museum. Natural History Museum (London) 6aa4-4349-956b-9bd253b1e036 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1897.10.30.6) from Natural VertNet database http://data.nhm.ac.uk/object/69828014- History Museum. Natural History Museum (London) 8f2c-4e67-8fc2-f1ac1e623b63 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1900.9.9.12) from Natural VertNet database http://data.nhm.ac.uk/object/e04fe7fe- History Museum. Natural History Museum (London) 46c3-43e2-89d4-cb106558181f Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1900.9.9.13) from Natural VertNet database http://data.nhm.ac.uk/object/1e5fed44- History Museum. Natural History Museum (London) 06b0-4cfe-959c-e1d3902947f3 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1908.1.5.20) from Natural VertNet database http://data.nhm.ac.uk/object/3ff29518- History Museum. Natural History Museum (London) a451-416a-9e56-42992ea304fb Collection Specimens (published on 2015-10-10) (Continued on next page) Current Biology 29, 340–349.e1–e7, January 21, 2019 e2 Continued REAGENT or RESOURCE SOURCE IDENTIFIER Locality data of crested ibis (1908.1.5.21) from Natural VertNet database http://data.nhm.ac.uk/object/1a3f8ba6- History Museum. Natural History Museum (London) 4d7f-45fe-b165-07fd1ef26207 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1908.1.5.22) from Natural VertNet database http://data.nhm.ac.uk/object/786144ad- History Museum. Natural History Museum (London) 774e-49d8-9ed9-b6378db8cdfe Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1908.1.5.23) from Natural VertNet database http://data.nhm.ac.uk/object/7a9df17a- History Museum. Natural History Museum (London) c7d3-4d73-8ddf-e72c3ece16e6 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1908.1.5.24) from Natural VertNet database http://data.nhm.ac.uk/object/365367ee- History Museum. Natural History Museum (London) fd32-4981-a62d-853df63b0a35 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1908.1.5.25) from Natural VertNet database http://data.nhm.ac.uk/object/35a1776c- History Museum. Natural History Museum (London) 7568-4876-b800-753f527569ec Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1908.1.5.26) from Natural VertNet database http://data.nhm.ac.uk/object/d580fe99- History Museum. Natural History Museum (London) 1e9d-4540-9ca0-2e10ab18dcc3 Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1912.9.23.14) from Natural VertNet database http://data.nhm.ac.uk/object/ec073508- History Museum. Natural History Museum (London) ceae-4688-880c-3fd00f855aaf Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1988.17.1) from Natural VertNet database http://data.nhm.ac.uk/object/141067d4- History Museum. Natural History Museum (London) 2105-438d-b128-2388ffd3670d Collection Specimens (published on 2015-10-10) Locality data of crested ibis (1988.17.2) from Natural VertNet database http://data.nhm.ac.uk/object/21cd5738- History Museum. Natural History Museum (London) ae0f-496c-b798-8fa4b60715f7 Collection Specimens (published on 2015-10-10) Locality data of crested ibis from University Museum VertNet database http://ipt.vertnet.org:8080/ipt/resource.do? of Zoology Cambridge (Zoology). UMZC Zoological r=umzc_vertebrates Specimens (published on 2015-02-16) Locality data of crested ibis from University of VertNet database http://ipt.vertnet.org:8080/ipt/resource.do? Michigan Museum of Zoology. UMMZ Birds r=ummz_birds Collection (published on 2015-10-28) WorldClim v1.4 [52] http://www.worldclim.org/current Software and Algorithms PALEOMIX [53] https://geogenetics.ku.dk/publications/ paleomix; RRID: SCR_015057 Picard v1.92 N/A http://broadinstitute.github.io/picard/; RRID:SCR_006525 ANGSD v0.615 [11] https://github.com/ANGSD/angsd RapidNJ v2.3.2 [54] http://birc.au.dk/software/rapidnj/ R package ‘pcadapt’ [13] https://cran.r-project.org/web/packages/ pcadapt/index.html STRUCTURE v2.3.4 [14] https://web.stanford.edu/group/pritchardlab/ structure.html; RRID:SCR_002151 PSMC [16] https://github.com/lh3/psmc PopSizeABC v2.1 [17] https://forge-dga.jouy.inra.fr/projects/ popsizeabc/ Maxent v3.4.0 [55] https://github.com/mrmaxent/Maxent GRASS GIS v7.2.1 N/A https://grass.osgeo.org/grass7/ R v3.2.2 [56] https://www.r-project.org NeEstimator v2.1 [57] http://www.molecularfisherieslaboratory. com.au/neestimator-software/ (Continued on next page) e3 Current Biology 29, 340–349.e1–e7, January 21, 2019 Continued REAGENT or RESOURCE SOURCE IDENTIFIER Genome-wide Complex Trait Analysis (GCTA) v1.91.4 [58] http://cnsgenomics.com/software/gcta/ PLINK v1.07 [59] http://zzz.bwh.harvard.edu/plink/; RRID:SCR_001757 GERMLINE v1.5.1 [60] http://gusevlab.org/software/germline/; RRID:SCR_001720 R package ‘pegas’ v0.10 [61] https://cran.r-project.org/web/packages/ pegas/index.html SnpEff v4.3 [62] http://snpeff.sourceforge.net/SnpEff.html; RRID:SCR_005191 GERBIL v1.1 [63] https://github.com/dice-group/gerbil R package ‘PopGenome’ v2.6.0 [64] https://cran.r-project.org/web/packages/ PopGenome/index.html CONTACT FOR RESOURCE SHARING Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Guojie Zhang (guojie.zhang@bio.ku.dk). EXPERIMENTAL MODEL AND SUBJECT DETAILS The current study compares whole genome sequencing data generated from 57 historic crested ibis (22 samples from Northeast Asia, 11 samples from East China, and 24 from Northeast China) against 8 previously published contemporary genomes, to charac- terize the temporal changes of genetic diversity for this iconic endangered species. All known information on the context and sequencing information of the samples is provided in Data S1. All researches are under the oversight of the institutional review board. METHOD DETAILS Whole-genome sequencing of museum samples DNA was extracted at the Centre for GeoGenetics from preserved historic skin and toepad samples from 57 crested ibis. Tissue was macerated in extraction buffer and digested overnight according to a modified silica extraction protocol of Rohland [65]. Briefly, buffer (pH8.0, 0.5M EDTA, 1% SDS, 50mM DTT, 1mg/ml, ProteinaseK) was used with samples and negative controls (in a ratio of 16:1) for overnight digestion at 55 C. Samples and controls were then spun down and the supernatant was transferred to filter columns (Amicon, 30KDa cutoff). After concentration to sub 1ml volumes, concentrates were passed through MinElute columns (QIAGEN, Venlo, Netherlands) according to the manufacturer’s protocol, and finally eluted in 60ml EB buffer. This DNA was con- verted into Illumina sequencing libraries through blunt-end ligation using the NEB E6070 kit and the slightly modified version of the protocol of Vilstrup et al. [66], reported by Paijmans et al. [67]. Each DNA tissue extract was used to create one library. Each unique library included individual barcodes and were sequenced by Illumina HiSeq2000 using a TruSeq SBS sequencing kit version 3 (Illumina, San Diego, CA) for 100 cycles from each end of the fragments. Reads were analyzed with Casava1.8.2 (Illumina, San Diego, CA). In summary, we obtained 1028 Gb of sequencing for the 57 historical samples (Data S1). Read alignment and SNP calling We used the identical bioinformatical procedures for the raw reads of all 65 study individuals, including 57 museum individuals in this study and 8 contemporary samples obtained from the previous study (PRJNA308878) [10] (Data S1). All raw sequencing reads were put into PALEOMIX [53] for the adaptor removal, the crested ibis reference genome mapping (PRJNA308878) [10], duplicate filtering, initial BAM rescaled, etc. In order to improve the quality of the initial alignments, the PALEOMIX pipeline incorporates several useful steps: Picard tools for filtering, mapDamage2.0 for rescaling (only applied for the museum samples), and local re-alignment. After that, the mean sequencing depth estimated from the BAM files, ranged from 3.83-10.13x (museum samples) to 20.05-26.71x (contemporary samples), with an average depth of 7x over museum individuals and 21.64x over contemporary individuals (Data S1). Further to avoid any bias in SNP calling due to the different sequencing depths, we randomly extracted data to a coverage of only 7x from the original BAM files of the contemporary samples. Current Biology 29, 340–349.e1–e7, January 21, 2019 e4 We identified SNPs using Analysis of Next Generation Sequencing Data (ANGSD) v0.615 [11] based on a BAM list containing 57 BAM files of the museum samples and 8 extracted BAM files of the contemporary samples with the following parameters: doMajorMinor 1 -GL 2 -doMaf 2 -doGeno 7 -SNP_pval 1e-6 -doPost 1 -doCounts 1 -doglf 4 -dumpCounts 4 -HWE_pval 1 -minMaf 0.01 Considering the low depth of sequencing data caused by the features of museum samples, we further added two steps to reduce the false positives: 1. ANGSD provides the genotypes for each individual. We removed the singleton SNPs. 2. If the genotype of one individual is assigned as major/major, we marked its genotype as NN when the depth of its major allele is less than 2. If genotype of one individual is assigned as minor/minor, we will mark its genotype to NN when the depth of its minor allele is less than 2. If the genotype of one individual is assigned as major/minor, we will mark its genotype to NN when neither the depth of its major nor minor allele is less than 2. Post-processing, a total of 5,268,206 SNPs were retained for further analysis. QUANTIFICATION AND STATISTICAL ANALYSIS Data pre-processing Since background linkage disequilibrium (LD) generated by genetic drift could generate spurious clustering [14, 68], we first need to avoid such bias. Given that the high missing ratio of the loci will influence the accuracy of the LD value estimation, in the following analyses we only used the loci for which the missing ratio is less than 20%. We first ran a pilot test to select the appropriate distance cutoff for pairwise LD calculation using PLINK v1.07 [59]. In our dataset, when the distance of neighboring SNPs is beyond 500k, 97.2% LD values are less than 0.2 (very weak LD values), which suggests there are very few localized SNP clusters on this scale. Furthermore, following published guidance [68], we set the cutoff for the high LD value to 0.5. We used the following steps to filter the data. On each chromosome/scaffold, we started at the first SNP locus, measured the pairwise LD successively, until we reached the first SNP for which LD was less than 0.5. We then retained the first and last SNP loci, then started this process again from the last SNP until all SNP loci had been scanned. We finally obtained a dataset of 3,246,559 SNPs to be used in the population structure analyses. Population structure We used the processed SNP set to construct a phylogenetic tree for all 65 sequenced samples using RapidNJ v2.3.2 [54] based on the neighbor-joining method [12]. The core of this method is to calculate the distance matrix of Dij between each pair of individuals (i and j) as the following formula: X M  Dij = dij L m=1 M is the number of segregating sites in i and j L is the length of regions dij is the distance between individuals i and j at site dij = 0, when individuals i and j are both homozygous for the same allele (AA/AA) dij = 0.5, when at least one of the genotypes of an individual i or j is heterozygous (Aa/AA, AA/Aa or Aa/Aa) dij = 1, when individuals i and j are both homozygous but for different alleles (AA/aa or aa/AA) We also identified the groups/subgroup structure of all 65 individuals by principal component analysis (PCA) of biallelic SNPs using the R package ‘pcadapt’ v3.0.4 [13] in R v3.2.2 [56]. We chose 4 as the optimal K from the PCA by scree plot, because the eigenvalue departs from the straight line. PCA result is consistent with that observed in the NJ tree. The individual ancestries of all 65 samples were inferred by the Bayesian inference program STRUCTURE v2.3.4 [14] after removing the missing loci. Analyses were run five times for each level of K (from 2 through 10), with a 40,000 Burn-in period and 15,000 MCMC repetitions using the admixture model without prior information about populations. The ad hoc statistic DK, based on the rate of change in the log probability of data between successive K values, was used to find the optimal number of subpop- ulations [15] (Data S2). Population demographic analysis To gain insights into the demographic history of the historical and contemporary populations, we first performed the PSMC analysis [16]. We used the heterozygous SNP loci with MAF > = 0.2 and sequencing depth > 2X and < 30X on autosomal sequences. Accord- ing to our previous experience [69], we set the parameters for the PSMC to be ‘‘N30 –t5 –r5 –p 4+30*2+4+6+10,’’ and performed bootstraps (100 times) to represent the variance of results. To further reconstruct recent demographic history (from the sampling time points to 10 kya), we ran an approximate Bayesian computation approach named PopSizeABC v2.1 [17] using the same input e5 Current Biology 29, 340–349.e1–e7, January 21, 2019 of the PSMC analysis, with the following parameters: mac (minor allele count threshold for AFS and IBS statistics computation) is 0; mac_ld (minor allele count threshold for LD statistics computation) are 2,3,5,9 separately; L (size of each segment, in bp) is 4000000; nb_rep (number of simulated datasets) is 500; nb_seg (number of independent segments in each dataset) is 30; all others as default. Ne of the three groups (EA, NW and contemporary group) at the sampling time point was estimated by NeEstimator v2.1 [57], and according to the NJ tree, the Ne of three NE subgroups were calculated respectively using the same method (Table S1). Ecological niche models analysis There are 35 occurrences of crested ibis assigned to clear locality in the 19th and 20th centuries in the VertNet database, which can be used for the breeding ranges predication. We put the latitude and longitude to these 35 occurrences based on their identified local- ities into this analysis. Subsequently these geographic coordinates were used to construct ecological niche models (ENMs) for the contemporary period (19th-20th centuries), the Last Glacial Maximum period (LGM, approximately 22 kya) and the Last Interglacial period (LIG, approximately 120-140 kya). To build the ENMs, we also used the bioclimatic variables from WorldClim v1.4 [52] as the environmental layers in Maxent v3.4.0 [55]. The spatial resolution for all bioclimatic variables was unified to 2.5 arc-minutes (approximately 4 km 3 4 km). For the LIG period, we aggregated the original 30 arc-second data to above unified resolution by func- tion r.resamp.stats in GRASS GIS v7.2.1 with new cell size 0.041666666667 and average aggregation method. For all 19 bioclimatic variables, the predictions were performed in 10 replicates using Maxent, with extrapolation disabled and the other parameters set to default. Temporal population polymorphism In order to comprehensively characterize the population polymorphism, we used three methods. First, we estimated nucleotide di- versity (p) for all 65 individuals in each 5Mb window across the chromosomes using the nuc.div function within the R package ‘pegas’ v0.10 [61]. Second, we measured the extent of the polymorphism fixation process in the contemporary group, when restricting the number of individuals with missing genotypes to be fewer than 5 in the historical groups, and no missing individuals in the contem- porary group. Third, we identified the Identity by Descent (IBD) regions using GERMLINE v1.5.1 [60] on each chromosome in the historical and contemporary groups respectively. Correlation between inbreeding coefficient and effective population size (Ne) As we know that populations with smaller Ne tend to have proportionally more individuals with higher inbreeding coefficients than populations with larger Ne, we explored whether the historical groups with the relative larger Ne had also suffered from inbreeding, and whether the inbreeding coefficient of the contemporary group is significantly higher than that of other species who have a similar Ne. We therefore collected 14 published population SNP datasets [43–51], estimated their inbreeding coefficients and Ne, and used these data to build the correlation relationship. The inbreeding coefficient (FUNI) was estimated based on the correlation between unit- ing gametes following Wright [24] using Genome-wide Complex Trait Analysis (GCTA) v1.91.4 [58]. The Ne was estimated by NeEstimator v2.1 [57]. All results can be found in Table S2. Furthermore, the linear regression function (lm) in R v3.2.2 was applied to build and examine the correlation between FUNI and Ne, and was also used to determine whether the values of the contemporary group deviated significantly from the regression equation. The F-statistic of the coefficient indicated that FUNI has a significant cor- relation with Ne (p value = 0.0000697), and the adjusted R squared (0.7891) also supported the high credibility of the regression equation. By using this regression equation, we found that when Ne was set to 10.7 (Ne of the contemporary group), the FUNI fit was 0.07 (95% CIs: 0.05-0.10). However, the FUNI of the contemporary group was 0.17, verifying its significant deviation from the regression equation. Patterns of deleterious mutations across populations We used two different measures for estimating the deleterious load: deleterious missense mutations, and loss of function (LOF) mu- tations. Since these two measures need the ancestral state for the crested ibis for comparison, we used the genotypes of major al- leles in the historical groups to represent the ancestral state. In each measure, we compared the counts of 1) only homozygous sites (2 per each site) and 2) both homozygous and heterozygous sites (2 per homo- and 1 per heterozygous site). The first measure we used to diagnose the deleteriousness of missense mutations is the Grantham Score [26], a measure of physical/chemical properties of amino acid changes. The score ranges from 5 to 215, with scores bigger than 150 designated as radical or deleterious [70]. We further categorized coding mutations into 1) deleterious, 2) benign, and 3) synonymous mutations as follows: 1) deleterious: missense mutations with Grantham Scores equal to or more than 150. 2) benign: missense mutations with Grantham Scores less than 150. 3) synonymous: synonymous mutations, which do not cause amino acid sequence change. The second measure is the LOF mutations inferred for the coding region SNPs of each sample using the SnpEff v4.3 algorithm [62]. We built a database from the crested ibis annotation and the reference crested ibis genome sequence, and used the input files in VCF format to annotate SNPs and assigned a SnpEff category to the input SNPs per individual. Briefly, two categories, stop codon (nonsense) and splice site disrupting single-nucleotide variants (SNVs), are considered as the LOF mutations. The ratio of homo- to Current Biology 29, 340–349.e1–e7, January 21, 2019 e6 homo- and heterozygous LOF sites in the contemporary population is significantly higher than that in the historic populations (0.336 ± 0.042 in contemporary group, 0.253 ± 0.038 in historical group, Welch Two Sample t test, p = 0.0005386, Figure S2F). MHC diversity analysis A previous study exploring the MHC genomic organization and structure showed that the crested ibis contains 54 genes within three regions spanning 500 kb [71]. We downloaded the following three region sequences from NCBI (Class I 150 kb KP182409.1, Class II 85 kb KP182408.1, and Extended Region 265 kb KP182407.1). We produced the MHC SNPs dataset for all 65 individuals by mapping the raw sequencing data against three regions using the same method as hereinbefore described. Since the genome is organized into blocks of haplotypes [72, 73], we compared the extent of the diversity changes in MHC region based on the haplotype block compo- sition. The software GERBIL v1.1, a tool in the GEVALT software suite [63], was used to capitalize on the greater resolution of genomic markers, and to aggregate SNPs into discrete haplotype blocks. Haplotype blocks were detected based on the MHC SNPs dataset of 57 historical samples after filtering the loci with more than half missing individuals. The historical MHC SNPs dataset was first con- verted to the format described in the GERBIL software manual. Briefly, individual tab-delimited files were produced for each of the scaffolds, with columns representing SNPs and every two consecutive rows corresponding to one sample. SNPs were represented in 0, 1 format or a ‘‘’’? if missing. Re-formatted files were analyzed in an iterative fashion using default software parameters (sample code: gerbil.exe < input filename > ). We also identified the haplotype blocks based on the SNPs dataset of 8 contemporary samples after filtering the loci with more than half missing individuals using the same way. In this way, haplotype struc- tures in two groups were independent, and we used the length and the allelic richness of haplotype blocks as the indexes to feature the haplotype structures and the genetic diversity, respectively. Since the different sample size should be taken into account when comparing the genetic diversity between the historical and contemporary groups, we picked up 8 samples from the historical groups 10 times randomly, counted the number of alleles for each time and used the average values to represent the genetic diversity of the historical groups (Figure 4, Figure S3). Natural selection analysis To measure the selection across the coding regions, we used the R package ‘PopGenome’ v2.6.0 [64] to calculate Tajima’s D value [74], since a positive value is indicative of balancing selection or population contraction. We divided all SNPs into the 3 historical groups (NW, NE and EA) and the contemporary group, and only chose the loci without missing data in each group as input data. To reduce the influence of demographic history when interpreting the Tajima’s D results, for each group we selected genes with the top 5% highest Tajima’s D values, and considered those genes that were present in at least two of the historic populations as candidate genes that were under balancing selection at least 100 years ago. DATA AND SOFTWARE AVAILABILITY Raw sequence read data have been deposited into the study accession number PRJNA481495 and are also available in the CNGB Nucleotide Sequence Archive under accession number CNSA: CNP0000084. e7 Current Biology 29, 340–349.e1–e7, January 21, 2019