Purpose and Scope

Although counts have been completed annually at most occupied wetlands since the 1980s, monitoring of movement and survivorship of released juvenile rails had not occurred until very recently (Zembal and others, 2017; Sawyer, 2024; Sawyer and Conway, in press). In addition, genetic monitoring of wild populations and genetic assessment of captive birds have been lacking until this study. Therefore, little is known about the cumulative effects of releases on population genetic structure and diversity of recipient populations. To address these uncertainties, we developed a set of microsatellite markers to allow for genetic monitoring of wild and captive rails. We evaluated the recent (2020–22) genetic population structure and diversity of rail populations throughout their U.S. range. We also compared the recent genetic structure to the pre-augmentation structure by comparing recent blood samples with blood samples available from the initial 1989 genetic surveys. Moreover, we examined genetic connectivity and diversity across the subspecies’ U.S. range to identify extant populations with high genetic diversity that could be considered for future captive-rearing sources.

Field Sampling

We visited and captured wild individuals at 17 wetlands throughout the U.S. range for genetic sampling and banding between 2020 and 2022. Sites were visited during the breeding season, roughly between April and September of each year. We used carpet traps (Harrity and Conway, 2020) with a broadcast of Ridgway’s Rail vocalizations to lure rails to the carpet traps (Pickens and King, 2013; Harrity and Conway, 2020). We removed rails from carpet traps immediately after capture, and we measured, weighed, photographed, and attached a federal leg band (Smith, 2013) to each rail. We collected blood samples from each captured rail via metatarsal venipuncture using a sterile, 26-gauge needle and transferred to a GenSaver 2.0 (AHLSTROM, Escondido, California, cat no. 8.566.0002.B-N) blood card with a non-heparinized capillary tube (Thermo Fisher Scientific, Waltham, Massachusetts, cat no. 22-260943). We then released rails at the capture location. All fieldwork was authorized following guidelines specified in Federal and State permits held by C. Conway (Federal Endangered Species Permit TE039466; Bird Banding Permit #22524; California Memorandum of Understanding (SCP-S-193610002-20008-001), and as approved by the University of Idaho Institutional Animal Care and Use Committee (2015-51).

Marker Development

Microsatellite libraries were developed for R. obsoletus at Cornell University’s Evolutionary Genetic Core Facility using genomic DNA extracted from four individuals. The Evolutionary Genetic Core Facility sequenced a tetrameric, enriched genomic library on an Illumina MiSeq with paired 250 base-pair reads (Nali and others, 2014), used SeqMan NGen (version 11, DNAStar, Madison, Wisconsin) to generate a de novo assembly from the paired fastq files (raw data), and used the program msatcommander 1.0.8_beta (Faircloth, 2008) to scan for candidate microsatellite loci and design primer pairs. To design a panel of highly multiplexed microsatellite markers, we randomly selected and evaluated approximately 500 candidate microsatellite loci (500 forward primers tagged at the 5-prime end with the sequence TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG, and 500 reverse primers, tagged at the 5’ end with the sequence GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG) for multiplex polymerase chain reaction (PCR) suitability using Multiple Primer Analyzer (MPA; Thermo Fisher Scientific, Sunnyvale, California). We used the MPA output in the package igraph (Csárdi and Nepusz, 2006) to cluster the loci into an arrangement that would minimize primer-dimer formation. This process resulted in four multiplexes, composed of 30–40 loci each. These loci were individually amplified in two individuals, and only loci with successful amplification in both samples (as confirmed by gel electrophoresis) were retained in the final multiplexes. All samples were genotyped using the resulting panel of 108 loci (table 1.3) arranged into four multiplexes (Mpx1-4), as described in the following section.

DNA Extraction, Amplification and Sequencing

We extracted genomic DNA from blood cards or capillary tubes using the Puregene kit (QIAGEN, Germantown, Maryland) according to the manufacturer’s protocol, with minor modifications including the addition of Proteinase K to cell lysis with an overnight incubation at 58 degrees Celsius (°C), and final resuspension in 100 microliter (µL) Tris Low ethylenediaminetetraacetic acid (EDTA; TLE) buffer (10 millimolar [mM] Tris, 0.1 mM EDTA, pH 8.0). We quantified extractions using Qubit Broad Range (Thermo Fisher Scientific) and standardized to 10–40 nanograms per microliter (ng/µL) before amplification with the Type-it Microsatellite PCR Kit (QIAGEN). We amplified loci by using four primer cocktails (Mpx1-4; table 1.3), with each primer at a 1.6 micromolar (µM) concentration in the primer cocktail. Each of four 10 µL PCR reactions contained 5 µL 2X Type-it Master Mix, 1 µL of Mpx1, Mpx2, Mpx3 or Mpx4, and 15–60 ng/µL DNA. Amplifications included 30 cycles of 95 °C for 5 minutes, 94 °C for 30 seconds, 56 °C for 1.5 minutes, 72 °C for 1.5 minutes, followed by a 12 °C hold. Upon completion, the four multiplexed PCRs per sample were pooled together, and the pooled PCR product was barcoded using Nextera N5/600 and N7/800 indexes to produce individual dual-indexed amplicon libraries for each sample. Individual sample libraries were then pooled together into one tube per 96-well plate and bead-cleaned to remove primer dimers. Pooled and bead-cleaned libraries from each plate of sample libraries were combined in equimolar proportions and sent for sequencing at MedGenome, Inc. (Foster City, California) on the NovaSeq 6000 (Illumina, San Diego, California), using the Illumina SP 300 cycle reagent kit v1.5.

Bioinformatics

We used the Python script amplicon (https://bitbucket.org/cornell_bioinformatics/amplicon) to extract reads from the Illumina runs and assign them to the appropriate locus and individual. Specifically, the script (1) trims adapters and low-quality reads, (2) creates contigs from overlapping reads (for paired-end sequencing), (3) identifies reads corresponding to each locus, (4) collapses identical reads for each individual, and (5) identifies the top two haplotypes for individuals at all loci (in other words, their diploid genotypes). We used the default options except for the following parameters: -c 1 (minimum number of samples per haplotype), -a 0.001 (minimum minor allele frequency), -l 75 (minimum haplotype length), -r 5 (maximum read count ratio between the two alleles in each sample). We then calculated the total number of reads per locus per individual, whether a locus was heterozygous or homozygous (and if heterozygous, the minimum number of minor allele reads per locus per individual). We scored loci as missing data if the total number of reads was less than (<) 200. Heterozygous loci were recoded as homozygous if the minor allele read count was low (<300), or if the total number of reads was low (<500). Before population genetic analyses, we used the R packages adegenet v. 2.1.10 (Jombart, 2008) and poppr v. 2.9.3 (Kamvar and others, 2014) to assess the quality of loci and samples using several filters. First, we removed any locus with greater than 10 percent missing data, and then removed any individual samples with greater than 10 percent missing data. Next, we applied a minor allele frequency cutoff (MAF=0.01) to identify monomorphic or uninformative loci. Once these loci and samples were removed, we used the R package genepop v. 1.2.2 (Rousset, 2008) to evaluate the dataset for linkage disequilibrium using the exact test for genotypic linkage disequilibrium with 10,000 dememorizations and 5,000 iterations; the significance of linkage disequilibrium was confirmed for loci with p-values below 0.0001. We also used the method described by Brookfield (1996) to estimate the frequency of null alleles for each locus with the R package popgenreport (Adamack and Gruber, 2014). We retained loci with null allele frequencies less than 0.2, following the recommendations of Dakin and Avise (2004).

Population Genetic Dataset

During field sampling, we captured and sampled hatch-year and adult rails. However, we removed hatch-year birds from the population structure and diversity analyses to avoid biases resulting from unequal sampling of family groups and to focus on the adult breeding populations present at the time of sampling. We included captive adults used in the breeding program to represent the captive “population” (hereafter “captive breeders”). We separated the captive breeders into two temporal groups: (1) parents of the captive offspring released before and during the wild sampling period (2019–21; group I), and (2) captive breeders held in the breeding program at the time of this report (2023–24; group II). Group I birds were included along with wild birds in structure and gene flow analyses to help evaluate the influence of the breeding program on genetic structure and diversity. Diversity metrics were calculated for groups I and II to provide information relevant to the breeders in captive breeding facilities at the time of this report. We analyzed the 1989 baseline samples separately from recent samples to compare population structure and diversity pre- and post-augmentation.

Loci were screened for deviations from Hardy-Weinberg equilibrium (HWE) at four wetland sites with greater than 20 samples (Newport Bay, Batiquitos Lagoon, San Elijo Lagoon, Tijuana Slough NWR) using an exact test based on 1,000 Monte Carlo permutations of alleles (Guo and Thompson, 1992) and applying the Benjamini and Yekutieli (2001) correction for multiple tests. We removed loci if they deviated significantly (corrected p-value less than 0.05) from HWE at three or more sites.

Population Structure and Gene Flow

We used multiple methods to assess population structure. First, we used STRUCTURE (Pritchard and others, 2000) to determine the supported number of genetic clusters (K) that conform to populations in genetic equilibrium. We specified a range for the maximum number of clusters that individuals could be assigned (K=1–10) and completed 10 replicate runs per K using 500,000 iterations of the Markov chain Monte Carlo algorithm following a burn-in of 500,000 iterations to verify consistency across chains. The optimal K was inferred by comparing the results from the maximum mean log-posterior probability for K estimated by STRUCTURE and the change in K (∆ K) criterion (Evanno and others, 2005). Second, we used principal component analysis (PCA) to visualize genotypes in multidimensional space with adegenet v2.1.10 (Jombart, 2008), in R v4.1.2 (R Core Team, 2018). We used the program PopCluster (Wang, 2022a) to estimate gene flow among populations. PopCluster provides estimates of recent gene flow rates (last 3 generations) from an admixture model. We first evaluated up to 10 clusters (K) with 20 replicate runs. After selecting the optimal K, we ran the PopCluster model with migration for 20 replicate runs to estimate gene flow rates among clusters from the individual admixture estimates.

We calculated allelic richness (Ar), private allelic richness (PAr), observed heterozygosity (Ho) and unbiased expected heterozygosity (He), and inbreeding coefficients across marsh sites, and clusters and groups of captive breeders. There was some geographic overlap between cluster assignments in Mission Bay (Kendall-Frost Mission Bay Marsh Reserve and San Diego River). For the purpose of reporting genetic diversity indices by cluster, we grouped these two wetlands in the south San Diego County cluster. The effective population size (N_e) was estimated in NeEstimator v2 (Do and others, 2014) for each cluster and period. We used the linkage disequilibrium method with monogamy and a minimum allele frequency of 0.02, and calculated 95-percent confidence intervals (CI) of point estimates by jackknifing across samples.

We compared genetic differentiation (F_ST), relatedness (R), allelic richness (Ar), and unbiased expected heterozygosity (He) between the baseline and recent sample periods by using group comparisons in FSTAT v2.9.4 (Goudet, 2001), with p-values derived from 10,000 permutations. We restricted our analysis to the three wetlands that were sampled in both periods (Mugu Lagoon, Newport Bay, Tijuana Slough NWR). During the time of our field sampling, only a handful of birds were observed in Seal Beach; we did not pursue sampling there to avoid disturbing the remaining rails. We also ran a PCA across paired wetlands to visualize any changes in genetic clustering over time.

Decision Framework for Genetic Rescue

Following the decision framework presented in Frankham and others (2017), we assessed whether populations met certain criteria indicating genetic erosion and whether genetic rescue could improve genetic diversity in local populations. We calculated the mean inbreeding coefficient (F):

F values greater than 0.1 indicate that genetic diversity is sufficiently higher in donor population(s) to benefit the receiver population (Frankham and others, 2017). We used our estimates of expected heterozygosity to calculate F for each regional cluster in relation to the following donors: (1) captive breeders, (2) Orange County cluster, (3) north San Diego County cluster, and (4) south San Diego County cluster.

Comparing Coancestry and Inbreeding Coefficients from Studbook and Genetic Data

We calculated pedigree-based coancestry and inbreeding coefficients for captive breeders using the R package kinship2 (Sinnwell and others, 2014). We then estimated the genetic-based relatedness and inbreeding coefficients for captive breeders in the software EMIBD9 (Wang, 2022b). The EMIBD9 software implements a likelihood expectation maximization (EM) method, updating allele frequencies and identity-by-descent coefficients for each pair of sampled individuals until convergence. The EM method estimates relatedness and allele frequencies simultaneously from a small sample of genotypes, in contrast to traditional methods, which rely on unbiased allele frequencies obtained from a large sample of unrelated genotypes (for example, relatedness presented in table 1). We then compared the productivity, pedigree-based metrics and genetic-based metrics for recent breeding pairs.

Recent Population Structure

Structure analyses best supported three genetic clusters across the range (fig. 4) that roughly corresponded to sampled regions ([1] Orange County plus Mugu Lagoon and the captive breeders, [2] north San Diego County plus San Diego River, and [3] Kendall-Frost Mission Bay Marsh Reserve and south San Diego County; fig. 5). Individuals of mixed assignment were reported in all clusters, indicating recent or ongoing dispersal and gene flow occur directly among wetlands, or that gene flow is facilitated through the efforts of the captive breeding program. Principal component analysis also grouped individuals into three regional clusters (fig. 6A) along axes 1 (7 percent of the total genetic variation) and 2 (6 percent of the variation), while axis 3 (5 percent of the variation) separated individuals within marsh sites, particularly within the Orange County cluster (fig. 6B). PopCluster also supported three clusters and estimated recent gene flow rates among clusters ranging from 4.1 percent (from north San Diego County to Orange County) up to 23.5 percent (from Orange County to north San Diego County; table 2). Recent gene flow estimates among clusters were high, especially into the north San Diego County cluster. However, natural versus augmented levels of gene flow are difficult to separate in this system given the approximately 20 years (10–20 generations) of captive breeding and releases before genetic monitoring efforts. Higher rates of recent gene flow (last 3 generations) from the Orange County cluster (the source of the captive program) into the other two clusters could reflect these captive release efforts. Recent telemetry data indicate rail movement is usually localized. In a group of transmittered wild (N=42) and captive released (N=46) hatch year rails, only one captive rail moved between wetlands (from Tijuana Slough NWR to the San Diego Bay NWR South Bay Unit, about 4 kilometers [km]); all other rails stayed close to the initial capture locations (Sawyer, 2024). Similar average distances are reported from earlier studies, although occasional long-distance movements of up to 258 km have been recorded (U.S. Fish and Wildlife Service, 2020). The structuring of individual wetlands into three broader genetic populations that appear to be connected by moderate levels of gene flow provides important context for population management, supporting an inclusive regional approach consistent with genetic structure, rather than focused on individual wetlands as independent populations.

Comparisons with Historical Samples

Comparing our recent samples to the 1989 baseline samples across the same set of wetlands suggests that genetic differentiation has declined slightly, relatedness has increased, and allelic richness has decreased over time (table 3). Although none of these changes were statistically significant, the direction of these respective measures is consistent with a small decline in overall genetic diversity in rail populations, which is noteworthy given the increase in sample size in two of three wetlands. The slight decline in F_ST may reflect the effect of the captive-bred individuals being sourced from a single site (Newport Bay) and released throughout the range or could reflect an increase in naturally occurring dispersal and gene flow among regions, facilitated by the increased population sizes in the center of the range. The PCA of baseline and contemporary samples separated sites spatially along the primary axis (fig. 7A). The second axis separated the temporal sampling periods at Tijuana Slough NWR and the third separated the temporal sampling periods within Mugu Lagoon and Newport Bay (figs. 7A, B). Mugu Lagoon seems to be the most distinctive over time, with non-overlapping point clouds (fig. 7B). Genetic differentiation over time is consistent with genetic drift (loss of genetic diversity over time), which is more extreme in smaller and more isolated populations.

Genetic Diversity and Effective Population Size

By all measures, north San Diego County has the highest genetic diversity of all sampled regions, followed by south San Diego County, Orange County and, having the lowest genetic diversity, Mugu Lagoon (table 1). Low genetic diversity at Mugu Lagoon could reflect its position at the northern range edge and consistently low survey numbers. Despite augmentation attempts with more than 100 captive-reared birds between 2001 and 2009, the maximum number of pairs observed at Mugu Lagoon during the last two decades was in the low 20s, and only a handful of pairs were observed during the past few years (fig. 2; Zembal and others, 2024). Newport Bay also appears to have relatively low levels of genetic diversity compared with its baseline sample. Counts during annual surveys have rapidly declined since 2017 in Newport Bay and in surrounding wetlands in Orange County. Declining numbers here have been attributed to increasing tidal inundation (Zembal and others, 2024). We could not estimate effective population size (N_e) at Mugu Lagoon due to low sample size. Among the other three regions, the contemporary N_e point estimate was lowest in Orange County and highest in north San Diego County (table 4), which is consistent with all other diversity metrics. Contemporary N_e point estimates were lower than baseline samples, although CIs overlapped (table 4). General guidelines suggest N_e should be greater than 50–100 to avoid inbreeding, and greater than 500–1,000 to preserve allelic richness and long-term adaptive potential (Frankham and others, 2014). Orange County may be at or below these lower thresholds (upper 95-percent CI=113), whereas north San Diego (upper 95-percent CI=486) and south San Diego County (upper 95-percent CI=559) may be at or below the upper thresholds.

Genetic Rescue

Frankham and others (2017) provides decision tables for determining whether a population could benefit from genetic rescue and state that appropriate source populations should have higher heterozygosity than the receiver population (F>0.1). Because it was estimated to have greater heterozygosity than the other genetic clusters, the north San Diego County cluster could be a genetically beneficial source for all other regions examined, producing F>0.1 in the receiver populations (table 5). Mugu Lagoon, with the lowest heterozygosity of any site, could benefit from genetic rescue from any other source (table 5). Finally, because the captive breeders have low diversity when compared to wild populations, augmentation results in negative F for all sites except for Mugu Lagoon (table 5). New source populations could improve diversity in the captive breeding program (see the “Managing Genetic Diversity in the Captive Program” section).

Whether or not wild regional clusters could benefit from genetic rescue can also be assessed with information about population size and isolation. Although Mugu Lagoon and the Orange County cluster have low or declining survey numbers, low effective population sizes, and are geographically more isolated, the north and south San Diego County clusters have larger survey numbers based on recent call-broadcast surveys (fig. 3) and higher effective population sizes. Although gene flow estimates among clusters were high, augmentation through the captive release program could account for some of this, and rates were lowest into the Orange County cluster.

Managing Genetic Diversity in the Captive Program

Adjustments to captive rearing source populations and release strategies, informed by new empirical estimates of population genetic diversity and structure, could help preserve genetic diversity. The pool of captive breeders has lower genetic diversity than all wild genetic clusters except for Mugu Lagoon. Therefore, the recently released hatch year birds likely added little or no genetic diversity benefit to the receiver populations into which they have been released (table 5). Two factors may contribute to this. First, breeding birds for the captive program have been consistently sourced from one wetland across the range (Newport Bay). However, this wetland has recently declined in size and has low genetic diversity and low effective population size, suggesting it may benefit from genetic rescue itself (tables 4, 5). Second, the captive breeding program is small, composed of up to six pairs annually. Some of the recent breeding birds have high inbreeding coefficients and some pairs have elevated (non-zero) genetic relatedness, despite efforts to minimize pairings between known relatives based on the pedigree (table 6). These genetic estimates could indicate non-zero relatedness among the wild ancestors. Given the small number of breeding birds in the captive program at any one time, efforts to rotate in wild birds more frequently could help to incorporate new genetic diversity. Retaining later generations of offspring in the breeding program could increase relatedness, depending on pairings. Large differences in productivity among breeding pairs may also skew the genetic makeup of captive-released cohorts. This could be reduced by limiting the number of clutches produced by each captive pair each season. Limiting breeding windows, especially to the beginning of the season, may also help increase the probability of survival for captive-released juvenile rails. Analysis of telemetry data indicated that captive rails released early in the summer had higher survival rates than those released later (Sawyer, 2024; Sawyer and Conway, in press). Finally, ensuring receiver sites receive a mix of clutches produced by unrelated pairs could decrease the overall relatedness of birds released at a single site and season.

Wetlands in north San Diego County have the highest heterozygosity, allelic richness, and private allelic richness across all surveyed regions and the lowest relatedness. Sourcing birds or eggs from the larger wetlands within the north San Diego County cluster could provide the greatest increase to the genetic diversity and representation within the captive breeding population for future population augmentation (table 4). Because relatedness values were generally higher within than among wetlands even within the same regional clusters, pairing birds sourced from different wetlands instead of a single wetland could also help reduce the chances of including closely related birds in the captive program. Finally, genotyping all candidate parents could directly estimate genetic relatedness and suggest pairings to minimize inbreeding.

Wetland Restoration

Given that wetlands in north San Diego County appeared largely unoccupied before the mid-2000s, it could be possible that a combination of habitat restoration coupled with captive releases (fig. 3) are responsible for the increase in numbers of pairs and high genetic diversity in north San Diego County. In addition, given estimated gene flow rates of 8–11 percent between south and north San Diego County genetic clusters, it is possible that natural dispersal of wild birds may be sufficiently high to maintain genetic diversity and connectivity across this part of the subspecies’ range.

Although opportunities to restore wetlands may be rare throughout the northern part of the subspecies’ range, restored wetlands could provide more stepping stones for increased connectivity. In the more immediate time frame, our genetic analyses suggest that Orange County and Mugu Lagoon clusters could benefit from augmentation and genetic rescue from a higher-diversity source population.

Another important factor in maintaining high diversity is retaining large populations to minimize the erosion of local genetic diversity. Habitat management and restoration can assist in maintaining large populations and could become even more critical given predicted sea-level rise, which may threaten wetland habitat in areas without sufficient upland habitat for marsh retreat (Osland and others, 2022), and may already be affecting the population at Newport Bay (Zembal and others, 2024). Models of California wetland vulnerabilities to sea-level rise, including three marshes occupied by rails (Newport Bay, Sweetwater, and Tijuana Slough NWR) predicted significant loss of high and middle marsh habitat by 2050 and between 50- and 100-percent conversion to bare mudflats by 2100 under moderate to high sea-level rise scenarios (Thorne and others, 2018). Survival of juvenile rails is affected by elevation, and the timing and water level at high tide (Sawyer, 2024; Sawyer and Conway, in press). The abundance of raptors may also have a negative effect on survival, especially for captive-released rails (Sawyer, 2024; Sawyer and Conway, in press). A recent 5-year study of mortality in California Ridgway’s Rails in San Francisco Bay indicated that avian predators accounted for most of the observed mortalities (Casazza and others, 2016).