Whole-genome analysis reveals the diversification of Galapagos rail (Aves: Rallidae) and confirms the success of goat eradication programs

Abstract

Similar to other insular birds around the world, the Galapagos rail (Laterallus spilonota Gould, 1841) exhibits reduced flight capacity following its colonization of the archipelago ~1.2 mya. Despite their short evolutionary history, rails have colonized seven different islands spanning the entire width of the archipelago. Galapagos rails were once common on islands with sufficiently high altitudes to support shrubs in humid habitats. After humans introduced goats, this habitat was severely reduced due to overgrazing. Habitat loss devastated some rail populations, with less than 50 individuals surviving, rendering the genetic diversity of Galapagos rail a pressing conservation concern. Additionally, one enigma is the reappearance of rails on the island of Pinta after they were considered extirpated. Our approach was to investigate the evolutionary history and geographic distribution of Galapagos rails as well as examine the genome-wide effects of historical population bottlenecks using 39 whole genomes across different island populations. We recovered an early divergence of rail ancestors leading to the isolated populations on Pinta and a second clade comprising the rest of the islands, historically forming a single landmass. Subsequently, the separation of the landmass ~900 kya may have led to the isolation of the Isabela population with more panmictic populations found on Santa Cruz and Santiago islands. We found that rails genomes contain long runs of homozygosity (>2 Mb) that could be related to the introduction of goats. Finally, our findings show that the modern eradication of goats was critical to avoiding episodes of inbreeding in most populations.

Graphical Abstract

Open in new tabDownload slide

Introduction

The Galapagos archipelago is an ideal system for studying patterns of evolution due to its well-defined geological history, the presence of several island replicates, and the marked isolation from mainland South America (Black 1974; Holland and Hadfield 2004; Grant and Grant 2008; Poulakakis et al. 2012; Garcia-R and Matzke 2021). While dispersion among islands is common, vicariance is also possible via the subdivision of widespread ancestral populations caused by an environmental barrier (Caccone et al. 2002; Arbogast et al. 2006; Parent and Crespi 2006; Steinfartz et al. 2009; Poulakakis et al. 2012, 2020). Specifically, less than 1 million years ago (mya), several islands in the central region of the archipelago were connected, constituting a single landmass (Geist et al. 2014; Karnauskas et al. 2017; Schwartz et al. 2018). Rising sea levels over time resulted in the isolation of several of these islands from one another and is suggested that this triggered the diversification of ancestors into current lineages of giant tortoises, mockingbirds, vermilion flycatchers, and others (Arbogast et al. 2006; Poulakakis et al. 2012; Carmi et al. 2016).

Among the avifauna that inhabit the archipelago, the Galapagos rail (Laterallus spilonota) represents an excellent model species for investigating the role of island diversification following the divergence of this landmass. First, the species originated in the Galapagos from a single volant ancestor that arrived from mainland South America ~1.2 mya (Chaves et al. 2020). Second, the presence of rails on several islands suggests a dynamic dispersal of its ancestral population across the archipelago (Poulakakis et al. 2012; Geist et al. 2014; Karnauskas et al. 2017; Schwartz et al. 2018). In particular, the breakup of the previous central landmass of the archipelago might have influenced the pattern of diversification of the Galapagos rail (Geist et al. 2014; Karnauskas et al. 2017; Schwartz et al. 2018). Rails are well documented to evolve flightlessness on islands (Slikas et al. 2002; Wright et al. 2016; Garcia-R and Matzke 2021). In the Galapagos islands, rails are not able to fly more than 30 m (Franklin et al. 1979). This may have facilitated further genetic differentiation and diversification (although Galapagos Rails are not believed to be fully flightless).

A major concern in isolated oceanic islands is the impact of invasive species on native biodiversity. Numerous flightless birds around the world, including rails, have gone extinct due to species introductions by humans (Olson and James 1982; Diamond 1984; Steadman 1995; Duncan and Blackburn 2004). The Galapagos rail, a species listed as Vulnerable in the IUCN Red List (Freile et al. 2019), was once common on islands characterized by humid environments and dense ground vegetation at high elevations (Franklin et al. 1979; Rosenberg 1990; Donlan et al. 2007). When goats were first introduced to the Galapagos islands in 1959 (Hoeck 1984), their population rapidly grew to at least 20,000 individuals by 1970 (Weber 1971). Consequently, goats drastically reduced the highland vegetation through overgrazing, leading rail populations from San Cristobal and Floreana islands to become extirpated (Hamann 1979; Rosenberg 1990; Gibbs et al. 2003; Campbell et al. 2004; Chaves et al. 2020). Overgrazing destroyed the vegetation that is used by rails for nesting and feeding (Franklin et al. 1979). Specifically, rails used herbaceous steam to build semi-domed-shaped nests covered by dense low vegetation like ferns typical of humid environments. Also, the main food resource (i.e. insects) of rails is found in vegetation that grows in the highlands. “Populations from Isabela and Santiago were highly impacted, with less than 25 individuals estimated to have survived on each island (Rosenberg 1990; Donlan et al. 2007). On Pinta, rails disappeared when goats peaked in population size around 1970 (Franklin et al. 1979). Only after the implementation of goat eradication programs starting in the 1970s did the population size of rails begin to increase (Gibbs et al. 2003; Donlan et al. 2007). The Santiago population increased from 23 to 233 individuals within a span of 18 yrs (Donlan et al. 2007), and rails on Pinta recovered following goat eradication (Franklin et al. 1979; Donlan et al. 2007). However, rail populations could still be at risk due to the erosion of genetic diversity and accumulation of harmful genomic mutations during periods of small population size (Lynch and Gabriel 1990; Lynch et al. 1995; Charlesworth and Willis 2009). It has been shown that such deleterious variants could have detrimental effects on the fitness of individuals, affecting reproductive success or compromising the immune system in wild populations (Wayne et al. 1991; Charlesworth and Willis 2009; Johnson et al. 2010; Dobrynin et al. 2015). Therefore, inbreeding remains a conservation concern, as it could lead to population decline, from which the species may never recover (Peterson et al. 2014; Kyriazis et al. 2021; Robinson et al. 2023).

