\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Not out of the woods yet: signatures of the prolonged negative genetic
consequences of a population bottleneck in a rapidly re-expanding wader,
the black-faced spoonbill \emph{Platalea minor}}



\author[1]{Shou-Hsien Li}%
\author[2]{Yang Liu}%
\author[1]{Chia-Fen Yeh}%
\author[3]{Yuchen Fu}%
\author[3]{Carol K. L. Yeung}%
\author[3]{Chun-cheng Lee}%
\author[1]{ Chi-Cheng Chiu}%
\author[3]{Tung Hui Kuo}%
\author[4]{FT Chan}%
\author[3]{Yu-Chia Chen}%
\author[3]{Wen-ya Ko}%
\author[3]{Chen-te Yao}%
\affil[1]{National Taiwan Normal University}%
\affil[2]{Sun Yat-sen University}%
\affil[3]{Affiliation not available}%
\affil[4]{Endemic Species Research Institute}%


\vspace{-1em}



  
  \date{May 16, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
The long-term persistence of a population which has suffered a
bottleneck partly depends on how historical demographic dynamics
impacted its genetic diversity and the accumulation of deleterious
mutations. Here we provide genomic evidence for the detrimental genetic
effect of a recent population bottleneck in the endangered black-faced
spoonbill (\emph{Platalea minor}) even after its rapid population
recovery. Our population genomic data suggest that the bird's effective
population size, \emph{N\textsubscript{e}}, had been relatively stable
(7,500-9,000) since the end of the last glacial maximum; however, a
recent brief yet severe bottleneck (\emph{N\textsubscript{e}}= 20)
around the 1940s wiped out more than 99\% of its historical
\emph{N\textsubscript{e}} in roughly three generations. By comparing it
with its sister species, the royal spoonbill (\emph{P. regia}) whose
conservation status is of lesser concern, we found that despite a more
than 15-fold population recovery since 1988, genetic drift has led to
higher levels of inbreeding (7.4 times more runs of homozygosity longer
than 100 Kb) in the black-faced spoonbill than in the royal spoonbill
genome. Although the two spoonbills have similar levels of genome-wide
nucleotide diversity and heterozygosity, because of relaxed purifying
selection, individual black-faced spoonbills carry 3\% more
nonsynonymous substitutions than royal spoonbills each of which is 7\%
more deleterious. Our results imply that the persistence of a threatened
species cannot be inferred from a recovery in its population. They also
highlight the necessity of continually using genomic indices to monitor
its genetic health and employing all possible measures to assure its
long-term persistence in the ever-changing environment.%
\end{abstract}%



\sloppy


\textbf{Introduction}

Human activities are driving up to one million species to
extinction(IPBES, 2019; Pimm et al., 2014). All threatened species are
characterized by small population sizes and a declining population trend
or have been through population bottlenecks. Theories suggest that when
a population becomes small or has gone through a population bottleneck,
the drastically augmented random fluctuation of allele frequencies over
time (genetic drift) could lead to the loss of its genetic variation,
i.e. genetic erosion. This would compromise a species' potential to
evolve in response to the ever-changing environment(Barrett \& Schluter,
2008; Lai et al., 2019; Lande \& Shannon, 1996; Bijlsma \& Loeschcke,
2012), and lower the efficacy of purifying selection in removing
deleterious genetic variants(Kimura, 1962; Kirkpatrick \& Jarne, 2000)
resulting in the accumulation of deleterious mutations in it (mutation
load(Kimura, 1962; Kimura, Maruyama, \& Crow, 1963). It can also elevate
levels of inbreeding that increase homozygosity and the expression of
deleterious recessive alleles, thereby reducing individual fitness
(Charlesworth \& Charlesworth, 1987; Keller \& Waller, 2002). Moreover,
theories suggest that such detrimental consequences may persist even
after a population re-expands (Kirkpatrick \& Jarne, 2000). Therefore,
further conservation measures would be required to assure the
persistence of a re-expanded endangered population. Although a low level
of genetic variation and accumulated deleterious mutations were found in
endangered or bottlenecked populations (Robinson et al., 2016; Grossen,
Guillaume, Keller, \& Croll, 2020), results of some genomic studies
suggest otherwise: bottlenecked non-African human populations do not
have lower genetic variation, and the effectiveness of purifying
selection to remove deleterious mutations is not compromised in the
European population (Do et al., 2015; Fu, Gittelman, Bamshad, \& Akey,
2014; Lohmueller, 2014; Simons, Turchin, Pritchard, \& Sella, 2014).
This implies that further conservation provisions might not be essential
for the long-term survival of a bottlenecked population. Therefore, the
ability of a threatened species to persist may partly depend on the
extent to which the historical bottleneck event caused a decay in its
genetic diversity and an increase in its mutation load.

The once critically endangered migratory wader, the black-faced
spoonbill (\emph{Platalea minor} )(BirdLife International, 2017), could
serve as an ideal system to assess the genetic legacy of a bottleneck
event in a re-expanding population. \emph{P minor} winters in the
coastal salty wetland habitats of the East Asian coast (Fig. 1A). Its
extant breeding colonies are mainly located on uninhabited rocky islets
along the western and eastern coast of the Korean Bay (Chong, Pak, Rim,
\& Kim, 1996; Ding, Lei, Yin, \& Liu, 1999; Ueta et al., 2002), and
recently expanded to the coast along the Northern Sea of the Japan Basin
(Litvinenko \& Shibaev, 2005; Shibaev, 2010). This species was described
as `common' in documents pre-dating the 1950s (La Touche, 1931; Austin,
1948). However, only 288 individuals were counted within its entire
range in 1988 (Kennerley, 1990). Since then, its population size has
increased remarkably to 4,864 individuals according to a census in the
winter of 2020 (Yu, Li, Tse, \& Fong, 2020). A study of mitochondrial
diversity suggested that the spoonbill had recently experienced a severe
bottleneck (Yeung, Yao, Hsu, Wang, \& Li, 2006). However, without a
documented history of the timing, magnitude and duration of the
bottleneck, we cannot rule out the possibility that its low population
size in the 1980s was part of a natural response to the drastic climate
changes since the end of the last glacial maximum (LGM), as in some
other endangered species (Mays et al., 2018).

To infer the timing and duration of its presumed recent bottleneck
event, we first sequenced and obtained a draft assembly of the
black-faced spoonbill genome, then re-sequenced the whole-genome of
multiple individuals collected from its two major wintering sites in
Taiwan and Hong Kong (Yu et al., 2020) (Fig. 1A). We evaluated the
extent of genetic diversity and deleterious mutations accumulated in the
extant black-faced spoonbill population by comparing its population
genomic data to that of its sister species the royal spoonbill (\emph{P.
regia} ) from which it diverged approximately half a million years ago (
Yeung et al., 2011). The royal spoonbill is commonly found in the
wetlands of Australia and nearby islands (Fig. 1A), and its conservation
status is in the `least concern' category; it does not have a documented
history of bottlenecks (Matheu \& del Hoyo, 1992). We specifically
addressed the following for the black-faced spoonbill: (1) whether the
drastic climate change since the end of the LGM significantly impacted
its population trajectory before the presumed recent population
bottleneck; (2) what the start time, duration and magnitude of the
presumed recent bottleneck were; (3) whether the recent bottleneck
event, if any, has led to a relaxation of selection and a higher level
of inbreeding, genetic drift and the accumulation of deleterious
mutations than in the royal spoonbill.

\textbf{Results}

\emph{De Novo} \textbf{assembly of a draft genome for the black-faced
spoonbill}

We sequenced and conducted a \emph{de novo} genome assembly to generate
a reference genome from a rescued male black-faced spoonbill. The
assembled genome is 1.19 Gb in length (genomic coverage = 42.5\selectlanguage{ngerman}×) and
contains 2,243 scaffolds (870 scaffolds \textgreater{} 10 Kb, total
1.186 Gb; scaffold N50 of 4.2 Mb; 34,176 contigs; contigs N50 of 71.0
Kb; Table S1). This draft genome includes a total of 78.2 Mb (6.6\% of
total genome) of repetitive sequences (Table S2). We found 25,315
predicted genes, 11,733 of which could be functionally annotated. In
total, we obtained 399.6 Mb of genic regions (33.6\% of the total
genome) consisting of 365.2 Mb of introns and 34.4 Mb of coding regions.
The completeness of our draft genome is high (Table S3): we identified
93.0\% (4,569) of the Benchmarking Universal Single-Copy Orthologs
v3.0(Waterhouse et al., 2018). Therefore, the completeness of our draft
genome is comparable to other published avian genomes, such as the draft
genome of horn lark (Mason, Pulgarin, Cadena, \& Lovette, 2020).

\textbf{Comparable genomic diversity between the black-faced spoonbill
and royal spoonbill}

We re-sequenced the genomes of 11 black-faced spoonbill individuals,
nine of which were from Tainan County, Taiwan, the wintering ground of
more than 50\% of the species' world population (Chesser, Yeung, Yao,
Tian, \& Li, 2010) and two from the Deep Bay Wetland, Hong Kong (Fig.
1A; Table S4). We also used the draft genome of the black-faced
spoonbill as the reference to re-sequence the whole genomes of nine wild
royal spoonbill individuals from Australia (Table S4). The genomic
diversities of the two spoonbill species are comparable: the nucleotide
diversity (\selectlanguage{greek}\emph{π} \selectlanguage{english}) and mean individual heterozygosity of the
black-faced spoonbill are 1.11\selectlanguage{ngerman}× 10\textsuperscript{-3} and 1.09×
10\textsuperscript{-3}, respectively. The comparable figures are 1.06×
10\textsuperscript{-3} and 1.05× 10\textsuperscript{-3} in the royal
spoonbill (Table S5). The genome-wide heterozygosities of both spoonbill
species are relatively low (Table S6) compared with those of most of the
32 avian species (Lai et al., 2019; Li et al., 2014) in the ``Least
Concern`` category in the IUCN Red List (mean± S.D.= 2.50± 1.07×
10\textsuperscript{-3}), but twice that of the critically endangered
crested ibis (\emph{Nipponia nippon} (S. Li et al., 2014); = 0.43×
10\textsuperscript{-3}), and close to that of seven vulnerable
species(S. Li et al., 2014) (1.18± 0.85× 10\textsuperscript{-3}). Our
results suggest that the genetic diversity of black-faced spoonbills is
still comparable to those of the un-threatened royal spoonbill and other
vulnerable but not endangered avian species. The genome-wide genetic
distance (\emph{D} \textsubscript{xy}) between the two spoonbills was
estimated to be 0.0020. This is consistent with the short mitochondrial
DNA distance (the Tamura-Nei distance= 0.015) observed between these two
sister species (Chesser et al., 2010).

\textbf{Genetic panmixia of black-faced spoonbills sampled in Taiwan and
Hong Kong}

Migratory birds with different wintering grounds might have fidelity to
their own breeding grounds which could lead to genetic differentiation
(Rubenstein et al., 2002). To test whether black-faced spoonbills
wintering in Taiwan and Hong Kong originated from different genetic
groups, we conducted a model-based admixture analysis (Alexander \&
Lange, 2011) based on 215,722 unlinked (\emph{r}
\textsuperscript{2}\textless{} 0.2) genome-wide autosomal SNPs (coverage
\textgreater{} 10×). The results of the admixture analysis (Fig. 1B) and
principal components analysis (Fig. 1C) suggest that all black-faced
spoonbill individuals sampled in the current study likely consisted of a
single breeding population, allowing us to combine the genomic data from
the individuals from the two wintering sites for the subsequent
analyses. Our results are consistent with satellite tracking records
which showed black-faced spoonbills tagged in both Taiwan and Hong Kong
returning to the Demilitarized Zone of Korea in the breeding season
(Ueta et al., 2002).

\textbf{Relatively stable post-LGM demographic histories of the
black-faced and royal spoonbills}

We reconstructed the historical demographic trajectories of the two
sister spoonbill species with \emph{SMC} ++ (Terhorst, Kamm, \& Song,
2017). The \emph{N} \textsubscript{e} of the black-faced spoonbill was
estimated to be about 8,900 individuals at the LGM (approximately 22,000
years ago), which then gradually declined to 7,432 individuals until
around 435 years ago (Fig. 2A). Conversely, the royal spoonbill had
an\emph{N} \textsubscript{e} that gradually increased from approximately
7,800 individuals at the LGM to 8,130 individuals around 500 years ago
(Fig. 2A). The similar post-LGM \emph{N} \textsubscript{e} trajectories
of the two species also make the royal spoonbill an ideal control group
to infer the effect of a recent bottleneck on the genome of the
black-faced spoonbill.

\textbf{A brief yet severe recent bottleneck for the black-faced
spoonbill}

To test whether the black-faced spoonbill experienced a recent
bottleneck event, we compared the maximum composite likelihood of a
recent bottleneck model and a constant population size model based on
215,722 unlinked autosomal SNPs using \emph{Fastsimcoal2} (Excoffier,
Dupanloup, Huerta-Sánchez, Sousa, \& Foll, 2013). The results
significantly supported the recent bottleneck model (maximum estimated
likelihood= -230,635.54, AIC= 421,681.07, AIC weight= 1.0) instead of a
constant population model (maximum estimated likelihood= -232,073.84;
AIC= 464,149.07, \selectlanguage{greek}\emph{Δ} \selectlanguage{english}AIC= 2,868.6, AIC weight= 0.0). The presumed
bottleneck event was estimated to be very recent and brief (Fig. 2B),
starting approximately six generations ago (95\% CI: 4-12 generations
ago) and ending about three generations ago (95\% credible interval, CI:
2-4 generations ago). Because most black-faced spoonbill individuals (9
of 12 samples) sequenced in this study were victims of a botulism event
in 2002, and assuming a generation time of 10 years (Yeung et al.,
2006), our results suggest that the recent bottleneck event occurred
around the 1940s.

\textbf{An elevated level of inbreeding in the black-faced spoonbill}

Small populations are particularly prone to an elevated level of
inbreeding, the mating between relatives, which can increase
homozygosity at loci that carry deleterious recessive mutations and thus
reduce the fitness of offspring (Crow \& Kimura, 1970) (i.e. inbreeding
depression). Two lines of evidence reveal the genomic footprints of
inbreeding during the bottleneck and possibly the post-bottleneck period
in the black-faced spoonbill population. A recent and severe bottleneck
event is expected to result in an increase in the number (NROH) and
length (SROH) of runs of homozygosity, which are identical homologous
genomic segments(McQuillan et al., 2008). We found that the black-faced
spoonbill had approximately 7.4 times the SROH (\emph{t} test=
-12.68,\emph{p} \textless{} 0.0001) and 6.2 times the NROH (\emph{t}
test= -22.78, \emph{p} \textless{} 0.0001) of the royal spoonbill (Fig.
3A; Table S7): the black-faced spoonbill had an average NROH of 426.42\selectlanguage{ngerman}±
49.65 and average SROH of 84,154.15± 19,442.36 Kb; by contrast, the
royal spoonbill had an average NROH of 68.89± 19.21 and an average SROH
of 11361.82± 3645.35 Kb. In both species, most of the ROHs were less
than 500 Kb long (84.6\% and 90.4\% of SROH for the black-faced and
royal spoonbill, respectively; 96.3\% and 97.6\% of NROH for the
black-faced and royal spoonbills, respectively), yet significantly
longer ROHs (\textgreater{} 2,000 Kb) were observed in two black-faced
spoonbills with the longest being 3,566.6 Kb. The genomic inbreeding
coefficient, \emph{F} \textsubscript{ROH}, which measures the fraction
of the autosomal genome in ROH ( Keller, Visscher, \& Goddard, 2011),
was approximately seven times higher in the black-faced spoonbill than
in the royal spoonbill (mean \emph{F} \textsubscript{ROH}= 0.071 and
0.001 for the black-faced and royal spoonbills, respectively), but lower
than those of other threatened populations (e.g. \emph{F}
\textsubscript{ROH}= 0.4 in the grey wolf (Gómez\selectlanguage{english}-Sanchez et al., 2018)).
This may reflect alleviated inbreeding due to the rapid population
recovery. Moreover, inbreeding could also lead to a higher level of
homozygosity and slower decay in linkage disequilibrium (Slatkin, 2008)
which can be measured by the association of occurrence of alleles at two
loci, expressed as a correlation coefficient (\emph{r}
\textsuperscript{2}). In this case,\emph{r} \textsuperscript{2} decay
was slower in the black-faced spoonbill than in the royal spoonbill
(Fig. 3B).

\textbf{An elevated level of genetic drift in the black-faced spoonbill}

A bottleneck event could intensify the effect of genetic drift,
resulting in an increased rate of extinction and therefore a deficiency
of low-frequency alleles as well as an over-abundance of medium
frequency alleles (Tajima, 1989; Gattepaille, Jakobsson, \& Blum, 2013;
Maruyama \& Fuerst, 1985). To measure the strength of genetic drift in
the black-faced spoonbill, we calculated the folded site frequency
spectrum (SFS) of unlinked SNPs of the autosomal chromosomes for the two
spoonbill species. The folded SFS of the two spoonbill species were
bimodal distributions (Fig. 3C) with a tall peak for singletons (minor
SNPs carried by only one haplotype) and a lower peak for
medium-frequency SNPs. However, the frequency of singletons in the
black-faced spoonbill (21.9\%; 47,227 of 215,722 SNPs) was significantly
lower than that in the royal spoonbill (30.7\%;\selectlanguage{greek}\emph{χ}
\selectlanguage{english}\textsuperscript{2}= 24.82, \emph{p} = 0.000). In theory, genome-wide
Tajima's \emph{D} values would be higher with fewer low-frequency SNPs
(Tajima, 1989) (e.g. singletons), as was the case in the black-faced
spoonbill (mean Tajima's \emph{D} = 0.24\selectlanguage{ngerman}± 0.98, the number of total
SNPs, \emph{n} = 4,522,032) compared with the royal spoonbill (mean
Tajima's \emph{D} = 0.05± 0.99, \emph{n} = 4,437,044;\emph{t} = 296.51,
\emph{p} \textless{} 0.001; Fig. 3D). The lower singleton frequency and
higher Tajima's \emph{D} values suggest that the black-faced spoonbill
has been subjected to stronger genetic drift caused by the recent
population bottleneck.

\textbf{Relaxation of selection in the black-faced spoonbill}

The elevated level of genetic drift in the black-faced spoonbill might
have compromised the efficacy of natural selection (Kimura, 1962). To
examine whether the elevated level of genetic drift had lowered
selection efficiency in the black-faced spoonbill, we used
\emph{RELAX}(Wertheim, Murrell, Smith, Kosakovsky Pond, \& Scheffler,
2015) to detect genes under relaxed selection in the two spoonbills. Of
the 5,657 well-aligned homologous genes between the two spoonbill
species and another Ciconiformes species, the crested ibis ( Li et al.,
2014),\emph{RELAX} analysis revealed significantly more genes in the
black-faced spoonbill (2,094, 37.0\%) to be under relaxed selection
(relaxation parameter \emph{k} \textless{}1, \emph{p} \textless{} 0.05)
than in the royal spoonbill (1,542, 27.2\%; \selectlanguage{greek}\emph{χ}
\selectlanguage{english}\textsuperscript{2}= 123.04, \emph{p} = 2.2e\textsuperscript{-16};
Figure 4A). Conversely, there was no evidence (\selectlanguage{greek}\emph{χ}
\selectlanguage{english}\textsuperscript{2}= 3.1671.04, \emph{p} = 0.075) to support the
hypothesis that selection was significantly intensified in the genome of
the royal spoonbill (33 and 55 genes with\emph{k} \textgreater{}1,
\emph{p} \textless{} 0.05 for the royal and black-faced spoonbill,
respectively; Figure 4A). Therefore, a higher level of genetic drift had
probably caused relaxed selection in the black-faced spoonbill.

\textbf{Higher mutation load in the black-faced spoonbill}

A consequence of relaxed selection is the accumulation of deleterious
mutations in the genome (Kimura, 1962; Kirkpatrick \& Jarne, 2000). We
found that black-faced spoonbill individuals carried about 3\%
(8,294.00\selectlanguage{ngerman}± 332.66 substitutions) more nonsynonymous substitutions (NS)
than royal spoonbill individuals (8,025.78± 211.77
substitutions;\emph{t} = -1.91, \emph{p} \textless{} 0.035; Fig. 4B;
Table S8). Moreover, the mean phenotypic effect, Grantham's score
(Grantham, 1974), for each deleterious NS (Grantham's
score\textgreater{} 50; Chun \& Fay, 2009)) was approximately 7\% higher
for the black-faced spoonbill (56.98± 0.69) than in the royal spoonbill
(53.25± 1.83; two-sample\emph{t} test= 6.5154, \emph{p} =
3.6e\textsuperscript{-6}; Table S9). As a result, the average
accumulated Grantham's score per individual was higher in the
black-faced spoonbill (527,350± 33,698.9) than in the royal spoonbill
(502,030± 9,571.9; \emph{t} = 2.47, \emph{p} = 0.014; Fig. 4C; Table
S9). Moreover, the ratio of homozygous to heterozygous NS in the
black-faced spoonbill (0.13± 0.04) was also significantly higher than
that in the royal spoonbill (0.11± 0.007; \emph{t} = -1.86,\emph{p}
\textless{} 0.04; Fig. 4D), suggesting that purifying selection has been
less effective in removing deleterious homozygotes in the bottlenecked
black-faced spoonbill. Therefore, we provide a real-world example of the
theory that mutation load, accumulation of deleterious mutations, is
insensitive to recent population growth (Simons et al., 2014).

\textbf{Discussion}

Relatively stable post-LGM demography

The historical demographic trajectory retrieved from population genomic
data (Fig. 2A) suggests that the \emph{N} \textsubscript{e} of the this
once critically endangered species is much larger than any census global
population sizes since 1988, and was stable for most of its recent
evolutionary histories, unperturbed by the drastic climatic changes
since the end of the LGM. The large and stable historical\emph{N}
\textsubscript{e} of the black-faced faced spoonbill is consistent with
the fact that its genome-wide genetic diversity (\selectlanguage{greek}\emph{π} \selectlanguage{english}and mean
individual heterozygosity) is similar to that of the royal spoonbill
which also had a similar post-LGM demographic trajectory (Fig. 2A).

Another line of evidence which supports the historical abundance of the
black-faced spoonbill as documented in early literature (La Touche,
1931; Austin, 1948) is the fact that that its estimated
long-term\emph{N\textsubscript{e}} could be translated, assuming a
simple one-to-ten rule for the ratio between \emph{N\textsubscript{e}}
and the census breeding population size (Frankham, 1995), into over
74,200 breeding individuals. This is comparable to the similarly deduced
long-term breeding population size of the royal spoonbill of 78,000
individuals and to the current population size of royal spoonbills
(25,000-100,000 individuals; Matheu \& del Hoyo, 1992). A stable and
large long-term \emph{N\textsubscript{e}} rules out the possibility that
the bottleneck observed in 1988 was the result of the long term
population decline driven by climate change since the LGM as in the case
of other endangered species (Mays et al., 2018; Li et al., 2014; Cho et
al., 2013; Zhao et al., 2013).