Another major puzzle is the reappearance of rails on the island of Pinta after they were extirpated in 1970 and rediscovered in 1973 (Franklin et al. 1979; Rosenberg 1990; Chaves et al. 2020). It is uncertain whether the current population represents a recolonization event from a nearby island or a relic population that survived undetected. Understanding the reappearance of rails on this island is critical for conservation management (Chaves et al. 2020), especially the effect of possible inbreeding depression after a sharp population bottleneck.

Here, we used whole-genome data from several rail populations across islands to investigate: 1) the biogeographic history of diversification of rails in the Galapagos Islands; 2) the effect of island isolation and goat introduction on genic diversity and genetic differentiation among rails; 3) the potential origins of the reappearing rails on the island of Pinta. To accomplish these goals, we analyzed phylogenetic relationships, patterns of genome-wide diversity, and deleterious variation by sequencing 39 whole genomes of endemic Galapagos rails.

Methods

Sequencing mapping and genotype calling

We extracted genomic DNA from 39 samples of Galapagos rails from four different islands (Isabela, Santiago, Santa Cruz, and Pinta) using the Qiagen DNEasy kit (Qiagen, USA). Samples were selected based on sufficient quality and quantity of genomic DNA with a DNA fluorometer (Qubit 2.0), a NanoDrop spectrophotometer (ThermoFisher, USA), and gel electrophoresis. Samples were collected in 2020 on three different islands (Pinta, Isabela, Santa Cruz, and Santiago) on specific localities within each of the islands (Supplementary Table S1). Details concerning the sampling process of Galapagos rails, handling of specimens as well as necessary permits are reported in Chaves and collaborators (2020) under the “Material and Methods” as well as the “Acknowledgments” sections.

We generated short insert (150 bp) libraries from these extractions using the TruSeq Nano PCR-Free kit (Illumina, USA) and then conducted whole-genome sequencing on these libraries using the Novaseq 6000 sequencer (Illumina, USA) with 150 bp paired-end reads at the Vincent J. Coates Genomics Sequencing Laboratory, University of California, Berkeley and the Fulgent Genetics laboratory at El Monte, California.

Raw sequencing reads were processed and filtered using a modified pipeline from the Genome Analysis Toolkit GATK v.3.8 best practice guide (Van der Auwera et al. 2013). We used minimap2 with options “ax” and “sr” (Li 2021) to map the reads of the Galapagos rail to the sex chromosomes of the chicken (GCF_016699485.2). We excluded scaffolds JAKCOX010000024.1 and JAKCOX010000010.1 which were identified as the sexual chromosomes W and Z, respectively. We then mapped reads with good quality to the Black rail (Laterallus jamaicensis) genome assembly (GCA_022605575.1; Hall et al. 2023), using BWA-MEM (Li and Durbin 2009). The black rail has an estimated divergence time of 1.2 mya with respect to the Galapagos rail (Chaves et al. 2020). This reference assembly had a mean coverage of 37.7×, a total length of 1.4 Gbp, and a contig and scaffold N50/L50 of 7.4 Mbp/44. We used GATK HaplotypeCaller to conduct joint genotype calling from sites of the Galapagos rail genomes that were mapped to the black rail reference genome. Then, we filtered the called genotypes for coverage and quality. We kept only genotypes that had a minimum of five reads at a given position and a high-quality score (Phred scores ≥ 20), and no more than the 99th percentile of coverage for each sample. Other variant filtering criteria followed GATK v.3.8 best practices guide (Van der Auwera et al. 2013). We filtered out CpG islands, indels, multi-nucleotide polymorphisms, and sites with more than one alternate allele. The command lines used for read mapping, variant calling, and filtering are available at https://github.com/dechavezv/2nd.paper.v2.

Phylogenetic analysis

We reconstructed a species tree representing the relationships among 39 Galapagos rail genomes using ASTRAL-III (Zhang et al. 2018). First, we extracted 6,540 alignments, each 25 kb in length, from the entire set of 35 autosomal chromosomes. Each independent window was aligned with PRANK v.150803 using iteration (-F once option). After trimming each of the multiple alignments with Gblocks (Castresana 2000), we calculated the maximum-likelihood gene tree phylogeny of the 6,540 alignments using IQ-TREE 2 (Minh et al. 2020) under the GTR model. For each alignment, the best tree was selected from the IQ-TREE 2 output, while the 100 bootstrap trees were merged into a single file per locus. We then created a consensus file with these trees with no more than 10% of missing data. The branch lengths that were shorter than 1e−05 substitution per site were collapsed. Also, we collapsed clades with posterior probabilities support lower than one using the SqCL pipeline (phylogeny_prep_astrid_astral.py). We used the best tree and consensus file of bootstrap trees to investigate the discordance between gene trees and the species tree using ASTRAL-III v.5.5 (Zhang et al. 2018). This resulted in maximum likelihood support values that were used to choose the best multi-locus tree. We ran 100 bootstrap replicates on this tree and further scored it to calculate the quartet support values and posterior probabilities for each node. We compiled 318,782 4-fold degenerate sites from 12,158 single-copy coding orthologs across 39 Galapagos rail genomes. Then we built a supermatrix and further used the MCMCTree tool from the PAML 4.8 package (Yang 2007), with the topology obtained from the ASTRAL-III analysis, to calculate divergence times (see Chavez et al. 2022 for details). We used 1 million years as a calibration prior for the divergence of Galapagos rails based on previous studies (Chaves et al. 2020), and a generation time of 2 yrs. Configuration file for MCMCTree analyses in the PAML 4.8 package, including parameter settings for the clock model (global), substitution model (HKY85), birth–death process, gamma priors on the transition/transversion rate ratio (kappa_gamma), and shape parameter for variable rates among sites (alpha_gamma) and the Dirichlet-gamma prior for the mean substitution rate (rgene_gamma). We ran the MCMC for 2,200,000 iterations, sampling every 2nd iteration, and discarding the first 200,000 iterations as burn-in.