Could DDT have been the cause of recent population bottleneck?

Although we cannot ascertain the cause of this recent bottleneck, its
estimated start coincides with the period of the widespread use of
organochlorine pesticides (e.g. DDT) around the world, including East
Asia (Wong, Leung, Chan, \& Choi, 2005), after World War II. DDT and its
metabolites are known to impair the reproduction of avian species by
thinning eggshells and acting as estrogen blockers (Cooke, 1973;
Risebrough, 1986). DDT-induced population crashes have been documented
in several raptors (Poole, 1989; Newton \& Wyllie, 1992); eggshells of
other Ciconiformes species (Faber \& Hickey, 1973; King, Meeker, \&
Swineford, 1980) were found to thin significantly in the presence of
even a miniscule amount of DDT. Furthermore, the black-faced spoonbill
is one of the top predator species in its major wintering habitat, the
estuaries, which are known to be repositories of contaminants (Ridgway
\& Shimmield, 2002) where biomagnification through the food chain would
expose the bird to a level of DDT and its metabolites higher than that
found in the environment. The adverse effects of DDT on wildlife
eventually led to its ban around the world (e.g. in 1972 by the United
States), and has been associated with the rapid recovery of some raptors
(Grier, 1982). The bottleneck event of the black-faced spoonbill was
estimated to have ended two to three generations ago, which is around
1970-1980s, consistent with the prohibition of the use of DDT in East
Asia in the 1980s (China in 1983, South Korea in 1986, and Taiwan in
1987) and the black-faced spoonbill's population growth trend since
1988.

The recent bottleneck event experienced by the black-faced spoonbill
was, although brief, extremely severe: the \emph{N}
\textsubscript{e}during the bottleneck period was estimated to be only
20 individuals (95\% CI: 12.5-47.5 individuals), roughly equivalent to
200 breeding individuals. By contrast, the pre-bottleneck \emph{N}
\textsubscript{e}was estimated to be 5491.5 (95\% CI: 4549.5-- 12310.5
individuals), close to the recent historical \emph{N} \textsubscript{e}
(7,432) estimated by \emph{SMC++} . In other words, more than 99\%
of\emph{N} \textsubscript{e} vanished during the bottleneck event.
Although the \emph{N} \textsubscript{e} of the black-faced spoonbill has
since increased to 236.5 (95\% CI: 57.6 -- 376.0), it is still only
about 4\% of its pre-bottleneck \emph{N} \textsubscript{e}. The
magnitude of\emph{N} \textsubscript{e} reduction would rule out the
Korean War as a candidate cause of population decline (Matheu \& del
Hoyo, 1992) since the uninhabited rocky islets where the bird breeds
would hardly have been major military targets.

The questionable long-term persistence of the black-faced spoonbill

The rapid population recovery of the black-faced spoonbill may preclude
the risk of immediate extinction, although whether its current\emph{N}
\textsubscript{e} of 405 could support the species' long-term
persistence remains questionable, as it has been suggested that a
minimum \emph{N} \textsubscript{e} of 1,000 is required to ensure the
evolutionary potential and long-term persistence of a population
(Frankham, Bradshaw, \& Brook, 2014). Moreover, the coastal wetlands of
East Asia are under unprecedented pressure of reclamation, urbanization
and industrialization (Larson, 2015). For example, a seawall in China
which encloses its coastal wetland to support rapid urbanization and
economic development reached 11,000 km (70\% of mainland China's
coastline) in 2010 (Ma et al., 2014), fundamentally altering the coastal
ecosystem and inevitably reducing the environmental carrying capacity
for the black-faced spoonbill to a level much lower than that in the
pre-bottleneck era. The deterioration of its wintering habitat is
unlikely to have been responsible for the past bottleneck event because
the habitat has deteriorated continually, which is inconsistent with the
recent population recovery (Baillie, Hilton-Taylor, \& Stuart, 2004).
Habitat deterioration would reduce the environmental carrying capacity
of the black-faced spoonbill and constrain the recovery potential and
persistence of this endangered species. Unless stringent conservation
measures are taken to improve or at least maintain its current carrying
capacity, the future of the black-faced spoonbill is still in doubt
despite its recent rapid population growth (Sung, Tse, \& Yu, 2018).

Although our results support the conclusion that the black-faced
spoonbill experienced a recent severe population bottleneck, its
genome-wide genetic diversity (\selectlanguage{greek}\emph{π} \selectlanguage{english}and mean individual
heterozygosity) is comparable to that of the royal spoonbill and other
non-vulnerable avian species (Table S5, S6). This is probably because
the short bottleneck duration and the rapid population recovery only led
to the loss of rare genetic variants. This is what is suggested by the
lower frequency of singleton SNPs in the folded site frequency spectrum
(Fig. 3C) and higher genome-wide Tajima's \emph{D} (\textgreater{} 0)
values in the black-faced spoonbill. This observation is consistent with
a simulation study: a short bottleneck is less likely to cause severe
genetic erosion (Allendorf, 1986). Since standing genetic variation with
variants present at high frequencies could allow a species to respond
more swiftly to selection pressures posed by changing environments
(Hermisson \& Pennings, 2005; Lande \& Shannon, 1996), our results
suggest that the recent drastic population bottleneck event might not
severely dampen the evolutionary potential of the black-faced spoonbill.
Therefore, if its further population recovery can be sustained, the
chance for this once critically endangered wader to cope with future
environmental changes should be approximately the same as for royal
spoonbills and other non-vulnerable species.

Although the extant black-faced spoonbill might not be suffering from
severe genetic erosion, it could still be the victim of an elevated
level of inbreeding (Fig. 3) and mutation load (Fig. 4). However,
despite the fact that these birds carry higher numbers of detrimental
mutations, its population has still grown steadily since 1988 (Yu et
al., 2020). One possible explanation for this is that the fitness effect
of mutations could be environment-dependent (Szafraniec, Borts, \&
Korona, 2001): the high genetic load in the black-faced spoonbill might
not necessarily confer lower population fitness in a benign environment.
However, when the environment becomes more stressful, the impact of
mutations' deleterious effects on fitness could increase (Martin \&
Lenormand, 2006), especially when exacerbated anthropogenic
environmental change is expected to further reduce the habitat quality
of wildlife for the foreseeable future. Therefore, even the black-faced
spoonbill population rebound led to its being down-listed from
critically endangered to endangered in the IUCN Red List, and is
considered a success example of international collaborative conservation
efforts (BirdLife International, 2000), our study suggests that these
birds are not totally out of the woods yet.

Furthermore, our results provide an example of modern genomic sequencing
technologies as a powerful tool to monitor the genetic health of and
predict the evolutionary potential for any species whose conservation is
of concern. In the case of black-faced spoonbills, genomic data suggest
that even after decades of rapid population rebound, its fate is still
in the mists and further conversation measures are needed to assure its
survival given that its habitat will inevitably be subject to tremendous
anthropogenic impacts. Given both theoretical (Robertson \& Waddington,
1960; Hill \& Robertson, 1966) and empirical (Weber, 2004) works
suggesting that selection is more effective in larger populations, rapid
population recovery after a bottleneck event may be critical to restore
the efficacy of selection to purge out accumulated detrimental
mutations. Aiding the re-expansion of the black-faced spoonbill should
be considered as an overarching objective of its conservation management
in near future. All conservation practices, such as habitat restoration
(Dobson, Bradshaw, \& Baker, 1997), supplementary feeding (Ewen, Walker,
Canessa, \& Groombridge, 2015), predator removal (Smith, Pullin,
Stewart, \& Sutherland, 2010), and population augmentation from \emph{ex
situ} population(Bowkett, 2009) should be seriously considered and
employed. In addition, we propose that the level of inbreeding and
mutation load of black-faced spoonbill individuals should be continually
monitored. Such information would provide temporal genomic indices (such
as ROH) to evaluate how its genetic health has been improved by current
conservation practices to assist its population recovery strategy. Taken
together, our results are especially relevant when anthropogenic
environmental changes are increasingly turning incidental natural
processes into catastrophic events which put many wildlife species in
great peril of population reduction.

\textbf{Acknowledgements}

We thank Terry Chesser and the Commonwealth Scientific and Industrial
Research Organization of Australia for providing samples of the royal
spoonbills and Carrie K.W. Ma for providing samplings of black-faced
spoonbills from Hong Kong. We are grateful for the valuable help and
comments by Chih-Ming Hong, Wei Liang, Teng-Chiu Lin, and Benyang Liao.
We are in debt to Alan Watson who has greatly improved the readability
of this manuscript. This work was supported by grants to C.- T. Yao from
the Taiwan Endemic Species Research Institute and to S.-H. Li from the
Ministry of Science and Technology, Taiwan, R.O.C., and Y. Liu by
Forestry Administration of Guangdong Province, China (DFGP Project of
Fauna of Guangdong-202115. It was also kindly supported by endorsement
from Mr. Y.-C. Chen.

\textbf{Author contributions} : S.-H.L., C.-T.Y. and Y.L. designed
research; C.-F.Y., and Y.C.C. performed research; W.-Y.K., T.H.K, and
F.-T. C. contributed new reagents/analytic tools; C.-F.Y., C.C.C., , and
Y.C.C., analyzed data; and S.-H.L., Y.L. C.K.L.Y., C.-C.L., Y.F., and
C.-T.Y. wrote the paper. All authors approved the final submission.

\textbf{Competing interests} : The authors declare that they have no
competing interests.

\textbf{Materials and Methods}

\textbf{Reference genome}

\emph{Sampling, DNA extraction, Genomic library preparation and
sequencing}

In 2015 we used a modified chloroform and LiCl precipitation
protocol(Gemmell \& Akiyama, 1996) to extract gross DNA from the blood
of a rescued male black-faced spoonbill (22.88\selectlanguage{ngerman}°N and 120.20°E) which was
kept in the Taiwan Endangered Species Research Institute.

Libraries with insert size 200 bp, 3 Kb and 5 Kb were constructed (with
the Paired-End DNA Sample Prep Kit for the 200 bp library and Mate Pair
Library Prep Kit v2 for the 3 Kb and 5 Kb libraries) and sequenced on an
Illumina HiSeq 2500 platform by Genomics BioSci \& Tech (Taipei,
Taiwan). In total, 124 Gb raw reads were generated (table S1).

\emph{Assembling the reference genome de novo}

We used \emph{ALLPATHS-LG} version-44099(Gnerre et al., 2011)
(parameters PLOIDY=2 and PHRED\_64=1) to assemble the draft genome of
the black-faced spoonbill. We ran \emph{ALLPATHS-LG} on a workstation
with 32 CPUs (2,199.882 MHz) and 377.8 Gb RAM. (The peak value of in-use
memory was 342.1 Gb.)

We used the correction steps \emph{CleanCorrectedReads}
and\emph{ErrorCorrectJump} in \emph{ALLPATHS-LG} to remove 1.7\% of
paired-end reads and to correct 69.3\% of mate-paired reads with the
criterion of low frequency k-mers (K=25 and 96 for paired-end and
mate-paired reads respectively). The raw reads for assembling were 50.5
Gb (67.8\%) and 28.8 Gb (57.5\%) of paired-end and mate-paired reads,
respectively (table S1). In total 34,176 contigs with N50 size of 71.0
kb were assembled with a total length of 1.18 Gb. Finally, 2,243
scaffolds (N50= 4.2 Mb) (table S1) were concatenated from these contigs.
The draft genome generated by \emph{ALLPATHS-LG} is a diploid genome.
For subsequent analyses, we randomly dropped one of the nucleotides from
each heterozygous SNP to generate a pseudo-haploid reference genome.

\emph{Masking repetitive sequences}

We masked repetitive sequences from the draft genome
with\emph{RepeatMasker} version open-4.0.7, with the options -gccalc and
-species ``aves.'' In total we masked 78.2 Mb (6.6\%) of the genome
sequence.

\emph{Removing contaminating sequences}

We removed potential contaminating sequences in the draft genome by
following the procedure in Jones et al.(Jones et al., 2012). On each
scaffold, a 1 Kb fragment without an unknown base ``N'' was sampled from
each 100 Kb window to check for the presence of contaminants. In total
we sampled 13,105 one kb fragments and blasted them against the
NCBI\emph{nt} database with \emph{BLASTN} (v. 2.6.0; e-value cut-off=
10\textsuperscript{-5}). Fragments of the \emph{BLASTN'} s top hit
species that were neither birds nor reptiles were considered
contaminants. We extracted contigs that contained these fragments from
the scaffolds. In total we found contaminated DNA in 9 scaffolds (153.8
Mb), which included bacteria (N= 2), fishes (N= 1) and mammals (N= 6).
Finally, 2,243 scaffolds were retained in the draft black-faced
spoonbill genome assembly.

\emph{Constructing the putative chromosome}

To infer the putative autosome of the black-faced spoonbill, we
used\emph{Satsuma} v2.1(Grabherr et al., 2010) to align our reference
genome assembly to the chicken genome (GRCg6a; GCA\_000002315.5(Warren
et al., 2017)) with option ``l (minimum alignment length) = 100''.
Scaffolds embedded within other scaffolds were removed. If the same
scaffold was repeatedly aligned to different chromosomes, we kept the
longest alignment. 486 scaffolds with a total length 1.10 Gb (93\% of
total draft genome) of the black-faced spoonbill genome were aligned to
seven macro-chromosomes, 30 intermediate and micro-chromosomes and the
mitochondrial DNA of the chicken genome. We classified unaligned
scaffolds as ``unlocalized'' (Table S4).

\textbf{Genome annotation}

\emph{Gene structure and functional annotation}

To annotate the draft genome of the black-faced spoonbill, we used
RNeasy Mini Kit to extract total RNA from a blood sample from the same
individual used for reference genome sequencing following the
manufacturer's instructions. RNA was sequenced by Genomics BioSci \&
Tech (Taipei, Taiwan) using a TruSeq Stranded mRNA Preparation Kit on an
Illumina NextSeq 500 platform with the average insert size of 180 bp. We
used \emph{Trinity} \emph{v} 2.4.0(Grabherr et al., 2011) to analyze 5.2
GB of RNA-seq reads. We used both \emph{de novo} and genome guide (using
the draft genome assembly as a reference) methods to generate 273 Mb and
171 Mb of RNA assembly, respectively. Then, the two assemblies were
merged with \emph{PASA} (Haas et al., 2003). \emph{TransDecoder}
in\emph{PASA} was used to identify likely coding regions (CDS) with the
CDS of the genome of the crested ibis \emph{Nippon
nippon}(GCF\_000708225.1\_ASM70822v1(S. Li et al., 2014)) as the
reference. Finally, we obtained 190 Mb of spoonbill RNA assembly as the
gene model to train \emph{Augustus v3.2.3} (Stanke \& Morgenstern, 2005)
for gene prediction. \emph{Augustus} predicted 11,733 genes totaling 411
Mb (35\% of the total genome) and 25 Mb of CDS in the spoonbill genome.
We extracted the coding regions of all annotated genes and blasted these
CDS to Gallus\_gallus-5.0 (GCF\_000002315.4) with \emph{BLASTP} (v.
2.6.0; e-value cut-off= 10\textsuperscript{-6}). There were only 10,558
genes that met the e-value cut-off criterion.

\emph{Evaluating the completeness of the draft genome}

We used 4,915 Benchmarking Universal Single-Copy Orthologs (BUSCO
v3.0.1, (Waterhouse et al., 2018)) to evaluate the completeness of the
genome assembly and annotation.

\textbf{Population re-sequencing}

\emph{Population sampling and whole genome resequencing}

In addition to the sample used for reference genome sequencing, we also
re-sequenced eleven black-faced spoonbills whose tissue was collected
from individuals that had been rescued or had accidentally died in
Tainan County, Taiwan (N=) and the Deep Bay Wetland, Hong Kong (N=2), as
well as tissue samples of royal spoonbills from Australia (N=9).

We used a protocol from Gemmel and Akiyam (2006) to extract gross DNA
from these tissue samples for resequencing. A whole genome re-sequencing
library with insert size 50 bp was constructed for each individual. Our
sequencing data was generated using an Illumina HiSeq 4000 platform
(BGI).

\emph{Re-sequencing mapping and variant calling}

A consensus sequence for each individual was generated
using\emph{Vcf-tools} v.0.1.13(Danecek et al., 2011). The average
coverage of the black-faced and royal spoonbills was 15.7× and 16.6×,
respectively.

\textbf{Population structure}

Two approaches were used to examine the population structure of the
black-faced spoonbill as sampled from Taiwan and Hong Kong. Firstly, we
used \emph{ADMIXTURE} v1.22(Alexander \& Lange, 2011), which is based on
a Bayesian clustering method, to infer the number of ancestors of the
black-faced spoonbills sampled in current study. Only unlinked autosomal
SNPs were used for the admixture analysis. We used \emph{PLINK} (Purcell
et al., 2007) to filter out linked autosomal SNPs. The window size to
calculate pair-wise SNP-SNP \emph{r} \textsuperscript{2} was set to 10
Kb, and the step of the sliding window was set to 2 Kb. The threshold
of\emph{r} \textsuperscript{2} was set to 0.2. SNPs that significantly
deviated from the Hardy-Weinberg expectation (\emph{p} \textless{}
0.001) were also filtered out. In total 25,172 unlinked autosomal SNPs
(\emph{r} \textsuperscript{2} \textless{} 0.2) with \textgreater{} 10×
coverage were used in the admixture analysis.

Secondly, we used \emph{GCTA} (Yang, Lee, Goddard, \& Visscher, 2011) to
carry out principal components analysis (PCA). The same set of 25,172
unlinked autosomal SNPs was analyzed. The score of the first component
(PC1; 1.33\% of variance) was plotted against the score of the second
component (PC2; 1.28\% of variance) of each individual to illustrate the
genetic affinity between the black-faced spoonbill samples.

\textbf{Historical demography}

\emph{Long-term historical demography}

To infer the historical demography of the black-faced spoonbill, we
applied \emph{SMC} ++(Terhorst et al., 2017) which employs a
coalescent-based method utilizing information on the site frequency
spectrum and linkage disequilibrium of a set of diploid genomes to infer
the historical demographic dynamics of the population concerned. To
enhance base calling accuracy, we used the VCF containing diploid
genotypes generated from all individuals (N=9) with coverage greater
than 15 for the \emph{SMC} ++ analysis. We used the command --CV to
perform cross-validation which yielded model parameters for calculating
the effective population size through time. The generation time of the
black-faced spoonbill was set to ten, and the mutation rate was set to
3.6× 10\textsuperscript{-8} per site per generation(Axelsson, Smith,
Sundström, Berlin, \& Ellegren, 2004). We also used \emph{SMC} ++ to
infer the historical demography of the royal spoonbill based on diploid
genotypes with coverage greater than 15 (N=7).

\emph{Detection of a recent bottleneck}

We used a model-selection approach to detect whether the black-faced
spoonbill had experienced a recent population bottleneck. Two
demographic models were constructed for the recent demography of the
black-faced spoonbill. The null model assumed that the effective
population size (\emph{N} \textsubscript{e}) of the black-faced
spoonbill had been constant (constant size model). The alternative model
assumed that the black-faced spoonbill had experienced a recent
bottleneck event (recent bottleneck model) and therefore contained the
following five parameters: \emph{N} \textsubscript{e} before the
bottleneck,\emph{N} \textsubscript{e} for the bottleneck, \emph{N}
\textsubscript{e}after bottleneck, the date of the initiation of the
bottleneck and the date of the end of the bottleneck. We used a method
implemented in\emph{Fastsimcoal2} (Excoffier et al., 2013) to conduct
model selection and parameter inference for the selected model. First,
we generated the observed folded (minor allele) site frequency spectrum
from the 215,722 unlinked autosomal SNPs of the black-faced spoonbill
with the script\emph{easySFS}
(https://github.com/isaacovercast/easySFS\#easysfs). For the constant
population size model, the prior for the\emph{N} \textsubscript{e} was
set to 100-1,000 haploid individuals. Assuming a generation time of ten
years and that the bottleneck event occurred between the 1950 and the
1980s(La Touche, 1931)\textsuperscript{,}(Austin, 1948), the priors for
the dates of bottleneck initiation (\emph{T} \textsubscript{bot}) and
termination (\emph{T} \textsubscript{endbot}) were set to 2-5 and 2-10
generations ago, respectively. Bracketing the long-term \emph{N}
\textsubscript{e}(1/2 to 2×) since the last glacial maximum estimated
from the\emph{SMC} ++, we set the prior of \emph{N} \textsubscript{e}
for the pre-bottleneck (\emph{N} \textsubscript{anc}) to 7,500-30,000
haploid individuals. Considering black-faced spoonbill's population size
of 288- 4000 since 1988, the priors for bottleneck (\emph{N}
\textsubscript{bot}) and post-bottleneck \emph{N}
\textsubscript{e}(\emph{N} \textsubscript{cur}) were set to 2-100 and
100-1,000 haploid individuals respectively. Using uniform random samples
from the priors of the two demographic models, we generated 100,000
folded SFSs for each model from the coalescent simulations and ran 40
optimization cycles to estimate each parameter and its composite
likelihood in each simulation. The set of parameters with the highest
likelihood was used for model selection. The simulation procedure
described above was repeated 1,000 times, and a total of 5,000 sets of
parameters were generated for each model. The parameters of the model
from the set of simulations with the highest estimated likelihood was
chosen as the best estimate of parameters for a given model. The Akaike
information criterion (AIC) was used to compare the best simulation of
the two demographic models to the observed folded site frequency
spectrum (SFS). We used \selectlanguage{greek}\emph{Δ} \selectlanguage{english}AIC and AIC weight (\emph{w} ) to
evaluate which model better fitted the observed folded SFS. Then we ran
the parameter estimation procedure under the best model to obtain 100
bootstrapped folded SFSs. Based on these bootstrapped folded SFS, we ran
the parameter estimate procedure for each bootstrapped folded SFS again
to compute the confidence interval of each parameter of the selected
model.

\textbf{Genetic consequences of the recent bottleneck}

\emph{Inbreeding}

In a small population random mating can result in a higher level of
inbreeding. To assess the level of inbreeding, we used \emph{PLINK} to
calculate the number and length (\textgreater{} 100 Kb) of ROHs through
the autosomal genome for each black-faced and royal spoonbill
individual.

We also calculated the genomic inbreeding coefficient,\emph{F}
\textsubscript{ROH}, which measures the proportion of ROHs in the
autosomal genome(M. C. Keller et al., 2011). To calculate\emph{F}
\textsubscript{ROH}, we divided the total length of ROH by the effective
autosomal genome size (the length of non-N bases of the non-sex
chromosomes).

In theory, a species with a higher inbreeding rate is expected to have a
lower linkage disequilibrium (LD) decay rate(Slatkin, 2008). Therefore
we also used \emph{PopLDdecay} v. 3.40(Zhang, Dong, Xu, He, \& Yang,
2019) to calculate the decay of LD, \emph{r} \textsuperscript{2}, for
the two spoonbill species. To calculate the pair-wise\emph{r}
\textsuperscript{2} between two SNPs, the maximum distance between them
was set to 150 Kb.

\emph{Genetic drift}

Genetic drift in a small population tends to cause the loss of rare
alleles, such as singletons(Motoo Kimura \& Ohta, 1969). To evaluate the
relative level of genetic drift in the black-faced spoonbill, we
compared the folded site frequency spectrum of the unlinked autosomal
SNPs between the two spoonbill species. We used the same folded SFS
for\emph{FastSimcoal2} analysis for the black-faced spoonbill. For the
royal spoonbill, we used the same filtering criteria as those for
evaluating population structure to generate a set of 232,409 unlinked
autosomal SNPs. Because the low frequency of rare alleles rendered by
genetic drift could lead to a higher Tajima's \emph{D} value(Tajima,
1989), we also used \emph{VCFtools} to calculate the Tajima's
\emph{D}value for all autosomal SNPs of both spoonbill species.

\emph{Relaxation of purifying selection}

Genetic drift can lead to relaxation of selection by lowering selection
efficacy. To evaluate the level of selection relaxation, we used the
descriptive model in \emph{RELAX} (Wertheim et al., 2015) to infer a
relaxation parameter \emph{k} for 7,835 orthologous genes between the
two spoonbill species and an outgroup, the crested ibis \emph{Nipponia
nippon} (S. Li et al., 2014). A relaxation parameter
\emph{k}\textgreater{} 1 indicates that the gene was under intensified
selection, and \emph{k~} \textless{} 1 implies relaxed selection. We
counted the number of genes with a \emph{k} \textless{} 1 for each
spoonbill species. First, we used the program \emph{Inparanoid} v.
4.1(Remm, Storm, \& Sonnhammer, 2001) to identify homologous genes
between the crested ibis and the black-faced spoonbill. Then, we
used\emph{gffread} v. 0.11.4 (https://github.com/gpertea/gffread) to
extract the coding region sequences of each homologous gene from each
black-faced and royal spoonbill individual. We used
\emph{MACSE}v0.9b1(Ranwez, Harispe, Delsuc, \& Douzery, 2011) to align
the homologous sequences from all three species. We removed sequences
that contained frameshifts and/or stop codons flagged by MACSE from the
subsequent analysis. Then we used \emph{RAxMAL8} v. 8.2.12(Stamatakis,
2014) to draw the maximum likelihood gene tree of all 5,657 well-aligned
homologous genes with the GTRGAMMA model. Gene trees of all orthologous
genes from all three species in focus were used for the
\emph{RELAX}analysis. The likelihood ratio test in \emph{RELAX} was used
to estimate\emph{k} for each homologous gene by comparing the model
which fixed\emph{k} =1 and the model which allowed \emph{k} to be
estimated. The MRCA branch and each branch of each spoonbill species
were alternatively set as the test branch in the \emph{RELAX} analysis.

\emph{Accumulation of deleterious mutations}

Relaxation of selection can result in an accumulation of deleterious
mutations in the genome. To evaluate whether the genome of the
black-faced spoonbill accumulated a significant number of deleterious
mutations, we counted the number of nonsynonymous substitutions in both
spoonbill species. First, we used \emph{ANNOVAR} (Wang, Li, \&
Hakonarson, 2010) to annotate substitutions in the coding region of both
spoonbill species to two categories, synonymous and non-synonymous
substitutions. Due to the lack of a suitable outgroup, we defined the
dominant allele in a species as the ancestral allele, and the minor
allele as the derived allele. Assuming that non-synonymous substitutions
are deleterious, we counted the number of non-synonymous substitutions
in each individual as accumulated deleterious mutations also known as
genetic load. To further quantify the deleterious effect of
non-synonymous substitutions, we used the Grantham score(Grantham, 1974)
to measure their potential phenotypic effect. The Grantham score
quantifies the differences in physiochemical similarity and evolutionary
exchangeability between two amino acids, where a higher score indicates
more physiochemical difference and lower evolutionary exchangeability
(i.e. higher deleterious phenotypic effect of nonsynonymous
substitution). Grantham scores can be classified into three
categories(Chun \& Fay, 2009), conservative changes (0-50), moderate
changes (51-100) and radical changes (\textgreater{}100). Assuming
conservative changes in amino acids to be neutral and different
deleterious non-synonymous substitutions to have an independent and
additive effect, we counted the accumulated deleterious effect of
non-synonymous substitutions by summing up the Grantham scores
categorized as moderate or radical changes for each individual.

\par\null

\textbf{References}

Alexander, D. H., \& Lange, K. (2011). Enhancements to the ADMIXTURE
algorithm for individual ancestry estimation. \emph{BMC Bioinformatics}
,\emph{12} (1), 246. doi: 10.1186/1471-2105-12-246

Allendorf, F. W. (1986). Genetic drift and the loss of alleles versus
heterozygosity. \emph{Zoo Biology} , \emph{5} (2), 181--190. doi:
10.1002/zoo.1430050212

Austin, O. L. Jr. (1948). The birds of Korea. \emph{Bull. Mus. Comp.
Zool.} , \emph{101} , 1--301.

Axelsson, E., Smith, N. G. C., Sundstr\selectlanguage{ngerman}öm, H., Berlin, S., \& Ellegren,
H. (2004). Male-Biased Mutation Rate and Divergence in Autosomal,
Z-Linked and W-Linked Introns of Chicken and Turkey. \emph{Molecular
Biology and Evolution} , \emph{21} (8), 1538--1547. doi:
10.1093/molbev/msh157

Baillie, J., Hilton-Taylor, C., \& Stuart, S. N. (2004). \emph{2004 IUCN
Red List of threatened species: A global species assessment} . IUCN.

Barrett, R. D. H., \& Schluter, D. (2008). Adaptation from standing
genetic variation. \emph{Trends in Ecology \& Evolution} , \emph{23}
(1), 38--44. doi: 10.1016/j.tree.2007.09.008

Bijlsma, R., \& Loeschcke, V. (2012). Genetic erosion impedes adaptive
responses to stressful environments. \emph{Evolutionary Applications}
,\emph{5} (2), 117--129. doi: 10.1111/j.1752-4571.2011.00214.x

BirdLife International. (2017). \emph{Platalea minor} . In \emph{The
IUCN Red List of Threatened Species 2017: E.T22697568A119347801} . IUCN.
Retrieved from
http://dx.doi.org/10.2305/IUCN.UK.2017-3.RLTS.T22697568A119347801.en.

BirdLife International 2000. (2000). \emph{Platalea minor} . In
\emph{The IUCN Red List of Threatened Species 2000:
E.T22697568A33426376.} IUCN.

Bolger, A. M., Lohse, M., \& Usadel, B. (2014). Trimmomatic: A flexible
trimmer for Illumina sequence data. \emph{Bioinformatics} ,\emph{30}
(15), 2114--2120. doi: 10.1093/bioinformatics/btu170

Bowkett, A. E. (2009). Recent Captive-Breeding Proposals and the Return
of the Ark Concept to Global Species Conservation. \emph{Conservation
Biology} , \emph{23} (3), 773--776. doi:
10.1111/j.1523-1739.2008.01157.x

Brondizio, E. S., Settele, J., Díaz, S., \& Ngo, H. T. (Series Eds.).
(2019). \emph{IPBES (2019): Global assessment report on biodiversity and
ecosystem services of the Intergovernmental Science-Policy Platform on
Biodiversity and Ecosystem Services} (p. 56). Bonn, Germany: IPBES
secetariat.

Charlesworth, D., \& Charlesworth, B. (1987). Inbreeding depression and
its evolutionary consequences. \emph{Annual Review of Ecology and
Systematics} , \emph{18} (1), 237--268.

Chesser, R. T., Yeung, C. K. L., Yao, C. T., Tian, X. H., \& Li, S.-H.
(2010). Molecular phylogeny of the spoonbills (Aves: Threskiornithidae)
based on mitochondrial DNA. \emph{Zootaxa} , (2603), 53--60.

Cho, Y. S., Hu, L., Hou, H., Lee, H., Xu, J., Kwon, S., \ldots{} Bhak,
J. (2013). The tiger genome and comparative analysis with lion and snow
leopard genomes. \emph{Nature Communications} , \emph{4} , 2433. doi:
10.1038/ncomms3433

Chong, J.-R., Pak, U.-I., Rim, C.-Y., \& Kim, T.-S. (1996). Breeding
biology of black-faced spoonbill \emph{Platalea minor} . \emph{Strix}
,\emph{14} , 1--10.

Chun, S., \& Fay, J. C. (2009). Identification of deleterious mutations
within three human genomes. \emph{Genome Research} , \emph{19} (9),
1553--1561. doi: 10.1101/gr.092619.109

Cooke, A. S. (1973). Shell thinning in avian eggs by environmental
pollutants. \emph{Environmental Pollution (1970)} , \emph{4} (2),
85--152. doi: 10.1016/0013-9327(73)90009-8

Crow, J. F., \& Kimura, M. (1970). \emph{An introduction to population
genetics theory.} New York: Harper \& Row, Publishers. Retrieved from
https://www.cabdirect.org/cabdirect/abstract/19710105376

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E.,
DePristo, M. A., \ldots{} Durbin, R. (2011). The variant call format and
VCFtools. \emph{Bioinformatics} , \emph{27} (15), 2156--2158. doi:
10.1093/bioinformatics/btr330

Ding, W., Lei, F., Yin, Z., \& Liu, R. (1999). The breeding sites of the
black-faced spoonbill \emph{Platalea minor} have been discovered in
Northern China. \emph{Bird Conservation International} , \emph{9} , 284.

Do, R., Balick, D., Li, H., Adzhubei, I., Sunyaev, S., \& Reich, D.
(2015). No evidence that selection has been less effective at removing
deleterious mutations in Europeans than in Africans. \emph{Nature
Genetics} , \emph{47} (2), 126--131. doi: 10.1038/ng.3186

Dobson, A. P., Bradshaw, A. D., \& Baker, A. J. M. (1997). Hopes for the
Future: Restoration Ecology and Conservation Biology. \emph{Science}
,\emph{277} (5325), 515--522. doi: 10.1126/science.277.5325.515

Ewen, J. G., Walker, L., Canessa, S., \& Groombridge, J. J. (2015).
Improving supplementary feeding in species
conservation.\emph{Conservation Biology} , \emph{29} (2), 341--349. doi:
10.1111/cobi.12410

Excoffier, L., Dupanloup, I., Huerta-Sánchez, E., Sousa, V. C., \& Foll,
M. (2013). Robust demographic inference from genomic and SNP
data.\emph{PLoS Genetics} , \emph{9} (10), e1003905--e1003905. doi:
10.1371/journal.pgen.1003905

Faber, R. A., \& Hickey, J. J. (1973). Eggshell thinning, chlorinated
hydrocarbons, and mercury in inland aquatic bird egg, 1969 and
1970.\emph{Pesticide Monitoring} , \emph{7} , 27--36.

Frankham, R. (1995). Effective population size/adult population size
ratios in wildlife: A review. \emph{Genetics Research} , \emph{66} (2),
95--107. doi: 10.1017/S0016672300034455

Frankham, R., Bradshaw, C. J. A., \& Brook, B. W. (2014). Genetics in
conservation management: Revised recommendations for the 50/500 rules,
Red List criteria and population viability analyses. \emph{Biological
Conservation} , \emph{170} , 56--63. doi: 10.1016/j.biocon.2013.12.036

Fu, W., Gittelman, R. M., Bamshad, M. J., \& Akey, J. M. (2014).
Characteristics of Neutral and Deleterious Protein-Coding Variation
among Individuals and Populations. \emph{The American Journal of Human
Genetics} , \emph{95} (4), 421--436. doi: 10.1016/j.ajhg.2014.09.006

Gattepaille, L. M., Jakobsson, M., \& Blum, M. G. (2013). Inferring
population size changes with sequence and SNP data: Lessons from human
bottlenecks. \emph{Heredity} , \emph{110} (5), 409--419. doi:
10.1038/hdy.2012.120

Gemmell, N. J., \& Akiyama, S. (1996). An efficient method for the
extraction of DNA from vertebrate tissues. \emph{Trends in Genetics:
TIG} , \emph{12} (9), 338--339. doi: 10.1016/s0168-9525(96)80005-9

Gnerre, S., MacCallum, I., Przybylski, D., Ribeiro, F. J., Burton, J.
N., Walker, B. J., \ldots{} Jaffe, D. B. (2011). High-quality draft
assemblies of mammalian genomes from massively parallel sequence
data.\emph{Proceedings of the National Academy of Sciences} , \emph{108}
(4), 1513--1518. doi: 10.1073/pnas.1017351108

Gómez\selectlanguage{english}-Sanchez, D., Olalde, I., Sastre, N., Ensenat, C., Carrasco, R.,
Marques\selectlanguage{english}-Bonet, T., \ldots{} Ramirez, O. (2018). On the path to
extinction: Inbreeding and admixture in a declining grey wolf
population. \emph{Molecular Ecology} , \emph{27} (18), 3599--3612. doi:
10.1111/mec.14824

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D.
A., Amit, I., \ldots{} Regev, A. (2011). Trinity: Reconstructing a
full-length transcriptome without a genome from RNA-Seq
data.\emph{Nature Biotechnology} , \emph{29} (7), 644--652. doi:
10.1038/nbt.1883

Grabherr, M. G., Russell, P., Meyer, M., Mauceli, E., Alfoldi, J., Di
Palma, F., \& Lindblad-Toh, K. (2010). Genome-wide synteny through
highly sensitive sequence alignment: Satsuma. \emph{Bioinformatics}
,\emph{26} (9), 1145--1151. doi: 10.1093/bioinformatics/btq102

Grantham, R. (1974). Amino acid difference formula to help explain
protein evolution. \emph{Science} , \emph{185} (4154), 862--864. Scopus.
doi: 10.1126/science.185.4154.862

Grier, J. W. (1982). Ban of DDT and subsequent recovery of Reproduction
in bald eagles. \emph{Science} , \emph{218} (4578), 1232--1235. doi:
10.1126/science.7146905

Grossen, C., Guillaume, F., Keller, L. F., \& Croll, D. (2020). Purging
of highly deleterious mutations through severe bottlenecks in Alpine
ibex. \emph{Nature Communications} , \emph{11} (1), 1--12. doi:
10.1038/s41467-020-14803-1

Haas, B. J., Delcher, A. L., Mount, S. M., Wortman, J. R., Smith, R. K.,
Hannick, L. I., \ldots{} White, O. (2003). Improving the Arabidopsis
genome annotation using maximal transcript alignment
assemblies.\emph{Nucleic Acids Research} , \emph{31} (19), 5654--5666.
doi: 10.1093/nar/gkg770

Hermisson, J., \& Pennings, P. S. (2005). Soft Sweeps: Molecular
Population Genetics of Adaptation From Standing Genetic
Variation.\emph{Genetics} , \emph{169} (4), 2335--2352. doi:
10.1534/genetics.104.036947

Hill, W. G., \& Robertson, A. (1966). The effect of linkage on limits to
artificial selection. \emph{Genetics Research} , \emph{8} (3), 269--294.
doi: 10.1017/S0016672300010156

Jones, F. C., Grabherr, M. G., Chan, Y. F., Russell, P., Mauceli, E.,
Johnson, J., \ldots{} Kingsley, D. M. (2012). The genomic basis of
adaptive evolution in threespine sticklebacks. \emph{Nature} ,\emph{484}
(7392), 55--61. doi: 10.1038/nature10944

Keller, L. F., \& Waller, D. M. (2002). Inbreeding effects in wild
populations. \emph{Trends in Ecology \& Evolution} , \emph{17} (5),
230--241. doi: 10.1016/S0169-5347(02)02489-8

Keller, M. C., Visscher, P. M., \& Goddard, M. E. (2011). Quantification
of Inbreeding Due to Distant Ancestors and Its Detection Using Dense
Single Nucleotide Polymorphism Data. \emph{Genetics} , \emph{189} (1),
237--249. doi: 10.1534/genetics.111.130922

Kennerley, K. J. (1990). A review of the status and distribution of
black-faced spoonbill. In V. Pickens (Ed.), \emph{The Hong Kong bird
report 1989} . Hong Kong: Hong Kong Bird Watching Society.

Kimura, M. (1962). On the probability of fixation of mutant genes in a
population. \emph{Genetics} , \emph{47} , 713--719.

Kimura, Motoo, Maruyama, T., \& Crow, J. F. (1963). The Mutation Load in
Small Populations. \emph{Genetics} , \emph{48} (10), 1303--1312.

Kimura, Motoo, \& Ohta, T. (1969). The Average Number of Generations
until Fixation of a Mutant Gene in a Finite Population. \emph{Genetics}
,\emph{61} (3), 763--771.

King, K. A., Meeker, D. L., \& Swineford, D. M. (1980). White\selectlanguage{english}-faced ibis
populations and pollutants in Texas, 1969--1976. \emph{Southwestern
Naturalist} , \emph{25} , 225--240.

Kirkpatrick, M., \& Jarne, P. (2000). The Effects of a Bottleneck on
Inbreeding Depression and the Genetic Load. \emph{The American
Naturalist} , \emph{155} (2), 154--167. doi: 10.1086/303312

La Touche, J. (1931). \emph{A handbook of birds of eastern China} .
London: Taylor \& Francis.

Lai, Y.-T., Yeung, C. K. L., Omland, K. E., Pang, E.-L., Hao, Y., Liao,
B.-Y., \ldots{} Li, S.-H. (2019). Standing genetic variation as the
predominant source for adaptation of a songbird. \emph{Proceedings of
the National Academy of Sciences} , \emph{116} (6), 2152--2157. doi:
10.1073/pnas.1813597116

Lande, R., \& Shannon, S. (1996). The role of genetic variation in
adaptation and population persistence in a changing
environment.\emph{Evolution; International Journal of Organic Evolution}
,\emph{50} (1), 434--437. doi: 10.2307/2410812

Larson, C. (2015). Hostile shores. \emph{Science} , \emph{350} (6257),
150--152. doi: 10.1126/science.350.6257.150

Li, H. (2011). A statistical framework for SNP calling, mutation
discovery, association mapping and population genetical parameter
estimation from sequencing data. \emph{Bioinformatics} , \emph{27} (21),
2987--2993. doi: 10.1093/bioinformatics/btr509

Li, H., \& Durbin, R. (2010). Fast and accurate long-read alignment with
Burrows--Wheeler transform. \emph{Bioinformatics} , \emph{26} (5),
589--595. doi: 10.1093/bioinformatics/btp698

Li, S., Li, B., Cheng, C., Xiong, Z., Liu, Q., Lai, J., \ldots{} Yan, J.
(2014). Genomic signatures of near-extinction and rebirth of the crested
ibis and other endangered bird species. \emph{Genome Biology} ,\emph{15}
(12), 557--557. doi: 10.1186/s13059-014-0557-1

Litvinenko, N. M., \& Shibaev, Y. V. (2005). Black-faced
Spoonbill\emph{Platalea minor} . In \emph{Red Data Book Primorsky krai.
Animals.}(p. 209). Vladivostok.

Lohmueller, K. E. (2014). The Impact of Population Demography and
Selection on the Genetic Architecture of Complex Traits. \emph{PLOS
Genetics} , \emph{10} (5), e1004379. doi: 10.1371/journal.pgen.1004379

Ma, Z., Melville, D. S., Liu, J., Chen, Y., Yang, H., Ren, W., \ldots{}
Li, B. (2014). Rethinking China's new great wall. \emph{Science}
,\emph{346} (6212), 912--914. doi: 10.1126/science.1257258

Martin, G., \& Lenormand, T. (2006). The Fitness Effect of Mutations
Across Environments: A Survey in Light of Fitness Landscape
Models.\emph{Evolution} , \emph{60} (12), 2413--2427. doi:
10.1111/j.0014-3820.2006.tb01878.x

Maruyama, T., \& Fuerst, P. A. (1985). Population bottlenecks and
nonequilibrium models in population genetics. II. Number of alleles in a
small population that was formed by a recent bottleneck.\emph{Genetics}
, \emph{111} (3), 675--689.

Mason, N. A., Pulgarin, P., Cadena, C. D., \& Lovette, I. J. (2020). De
Novo assembly of a high-quality reference genome for the horned lark
(\emph{Eremophila alpestris} ). \emph{G3
Genes\textbar{}Genomes\textbar{}Genetics} , \emph{10} (2), 475--478.
doi: 10.1534/g3.119.400846

Matheu, E., \& del Hoyo, J. (1992). Family Threskiornithidae (Ibises and
Spoonbills). In \emph{Handbook of the birds of theworld} (Vol. 1, pp.
472--507). Barcelona: Lynx Edicions.

Mays, H. L., Hung, C.-M., Shaner, P.-J., Denvir, J., Justice, M., Yang,
S.-F., \ldots{} Primerano, D. A. (2018). Genomic Analysis of Demographic
History and Ecological Niche Modeling in the Endangered Sumatran
Rhinoceros Dicerorhinus sumatrensis. \emph{Current Biology} ,\emph{28}
(1), 70-76.e4. doi: 10.1016/j.cub.2017.11.021

McQuillan, R., Leutenegger, A.-L., Abdel-Rahman, R., Franklin, C. S.,
Pericic, M., Barac-Lauc, L., \ldots{} Wilson, J. F. (2008). Runs of
Homozygosity in European Populations. \emph{The American Journal of
Human Genetics} , \emph{83} (3), 359--372. doi:
10.1016/j.ajhg.2008.08.007

Newton, I., \& Wyllie, I. (1992). Recovery of a Sparrowhawk Population
in Relation to Declining Pesticide Contamination. \emph{Journal of
Applied Ecology} , \emph{29} (2), 476--484. JSTOR. doi: 10.2307/2404515

Pimm, S. L., Jenkins, C. N., Abell, R., Brooks, T. M., Gittleman, J. L.,
Joppa, L. N., \ldots{} Sexton, J. O. (2014). The biodiversity of species
and their rates of extinction, distribution, and
protection.\emph{Science} , \emph{344} (6187). doi:
10.1126/science.1246752

Poole, A. F. (1989). \emph{Ospreys: A natural and unnatural history} .
Cambridge: Cambridge University Press.

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R.,
Bender, D., \ldots{} Sham, P. C. (2007). PLINK: A Tool Set for
Whole-Genome Association and Population-Based Linkage Analyses.\emph{The
American Journal of Human Genetics} , \emph{81} (3), 559--575. doi:
10.1086/519795

Ranwez, V., Harispe, S., Delsuc, F., \& Douzery, E. J. P. (2011). MACSE:
Multiple Alignment of Coding SEquences Accounting for Frameshifts and
Stop Codons. \emph{PLOS ONE} , \emph{6} (9), e22594. doi:
10.1371/journal.pone.0022594

Remm, M., Storm, C. E. V., \& Sonnhammer, E. L. L. (2001). Automatic
clustering of orthologs and in-paralogs from pairwise species
comparisons. \emph{Journal of Molecular Biology} , \emph{314} (5),
1041--1052. doi: 10.1006/jmbi.2000.5197

Ridgway, J., \& Shimmield, G. (2002). Estuaries as Repositories of
Historical Contamination and their Impact on Shelf Seas.\emph{Estuarine,
Coastal and Shelf Science} , \emph{55} (6), 903--928. doi:
10.1006/ecss.2002.1035

Risebrough, R. W. (1986). Pesticides and Bird Populations. In R.
Johnston (Ed.), \emph{Current Ornithology} (Vol. 3, pp. 397--427).
Springer.

Robertson, A., \& Waddington, C. H. (1960). A theory of limits in
artificial selection. \emph{Proceedings of the Royal Society of London.
Series B. Biological Sciences} , \emph{153} (951), 234--249. doi:
10.1098/rspb.1960.0099

Robinson, J. A., Ortega-Del Vecchyo, D., Fan, Z., Kim, B. Y., vonHoldt,
B. M., Marsden, C. D., \ldots{} Wayne, R. K. (2016). Genomic Flatlining
in the Endangered Island Fox. \emph{Current Biology} , \emph{26} (9),
1183--1189. doi: 10.1016/j.cub.2016.02.062

Rubenstein, D. R., Chamberlain, C. P., Holmes, R. T., Ayres, M. P.,
Waldbauer, J. R., Graves, G. R., \& Tuross, N. C. (2002). Linking
Breeding and Wintering Ranges of a Migratory Songbird Using Stable
Isotopes. \emph{Science} , \emph{295} (5557), 1062--1065. doi:
10.1126/science.1067124

Shibaev, Y. V. (2010). Breeding of the Black-faced Spoonbill
(\emph{Platalea minor} ) in Peter the Great Bay (Primorye, Russia).
Situation and prospects. In \emph{Annual Report of Pro Natura Fund}(Vol.
19, pp. 151--164).

Simons, Y. B., Turchin, M. C., Pritchard, J. K., \& Sella, G. (2014).
The deleterious mutation load is insensitive to recent population
history. \emph{Nature Genetics} , \emph{46} (3), 220--224. doi:
10.1038/ng.2896

Slatkin, M. (2008). Linkage disequilibrium---Understanding the
evolutionary past and mapping the medical future. \emph{Nature Reviews
Genetics} , \emph{9} (6), 477--485. doi: 10.1038/nrg2361

Smith, R. K., Pullin, A. S., Stewart, G. B., \& Sutherland, W. J.
(2010). Effectiveness of Predator Removal for Enhancing Bird
Populations. \emph{Conservation Biology} , \emph{24} (3), 820--829. doi:
10.1111/j.1523-1739.2009.01421.x

Stamatakis, A. (2014). RAxML version 8: A tool for phylogenetic analysis
and post-analysis of large phylogenies. \emph{Bioinformatics} ,\emph{30}
(9), 1312--1313. doi: 10.1093/bioinformatics/btu033

Stanke, M., \& Morgenstern, B. (2005). AUGUSTUS: A web server for gene
prediction in eukaryotes that allows user-defined
constraints.\emph{Nucleic Acids Research} , \emph{33} (Web Server
issue), W465--W467. doi: 10.1093/nar/gki458

Sung, Y.-H., Tse, I. W.-L., \& Yu, Y.-T. (2018). Population trends of
the Black-faced Spoonbill \emph{Platalea minor} : Analysis of data from
international synchronised censuses. \emph{Bird Conservation
International} , \emph{28} (1), 157--167. doi: 10.1017/S0959270917000016

Szafraniec, K., Borts, R. H., \& Korona, R. (2001). Environmental stress
and mutational load in diploid strains of the yeast \emph{Saccharomyces
cerevisiae} . \emph{Proceedings of the National Academy of Sciences of
the United States of America} , \emph{98} (3), 1107--1112.

Tajima, F. (1989). The effect of change in population size on DNA
polymorphism. \emph{Genetics} , \emph{123} (3), 597--601.

Terhorst, J., Kamm, J. A., \& Song, Y. S. (2017). Robust and scalable
inference of population history from hundreds of unphased whole
genomes.\emph{Nature Genetics} , \emph{49} (2), 303--309. doi:
10.1038/ng.3748

Ueta, M., Melville, D. S., Wang, Y., Ozaki, K., Kanai, Y., Leader, P.
J., \ldots{} Kuo, C.-Y. (2002). Discovery of the breeding sites and
migration routes of Black-faced Spoonbills \emph{Platalea minor}
.\emph{Ibis} , \emph{144} (2), 340--343. doi:
10.1046/j.1474-919X.2002.00037.x

Wang, K., Li, M., \& Hakonarson, H. (2010). ANNOVAR: Functional
annotation of genetic variants from high-throughput sequencing
data.\emph{Nucleic Acids Research} , \emph{38} (16), e164--e164. doi:
10.1093/nar/gkq603

Warren, W. C., Hillier, L. W., Tomlinson, C., Minx, P., Kremitzki, M.,
Graves, T., \ldots{} Cheng, H. H. (2017). A New Chicken Genome Assembly
Provides Insight into Avian Genome Structure. \emph{G3: Genes, Genomes,
Genetics} , \emph{7} (1), 109--117. doi: 10.1534/g3.116.035923

Waterhouse, R. M., Seppey, M., Simao, F. A., Manni, M., Ioannidis, P.,
Klioutchnikov, G., \ldots{} Zdobnov, E. M. (2018). BUSCO Applications
from Quality Assessments to Gene Prediction and
Phylogenomics.\emph{Molecular Biology and Evolution} , \emph{35} (3),
543--548. doi: 10.1093/molbev/msx319

Weber, K. (2004). Population Size and Long-term Selection. In J. Janick
(Series Ed.), \emph{Plant Breeding Reviews, Part1: Long-term Selection:
Maize} (Vol. 24, pp. 249--268). Hoboken, New Jersey: John Wiley \& Sons.

Wertheim, J. O., Murrell, B., Smith, M. D., Kosakovsky Pond, S. L., \&
Scheffler, K. (2015). RELAX: Detecting Relaxed Selection in a
Phylogenetic Framework. \emph{Molecular Biology and Evolution}
,\emph{32} (3), 820--832. doi: 10.1093/molbev/msu400

Wong, M. H., Leung, A. O. W., Chan, J. K. Y., \& Choi, M. P. K. (2005).
A review on the usage of POP pesticides in China, with emphasis on DDT
loadings in human milk. \emph{Chemosphere} , \emph{60} (6), 740--752.
doi: 10.1016/j.chemosphere.2005.04.028

Yang, J., Lee, S. H., Goddard, M. E., \& Visscher, P. M. (2011). GCTA: A
Tool for Genome-wide Complex Trait Analysis. \emph{American Journal of
Human Genetics} , \emph{88} (1), 76--82. doi: 10.1016/j.ajhg.2010.11.011

Yeung, C. K. L., Tsai, P.-W., Chesser, R. T., Lin, R.-C., Yao, C.-T.,
Tian, X.-H., \& Li, S.-H. (2011). Testing Founder Effect Speciation:
Divergence Population Genetics of the Spoonbills \emph{Platalea
regia}and \emph{Pl. minor} (Threskiornithidae, Aves). \emph{Molecular
Biology and Evolution} , \emph{28} (1), 473--482. doi:
10.1093/molbev/msq210

Yeung, C. K.-L., Yao, C.-T., Hsu, Y.-C., Wang, J.-P., \& Li, S.-H.
(2006). Assessment of the historical population size of an endangered
bird, the black-faced spoonbill (\emph{Platalea minor} ) by analysis of
mitochondrial DNA diversity. \emph{Animal Conservation} , \emph{9} (1),
1--10. doi: 10.1111/j.1469-1795.2005.00007.x

Yu, Y. T., Li, C. H., Tse, I. W. L., \& Fong, H. H. N. (2020). \emph{The
International Black-faced Spoonbill Census 2020} . Hong Kong:
Black-faced Spoonbill Research Group, The Hong Kong Bird Watching
Society.

Zhang, C., Dong, S.-S., Xu, J.-Y., He, W.-M., \& Yang, T.-L. (2019).
PopLDdecay: A fast and effective tool for linkage disequilibrium decay
analysis based on variant call format files. \emph{Bioinformatics
(Oxford, England)} , \emph{35} (10), 1786--1788. doi:
10.1093/bioinformatics/bty875

Zhao, S., Zheng, P., Dong, S., Zhan, X., Wu, Q., Guo, X., \ldots{} Wei,
F. (2013). Whole-genome sequencing of giant pandas provides insights
into demographic history and local adaptation. \emph{Nature Genetics}
,\emph{45} (1), 67--71. doi: 10.1038/ng.2494

\textbf{Data Accessibility:}

DNA sequences: NCBI BioProject ID: PRJNA597661



\textbf{Figure legends}

Figure 1. A. Distribution map of the black-faced (red: breeding range;
gold: wintering range) and royal spoonbills (blue; modified from Hancock
et al. 1992). The arrows indicate the two sampling sites, Tainan, Taiwan
(blue) and Hong Kong (orange), of the black-faced spoonbill used in
current study. B. The cross-validation (CV) error of the\emph{ADMIXTURE}
analysis (\emph{K} : 1- 5) of the black-faced spoonbill from Taiwan and
Hong Kong was smallest when \emph{K} was equal to 1. The origin of
samples is labeled at the bottom: T: Tainan, Taiwan; HK: Hong Kong. C.
Samples of the black-faced spoonbill from Tainan, Taiwan and Hong Kong
did not form distinct clusters in the plot of scores of PC1 (1.33\% of
genetic variance) versus PC2 (1.28\% of genetic variance). The origin of
the samples is indicated by colors: blue: Tainan, Taiwan; orange: Hong
Kong.

Figure 2. A. The historical effective population size (\emph{N}
\textsubscript{e}) trajectories estimated by \emph{SMC} ++ for both the
black-faced and the royal spoonbill since the last glacial maximum (c.a.
22,000 years before the present); the full\emph{N} \textsubscript{e}
trajectories of the two spoonbill species estimated by \emph{SMC} ++ are
shown within the inset. The black-faced spoonbill in magenta and the
royal spoonbill in apricot. B. Demographic parameters of the recent
bottleneck model for the black-faced spoonbill estimated from
simulations of \emph{Fastsimcoal2} based on the folded site frequency
spectrum of the 215,722 unlinked autosomal SNPs. The solid line is the
maximum likelihood estimate of various\emph{N} \textsubscript{e}; The
dashed line is the 95\% CI of\emph{N} \textsubscript{anc} and \emph{N}
\textsubscript{cur}.\emph{T} \textsubscript{bot}: initiation time of the
bottleneck;\emph{T} \textsubscript{endbot}: termination time of the
bottleneck;\emph{N} \textsubscript{anc}: long term \emph{N}
\textsubscript{e} before the bottleneck; \emph{N} \textsubscript{bot}:
\emph{N} \textsubscript{e}during the bottleneck; and \emph{N}
\textsubscript{cur}: the post-bottleneck \emph{N} \textsubscript{e}.

Figure 3. A. The black-faced spoonbills sequenced had more runs of
homozygosity (\textgreater{} 100 Kb; NROH) and longer total runs of
homozygosity (SROH) than the royal spoonbills sequenced. The black-faced
spoonbills are shown in circles, and the royal spoonbills in triangles.
B. The linkage (\emph{r} \textsuperscript{2}) between SNPs on the same
scaffold decayed more slowly in the black-faced spoonbill genome than in
the royal spoonbill genome. C. The proportion of singletons in the
folded site frequency spectrum for unlinked (\emph{r}
\textsuperscript{2}\textless{} 0.2) autosomal SNPs was significantly
lower in the black-faced spoonbill genomes (21.9\%) than in the royal
spoonbill genome (31.7\%; \selectlanguage{greek}\emph{χ} \selectlanguage{english}\textsuperscript{2}= 24.82, \emph{p}
= 0.000). D. The Tajima's \emph{D} value for all autosomal SNPs of the
black-faced spoonbill (0.24\selectlanguage{ngerman}± 0.98, \emph{n} = 4,522,032) is higher than
that of the royal spoonbill (0.05± 0.99,\emph{n} = 4,437,044; \emph{t} =
296.51, \emph{p} \textless{} 0.001). The black-faced spoonbill is shown
in magenta and the royal spoonbill in apricot.

Figure 4. A. Significantly more genes in the black-faced spoonbill
(2,094, 37.0\%) are under relaxed selection (relaxation parameter
\emph{k} \textless{}1,\emph{p} \textless{} 0.05) than in the royal
spoonbill (1,542, 27.2\%;\selectlanguage{greek}\emph{χ} \selectlanguage{english}\textsuperscript{2}= 123.04, \emph{p}
= 2.2e-16; Figure 4A); but not for genes under intensified selection
\selectlanguage{greek}\emph{χ} \selectlanguage{english}\textsuperscript{2}= 3.17, \emph{p} = 0.075). B. The mean
number of nonsynonymous substitutions found in each individual of the
black-faced spoonbill genome (8294.00\selectlanguage{ngerman}± 332.66 substitutions) is higher
than that found in the royal spoonbill genome (8025.78± 211.77
substitutions; \emph{t} = -1.91,\emph{p} \textless{} 0.035). C. The
average Grantham's score for each deleterious nonsynonymous substitution
(Grantham's score \textgreater{} 50) is significantly higher for the
black-faced spoonbill (56.98± 0.69) than for the royal spoonbill (53.25±
1.83; two-sample \emph{t} test= 6.5154, \emph{p} = 3.6e-6). D. The ratio
of homozygous to heterozygous nonsynonymous substitutions in the genome
of the black-faced spoonbill (0.13± 0.04) is higher than that in the
royal spoonbill genome (0.11± 0.007; \emph{t} = -1.86, \emph{p}
\textless{} 0.04). The black-faced spoonbill is shown in magenta and the
royal spoonbill in apricot.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig1-sp/Fig1-sp}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig2-sp/Fig2-sp}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig3-sp/Fig3-sp}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig4-sp/Fig4-sp}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