BioGeoBEARS

We investigated the biogeographic history of Galapagos rails using the R package BioGeoBEARS (Matzke 2013). This tool estimates the maximum-likelihood distribution of hypothetical ancestors (internal nodes) by modeling shifts between different geographical ranges along the phylogeny as a function of time. First, we tested three different models: dispersal-extinction-cladogenesis (DEC), dispersal vicariance analysis (DIVALIKE), and Bayesian analysis of biogeography (BAYAREALIKE). Additionally, we tested the same three models plus a founder effect parameter added to each model named “J”: DEC + J, DIVALIKE + J, and BAYAREALIKE + J (Matzke 2013). The estimation to determine the best-fitting model was conducted by comparing model likelihoods through the Akaike information criterion (AIC) scores and AIC weights (Supplementary Table S2). We used default parameters and set the max_range_size = 4. The scripts and files used to run these models can be found at https://github.com/dechavezv/2nd.paper.v2/blob/main/3-Phylogenomics/BioGeoBEARS/. This model estimates the most probable ancestral geographical range in a particular node in the tree. Multiple ranges can be illustrated in a particular node with their probabilities shown by different colors in pie charts.

F_ST, FastStructure, and PCA

We calculated F_ST for every island-pair population with SNPRelate (Zheng et al. 2012). We used the method “W&C84” on .gds files and removed sites that were monomorphic and with no more than 20% missing data. The geographic distances were obtained using the geosphere package in R (https://CRAN.R-project.org/package=geosphere). Then, we created a plot of F_ST vs. distance.

To estimate the level of relatedness among samples, we converted .vcf files to the gds inputs required by PLINK v.2.0 using SNPRelate v.3.18 (Zheng et al. 2012). This tool was used on a subset of ~10k high-quality SNPs pruned for LD, with a r² threshold of 0.2, and a maximum missing rate greater than 0.05. Then, we assessed the level of relatedness among rail samples from different islands by estimating identical-by-descent (IBD) using PLINK’s method of moments approach (snpgdsIBDMoM; Purcell et al. 2007), with a minor allele frequency cutoff of 0.05. When the resulting kinship scores between a pair were greater than 0.2, we randomly removed one sample from this pair. As a result, we removed four samples from Santa Cruz and one sample from Isabela.

To evaluate the existence of genetic structure among rail populations, we used SNPrelate to estimate the principal component analysis (PCA) on 10,274 SNPs pruned for LD (threshold of 0.2) and a minor allele frequency of 0.1. We used Windows 500 kb in length. Among the different principal components (PC) estimated, we plotted PC1 and PC2 with ggplot2 (Wickham 2016), as these components had most of the variation explained by the data. We further conducted a Bayesian clustering analysis using fastStructure (Raj et al. 2014). For each model, we tested 10 different partitions from k = 2 to k = 10 and selected the partition that had the optimal k value based on the “chooseK” criteria implemented in fastStructure as well as their marginal likelihoods.

Treemix

We used TreeMix (Pickrell and Pritchard 2012) to further evaluate the relationship among island populations and evaluate levels of admixture among Galapagos rails. Specifically, we used the whole-genome allele frequencies of 39 samples to discover the best population tree. This tree was obtained by comparing independent maximum-likelihood scores. TreeMix then compared the covariance of the depicted tree to the observed covariance between populations to identify possible admixture events (Pickrell and Pritchard 2012). Population pairs that fail to fit the modeled tree could be subject to admixture. We tested several models in our whole-genome data for the subset of rails without related individuals as identified by SNPrelate (Zheng et al. 2012). First, we tested a model with no migration and blocks of 500 SNPs (-k 500). Then, we added eight different migration edges (m = 1 to m = 10) to the phylogenetic model in TreeMix. The migration score is obtained by comparing the standard error of a model without gene flow (migration edges = 0) with a model that allows for admixture (migration edges > 0). A positive migration weight suggested that adding migration to the model reduced the stand error and improved the fit to the data.

To determine the best-fitting model, we used the -global option and -se option to calculate the standard errors of the residual covariance. Finally, we considered the residual covariance between different population pairs and divided this difference by the standard error from all different pair populations.

Genomic diversity

We examined the site heterozygosity in non-overlapping 100 kb windows across the genome of all 39 Galapagos rails. We defined heterozygosity as the number of heterozygous genotypes divided by the total number of sites that were called. The total genotypes called within each window included: the sum of heterozygous, homozygous-derived, and homozygous reference genotypes. We kept only windows with no more than 20% of missing data. The script used to calculate 100 kb windows heterozygosity was modified from Robinson et al. (2019) and is available at https://github.com/dechavezv/2nd.paper.v2/tree/main/4-Demography/Heterozygosity/WindowHet. We then quantified the extent of Runs of Homozygosity (ROH) in rails using PLINK v.2.0 (Purcell et al. 2007). The parameters chosen to calculate ROH were: SNPs within a window = 200, heterozygotes allowed within a window = 3, and missing sites within a window = 50. We binned these segments into three different size categories using PLINK v.2.0 (Purcell et al. 2007). The size categories were: short (0.5 MB <), medium ROH (0.5 to 2 MB), and long ROH (> 2 Mb).

Deleterious variation

To investigate the potential consequences of past population declines on the Galapagos rail we analyzed protein-coding variants. First, we mapped individual raw reads to the closet species with full annotation, the Inaccessible Island rail (Atlantisia rogersi) which is 7 million apart from the Galapagos rails (Chaves et al. 2020) and located at https://www.ncbi.nlm.nih.gov/assembly/GCA_013401215. Then we conducted a joint variant calling with the tool haplotype caller from GATK v.3.8 (Van der Auwera et al. 2013). We evaluated the effects of the joint variant call from each of the 39 Galapagos rails on their associated protein-coding genes using the tool SnpEff (Cingolani et al. 2012). We used the gff file of the Inaccessible Island rail with coordinates for 13,465 genes to annotate mutations as synonymous, nonsynonymous, and loss of function (LOF) (e.g. stop-gain or frameshift variant). When more than one transcript was available for a single gene, we used the longest transcript for further analysis. The script to annotate and identify deleterious variants can be found at: https://github.com/dechavezv/2nd.paper.v2/tree/main/4-Demography/DeltVariation.

Mitochondrial genomes

We mapped the 39 raw reads to the Galapagos rail mitochondrial genome assembly (MW067132) reported by Chaves (Chaves et al. 2020). From these mitogenomes, we represented population relationships using a median-joining haplotype network with the haploNet function in the R package Pegas v.1.3 (Paradis 2010). Additionally, we aligned the Galapagos rails mitogenomes as well as one black rail with the tool PRANK v.150803 (Loytynoja and Goldman 2008) using iteration (-F once option) and built a phylogenetic tree using RAxML v8 (Stamatakis 2014). For the outgroup, we used the mitochondrial genome from the black rail (CM040151.1) reported by Hall et al. (2023).

Demographic inference.

To infer rails, recent demographic model we used the program GONE (Santiago et al. 2020) to infer demographic history in the recent past from patterns of linkage disequilibrium (LD) for the four populations of Galapagos rails. We used 150 scaffolds of the reference genome representing the N90 (1.4 Mb) of the black rail reference genome. We ran the model with the following parameters based on recombination rate in birds cMMb = 0.14 (Singhal et al. 2015), a minor allele frequency of 0.01 with 200 repetitions. Finally, to infer NE trajectories earlier in rails evolutionary history we used the pairwise and multiple sequentially Markovian coalescent (PSMC’) model within the program MSMC (Schiffels and Durbin 2014) to calculate the instantaneous inverse coalescence rates (IICR) among rail lineages. We ran MSMC on 39 genomes independently. We scaled the IICR from the MSMC model by 2µ to use it as a proxy of the effective population size (Ne). We consider a generation time of 2 yrs (Bird et al. 2020) and a mutation rate of 4.6e−09 (Smeds et al. 2016).

Results

Biogeographic and demographic history

To provide an evolutionary framework for the biogeographic history of Galapagos rails, we reconstructed a phylogenetic hypothesis from 39 individuals representing all extant island populations (Supplementary Table S1). We first mapped reads from Galapagos rail genomes to the Black rail (L. jamaicensis) genome assembly (GCA_022605575.). For all Galapagos rail samples, we found that more than 97% of reads were successfully mapped to the black rail (Supplementary Table S1). We then extracted 6, 540 alignments of 25 kb windows from these mapped genomes. and used maximum likelihood to construct independent phylogenetic trees from each of these windows. To account for phylogenetic discordance among these trees, we generated a consensus phylogenetic tree with ASTRAL-III (Zhang et al. 2018). Our species tree and divergence time estimates suggest that the ancestral population to all rails colonized the Galapagos Islands around 1 mya. The population on Pinta island diverged from other populations shortly after this time, ca. 0.9 mya (Fig. 1a), and has remained isolated ever since as suggested by the absence of gene flow (see Fig. 2b–d). The Isabela population diverged around 0.81 mya and was sister to a third clade that included samples from Santiago and Santa Cruz. These two island populations diverged 0.89 kya (Fig. 1a).

To elucidate the biogeographical history of rails in the Galapagos Islands, we used the R package BioGeoBEARS (Matzke 2013). This tool estimates the maximum-likelihood distribution of hypothetical ancestors (internal nodes) by modeling shifts between different geographical ranges along the phylogeny. Among six different models tested in BioGeoBEARS (see the “Methods” section for details), the BAYAREALIKE fit best data with a corrected AIC 26% higher than the other five models (Supplementary Table S2). This model assumes episodes of geographical dispersion, contraction, and changes of an ancestral distribution without speciation. According to this model, rails had a relatively large geographical range that contracted into smaller isolated areas. Specifically, there is a 55% probability that Galapagos rails originally occurred across a geographical range that included Isabela, Santa Cruz, and Santiago (see pie charts in Fig. 1a). From this ancestral distribution, there was a 98% probability that a population dispersed to the geographically isolated island of Pinta (Fig. 1a and b). Following this dispersion event, the ancestral distribution including Isabela, Santa Cruz, and Santiago then split into two different areas. One population occurred in Santa Cruz and Santiago with 48% probability and the other area occurred in Isabela with 98% probability (Fig. 1a).

To evaluate changes in effective population size (Ne) we used the multiple sequentially Markovian coalescent (MSMC) method (Schiffels and Durbin 2014). We applied the PSMC’ model by running MSMC separately on individual genomes. Island populations resulted in similar demographic trajectories during the past 1 million years. Two samples per population were plotted in Fig. 1c. Specifically, we found a sharp population size decline starting at 1 mya followed by a possible recent recovery between 100 and 10 kya (Fig. 1c). We further inferred the recent demography history of 39 high-coverage samples from patterns of LD for Pinta, Isabela, Santa Cruz, and Santiago populations. These four populations experienced a population size increase at different times starting with Pinta around 600 yrs ago followed by Isabela around 200 yrs (Fig. 1d). Santa Cruz and Santiago islands showed a Ne trajectory that increased around 150 yrs ago. After this period, the Ne trajectories of all four rail populations remained stable until 30 yrs ago when they experienced a rapid population decline.

Population structure and admixture

We investigated the amount of genetic differentiation among Galapagos rail populations by calculating the mean pairwise F_ST between islands using 10,274 SNPs pruned for LD from an original set of 651,558 SNPs. We found that rails from Pinta were the most differentiated population relative to other islands, with a pairwise F_ST ranging from 0.108 with Santa Cruz to 0.115 with Isabela. Individuals from Isabela were the second-most differentiated population with an F_ST ranging from 0.035 with Santa Cruz to 0.115 with Santiago. Rail populations from Santiago and Santa Cruz islands had the lowest F_ST value of 0.001 We found some degree of association between F_ST and geographical distance between islands. Specifically, values of F_ST increased with distance between populations (Supplementary Fig. S1).

To assess the level of population structure among island populations, we used the 10,274 SNPs to conduct a PCA. To avoid including related individuals in the PCA analysis, we used an identity-by-descent method with PLINK (Purcell et al. 2007) and removed one individual from a pair that had a kinship coefficient >0.2 from Isabela and four related individuals from Santa Cruz (kinship scores >0.2; see the “Methods” section for details). We found that the Pinta population was distinguished from the other three populations (Santiago, Santa Cruz, and Isabela) along PC1. This component represented 7.7% of the variance, separating samples between the northern and southern islands (Fig. 2a). PC2 separated the Isabela population, located in the west region of the Galapagos, from islands in the east (Santiago, Santa Cruz, and Pinta) with 6% of the variance explained (Fig. 2a). Populations from Santiago and Santa Cruz cluster together on both PCs (Fig. 2a), indicating low levels of population structure between these islands. To better elucidate any possible pattern of genetic differentiation between Santiago and Santa Cruz, we excluded individuals from Pinta from the PCA, with similar outcomes, confirming the lack of population structure among these populations (Supplementary Fig. S2).

To assess genetic connectivity and admixture among island populations, we used TreeMix (Pickrell and Pritchard 2012). This approach incorporates whole genome allele frequencies to find the best population tree as well as to infer gene flow between different populations. Our TreeMix results showed a similar topology to the species tree in Fig. 1a, with Santiago and Santa Cruz rail populations being closely related and sister to the Isabela population (Fig. 2b). In contrast, the Pinta population was found to be the most divergent as evidenced by a relatively long branch in the tree. Notably, this branch reached a drift score of 0.05 which was five times greater than the score of other populations at ~0.01 (see branch length in Fig. 2b). Importantly, we found evidence of admixture in the direction from the Isabela to Santiago population with a migration weight value of 0.5 (Fig. 2b and Supplementary Fig. S3).

To further investigate the genetic clustering and admixture of Galapagos rails, we used the program fastStructure (Raj et al. 2014) using values of k from 1 to 10. The different k values tested had similar marginal likelihoods (Supplementary Table S3). At k = 2, we found individuals from Pinta separating from another group composed of individuals from Isabela, Santiago, and Santa Cruz (Fig. 2c). At k = 3 one individual from Santiago and another from Santa Cruz formed an additional cluster (Fig. 2c and Supplementary Fig. S3). The other individuals have a relatively high assignment (>20%) to a cluster that includes Santa Cruz, Santiago, and Isabela (see k = 3 in Fig. 2c). These genetic clusters were observed across the different k’s tested (from k = 4 to k = 10).

Finally, we used a total of 39 mitochondrial genomes to further elucidate relationships among Galapagos rails populations. The haplotype network for the populations of rails resulted in three main dominant haplotypes (see Fig. 1d). The XI haplotype belonged to individuals from the Pinta island while the remaining two haplotypes II and III showed a relatively high proportion of haplotype sharing among the islands of Isabela, Santa Cruz, and Santiago, which suggests common ancestry. Our mitogenome-based phylogeny showed a similar pattern with samples from Pinta recovered as a monophyletic clade. In contrast, samples from Isabela, Pinta, and Santiago were paraphyletic (Supplementary Fig. S1).

Genetic diversity

To understand how demographic history has shaped patterns of genetic variation among Galapagos rail populations, we examined heterozygosity in non-overlapping 100 kb windows across the genome of every sequenced individual (Fig. 3 and Supplementary Fig. S4). We observed that genomes from different islands had regions of high heterozygosity alternating with homozygous stretches (Fig. 3b). The level of genome-wide heterozygosity was similar among individuals from different islands and ranged from 0.7 (heterozygosity/kb) on Pinta to 0.9 on Santa Cruz (Fig. 3b and Supplementary Table S4).

To more precisely examine the regions across the genomes depleted of heterozygosity, we quantified the extent of runs of homozygosity (ROH) grouped into three different size categories using PLINK (Purcell et al. 2007). Long ROH (>2 Mb) are a likely consequence of recent close inbreeding due to sharp population declines, whereas ROH of a short (<0.5 Mb) and medium (0.5 to 2 Mb) length could reflect ancient population declines as there would be more time for recombination to break up long segments (Ceballos et al. 2018). ROHs were relatively high in medium and long categories in the Pinta population (Fig. 3a and c), consistent with Pintas’ lower levels of heterozygosity (~0.7 heterozygosity/kb). The average total length of short ROH in the genomes from Pinta was similar to other islands with 50.8 Mb, which represents 3% of the genome, whereas the average total length of medium ROH was 90 Mb which represents 6% of the genome (Fig. 3a and c; Supplementary Table S4). Similarly, the total length of long ROH fragments was ~105 Mb on average (Supplementary Table S4), which is two times longer than the ROH of some rails from other islands and represents 8% of the genome (Fig. 3a and c). Individuals with relatively higher heterozygosity (~0.8 heterozygosity/kb) included samples from Isabela, Santa Cruz, and Santiago (Fig. 3a). These island populations comprise two categories based on ROH lengths: The first group includes samples with a few long ROH (>2 Mb) ranging from 6 Mb in genomes from Santa Cruz to 47 Mb in Santiago genomes, which represented less than 3% of the genome (Fig. 3a and c; Supplementary Table S4). The second group includes individuals with considerably long ROH (2 to 10 Mb; see red-colored labels in Fig. 3a) that represented up to 10% of the genome with values ranging from 53 Mb in Santa Cruz genomes to 140 Mb in Santiago genomes (Fig. 3a and c; Supplementary Table S4). Overall, our results suggest that recent population declines and inbreeding may have been pronounced in island populations, especially in Pinta.

Deleterious variation

To investigate the effects of demography (e.g. population declines) on deleterious variation in rails, we analyzed evidence for the accumulation of deleterious variants using SnpEff (Cingolani et al. 2012). We annotated mutations in protein-coding regions of the genome as synonymous (as a proxy of neutrality), missense (nonsynonymous amino acid substitution), or loss of function (LOF; e.g. premature stop codon). The latter category is expected to have a higher impact on the phenotype since the function of the protein could be affected, so we regarded LOF mutations as probable deleterious mutations. We found that each population had similar proportions of homozygous-derived genotypes for each mutation type (synonymous, missense mutation, and LOF; Fig. 4). For instance, the proportion of homozygous-derived genotypes under the LOF category was ~0.55 for Isabela, Pinta, Santiago, and Santa Cruz populations. These results suggest that inbreeding in Pinta, as suggested by their relatively long stretches of ROH that total ~140 Mb, have not increased their proportion of homozygous-derived genotypes with respect to the proportion of such genotypes from other populations (Fig. 4).

Discussion

Pattern of island diversification

Our biogeographic reconstruction supports the recurrent pattern of diversification shaped by several interisland dispersion events and the effect of geographic vicariance resulting from the flooding of the low island passes by seawater of a larger central landmass (Caccone et al. 2002; Arbogast et al. 2006; Parent and Crespi 2006; Steinfartz et al. 2009; Poulakakis et al. 2012, 2020). Notably, our genomic analyses suggest that the ancestral geographical distribution of rails comprised a single area that included the islands of Santiago, Santa Cruz, and Isabela (Fig. 1b; Black 1974). Importantly, the Isabela population was the first to diverge genetically from the rest (Fig. 1a and b), with populations from Santa Cruz and Santiago islands showing higher levels of connectivity. This result is supported by the absence of population structure as well as the lack of reciprocal monophyly between genomes from these two islands (Fig. 2, Supplementary Figs. S1 and S3). This biogeographic model coincides with the historical geographical distribution of these islands and their geologic history at the time of rail’s diversification (Black 1974; Geist et al. 2014; Karnauskas et al. 2017; Schwartz et al. 2018). A single land mass became divided into two or more individual islands resulting from the rising of sea levels around <1 mya (Fig. 1a and b; Poulakakis et al. 2012; Schwartz et al. 2018), where a widely distributed and connected ancestral lineage became the present-day populations on Isabela relative to Santiago, and Santa Cruz. Thus, the phylogenetic history of Galapagos rails indicates the role of vicariant events that resulted in population fragmentation and disruption of gene flow (Fig. 1a).

This model of lineage diversification is not uncommon as it has been recovered for Galapagos tortoises (Chelonoidis niger complex) (Poulakakis et al. 2012, 2020) and for vermilion flycatchers (Pyrocephalus nanus) with genetic clustering patterns remarkably similar to the grouping reported here for rails (Fig. 1a). In both vermilion flycatchers and rails, there is evidence that populations from Isabela diverged first, followed by Santa Cruz and then Santiago populations in the east of the archipelago (Carmi et al. 2016). This recurrent branching pattern across endemic Galapagos species suggests that Isabela may have been the first island to complete its isolation from the central region (Black 1974; Poulakakis et al. 2012, 2020; Carmi et al. 2016; Schwartz et al. 2018).

We found a lack of genetic differentiation among populations from the Isabela, Santiago, and Santa Cruz islands. Interestingly, these islands are currently isolated by water barriers, so population structure was expected. One possibility, as reported in other island birds (Reeve et al. 2023), is that past climatic oscillations promoted the formation of land bridges between nearby islands. Indeed, underwater mountains between Santa Cruz and Santiago arose from the current sea level during the cooling periods of the Pleistocene (Schwartz et al. 2018). The connection between these islands would have contributed to population gene flow consequently reducing genetic differentiation.

An isolated population on Pinta island

Our biogeographical reconstruction model suggests a dispersal event from the central landmass (e.g. I-C-S) to Pinta island (Fig. 1a and b), in contrast with the suggested origin and colonization of other Pinta endemic species. Pinta giant tortoises (Chelonoidis abingdoni) are reported to have originated from a dispersal event from Española Island located at the extreme south of the archipelago (Poulakakis et al. 2012). Pinta lava lizards (Microlophus pacificus) were likely the result of a dispersal event from Floreana via Isabela Island (Benavides et al. 2009) as Pinta Island never shared a terrestrial surface with any of these islands (Geist et al. 2014). The incorporation of historical rail samples from the islands of San Cristobal and Floreana, where rails have been extirpated, (see Carmi et al. 2016) could provide vital information to fully investigate island colonization hypotheses with a higher level of precision.

Rails on Pinta island showed unique genetic signatures supporting an independent evolutionary trajectory compared to the other islands. If our temporal estimation is correct, we hypothesize that the origin of this population is the result of a dispersal event from the central landmass after the original colonization of rail ancestors at 1 mya. At this time, Pinta Island was not formed and was only available for colonization around 700 kya (Geist et al. 2014). To date, it is uncertain how rails with reduced flight capacity were able to colonize the island of Pinta, which is ~80 km from the central archipelago (Black 1974). The loss of flight of volant colonizers after island colonization is a common pattern in insular rail species around the world (Slikas et al. 2002; Wright et al. 2016; Garcia-R and Matzke 2021). However, the process of flightlessness is not instantaneous and includes both the progressive reduction of flight muscles followed by an enlargement of hind limbs (Olson 1973; Wright et al. 2016) as well as changes in a variety of genes involved in developmental pathways (Burga et al. 2017). The Galapagos rail has been qualified as an “intermediate flyer” with morphological traits somewhere between the long-distance migrant relative (Black rail: L. jamaicensis) and its flightless relative (Inaccessible Island rail: Laterallus rogersi; Chaves et al. 2020; Durham 2021). It is possible that Galapagos rail individuals were fully volant for most of their early evolutionary history, potentially facilitating the movement across neighboring islands as well as the long-distance colonization of Pinta. This could explain the absence of population structure among the Santiago and Santa Cruz individuals (Fig. 2a) as well as the observed admixture between the Isabela and Santiago genomes (Fig. 2b and d).

On a more recent time scale, the ground vegetation of the highland region of Pinta, which rails previously inhabited, was cleared due to overgrazing by goats (Weber 1971), and no rails were detected on this island but rails reappeared shortly after a goat eradication program in 1971 (Kramer and Black 1970; Franklin et al. 1979). It was hypothesized that rail individuals on Pinta either recolonized this island from nearby islands or alternatively, it corresponds to a relic population that managed to survive this intense habitat modification and went undetected for years, most likely in small numbers (Kramer and Black 1970; Chaves et al. 2020). Rails are elusive small birds and even big-sized species such as Galapagos tortoises had remained undetected on Pinta (Snow 1964; Castro 1969; Franklin et al. 1979). Tortoises were thought to be extinct in the 1960s after multiple attempts to find them (Snow 1964; Castro 1969). Ten years later, however, one specimen named “Lonesome George” was found by accident (Vagvolgyi 1974). Our findings suggest that the current population of rails on Pinta is a relic population that went undetected for years and most likely in small numbers. Specifically, our analyses (PCA, Structure, and TreeMix) suggest that rails from Pinta are genetically divergent and have an independent evolutionary history compared to those on the other islands (Figs. 1 and 2). This finding is supported by the mitochondrial haplotype network as well as the phylogenetic tree (Fig. 2d and Supplementary Fig. S1) showing an independent origin of the Pinta. Our findings also suggest that rails remained unnoticed during previous surveys conducted by Kramer and Black (1970), suggesting that this population survived in small numbers with the current population rebound in large numbers that benefited from the intensive goat eradication program (Campbell et al. 2004).

Finally, our demographic model depicted by GONE (Fig. 1d) showed a population size increase in Galapagos rail populations at different time scales. For instance, the population size of Pinta had increased around 600 yrs ago (Fig. 1d). At this time, humid areas may have shifted to lower elevations (Restrepo et al. 2012), allowing rails to reach vegetation suitable for breeding and nesting (Franklin et al. 1979). During the last 30 yrs, all four rail populations showed a signal of drastic decline, presumably related to multiple events of goat introduction into the islands (Hoeck 1984).

Eradication of goats from Galapagos Islands: Did it work?

We found evidence of long ROH within populations which constituted 10% of the entire genome, indicative of recent population decline and inbreeding probably associated with the invasion of goats into the islands (Fig. 3b). Similar to other species with inbreeding, the Galapagos rail population, especially the Pinta population, showed a dominance of long blocks of ROH across their genomes (Kardos et al. 2018; Robinson et al. 2019).

Around 1970, the ground vegetation of the highland region of Pinta, which rails previously inhabited, was cleared due to overgrazing by goats (Weber 1971), and no rails were detected on this island but rails reappeared shortly after a goat eradication program in 1971. Rail populations that experienced sharp population declines due to habitat lost (e.g. Pinta) might be expected to have an overrepresentation of deleterious variants as observed in some populations of sea otters and wolves (Beichman et al. 2019; Robinson et al. 2019). However, the proportion of homozygous-derived mutations, including LOF mutations, across the different island populations was similar (Fig. 4). One possibility, like other island birds (Dussex et al. 2021), could be that individuals have a reduced mutational load that results from a combination of genetic drift and purging of deleterious mutations. However, to test the possibility of purging selection, the proportion of harmful variation of Galapagos rails should be compared to mainland rail closet species to demonstrate purging selection (Dussex et al. 2021). Alternately, our current approach to identifying deleterious mutations may lack the power to detect harmful mutations. Given that the distribution of fitness effects (DFE) of the Galapagos population is unknown our current method may not be able to discriminate strongly deleterious mutations from weakly deleterious mutations or detect the dominance coefficient of the different mutations (Robinson et al. 2023). This could mislead our interpretation of LOF mutations as being damaging mutations. For instance, some LOF mutations in some island populations could be subject to positive selection or missense mutations could be highly deleterious. Despite these caveats, growing evidence of empirical data has proven that our approach is useful in detecting genetic load and purging (Robinson et al. 2016, 2019). Like previous studies, we recommend incorporating historic samples (Robinson et al. 2018) before goat introduction and using computer simulations (Kyriazis et al. 2023) to better estimate the effect of demography on the Galapagos rail genetic diversity.

Another hypothesize is that the lack of genetic load in rails on Pinta relative to other islands was due to goats remaining in low-elevation regions for approximately 8 yrs before their incursion to higher grounds (Hamann 1979; Campbell et al. 2004), the time by which the eradication program had already been established on Pinta (Franklin et al. 1979). Another factor that could have mitigated severe inbreeding in rails could be the rapid recovery of their habitat (Hamann 1979). Although goats completely cleared the ground vegetation in the highlands, the vegetation was reestablished in less than 1 yr due to vegetative reproduction, as opposed to plant species found in the lowlands which recovered relatively slowly due to phases of growing, flowering, and fruiting (Hamann 1979). Thus, both the rapid eradication of goats and the re-establishment of highland vegetation on Pinta most likely allowed rails to recover without experiencing the genomic consequences of inbreeding (Hamann 1979; Campbell et al. 2004).

The goat eradication programs on Santiago and Isabela were established 26 yrs after the program in Pinta. However, the relatively large habitat of these islands coupled with multiple rail subpopulations may have mitigated the accumulation of deleterious mutations (Rosenberg 1990; Gibbs et al. 2003; Donlan et al. 2007). Particularly, Isabela is 15 times larger in area than Pinta (Black 1974). There are six different volcanoes across Isabela with enough elevation to support suitable habitats for rails. Similar to Galapagos tortoises, rails living on different volcanoes could be genetically differentiated (Beheregaray et al. 2004). Therefore, gene flow from other populations might have restored the genetic variation to small populations most affected by goat overgrazing (Rosenberg 1990; Donlan et al. 2007; Carrion et al. 2011). The elimination of approximately 140,000 goats from Isabela may have also mitigated population decline and inbreeding in rails on this island (Carrion et al. 2011). This outcome is contrary to the marked population structure of Galapagos weevils (Galapaganus conwayensis) resulting from host plant habitat fragmentation by grazing goats on Isabela (Sequeira et al. 2016). In the case of Santiago Island, the elimination of 79,000 goats led to a 16-fold increase in rail population size (Donlan et al. 2007; Carrion et al. 2011). The extirpation of goats from this island could explain the lack of a considerable proportion of harmful mutations among its rails. In contrast to rails from Santiago Island, the eradication of goats in Santa Cruz Island was minimal, with just 1,700 goats eliminated. Among the 11 individuals captured and sequenced from Santa Cruz in the presence study, 40% had a kinship coefficient of 0.25. This score is usually expected to be found between siblings, which could be an indication of a current and ongoing episode of inbreeding. However, high levels of kinship coefficient may also reflect the nature of the geographic sampling which could disproportionately sample members of the same pedigree (Chaves et al. 2020). Most samples were restricted to one valley near Cerro Crocket and corresponded to related individuals holding neighboring territories.

Conclusion

Our results indicate that the separation of the Galapagos central landmass <1 mya played a critical role in the diversification of endemic species in the Galapagos, including the Galapagos rail (Poulakakis et al. 2012, 2020; Carmi et al. 2016; Schwartz et al. 2018). We have shown that rails from Pinta Island have persisted as an isolated population despite the almost complete loss of habitat due to overgrazing by goats. Conservation efforts should focus on this population possessing a unique genomic signature and independent evolutionary history (Smith et al. 2014). The extinction of this small population would represent a significant loss of potential locally adapted alleles. Additionally, we found that the timely eradication of goats mitigated episodes of inbreeding among rail populations, with the one exception of the rails from Santa Cruz Island. Goats have not been completely eradicated from this island and the invasion of Cinchona trees may further affect rail populations on this island (Shriver et al. 2011). Different islands have a similar proportion of harmful variants which suggests that mixing distinct genetic lineages from different islands could result in genetic rescue. However, careful consideration must be taken into the estimation of DFE and simulations should be conducted to predict the outcome of moving individuals from one island to another (Dussex et al. 2021).

Supplementary Material

Supplementary material is available at Journal of Heredity Journal online.

Acknowledgments

We thank Annabel Beichman, Pooneh Kalhori, Klaus-peter Koepfli, Tom Smith, Kirk Lohmueller, and Chris Kyriazis for guidance in data analysis and improvement of the manuscript. We thank Audra Huffmeyer, Gaurav Kandlika, Laila Hualpa John Bates, and Elana Bachrach for proofreading the manuscript. We acknowledge the Vincent J. Coates Genomics Sequencing Laboratory at the University of California, Berkeley, and Fulgent Genetics in El Monte, California for DNA sequencing of samples. We thank SENESCYT for the financial support received during D.E.C.’s doctoral studies at UCLA.

Funding

This work was supported by the Lida Scott Brown Fellowship from the University of California, USA. Also, genome sequencing was possible thanks to the start-up funding from San Francisco State University.

Author contributions

Daniel Chavez (Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing), Taylor Taylor Hains (Conceptualization, Formal analysis, Methodology), Sebastian Espinoza-Ulloa (Formal analysis, Methodology), Robert Wayne (Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Visualization, Writing – original draft, Writing – review & editing), and Jaime Chaves (Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing)

Data availability

The raw sequencing reads of Galapagos rails generated for this study are available under BioProject PRJNA1057161. The Black rail (Laterallus jamaicensis) genome assembly (GCA_022605575.1; Hall et al. 2023) was obtained from https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_022605575.1/. The Inaccessible Island rail reference genome (Atlantisia rogersi) was downloaded from https://www.ncbi.nlm.nih.gov/assembly/GCA_013401215.1.The scripts used in this work are available on Github (https://github.com/dechavezv/2nd.paper.v2).

References