Understanding fatality patterns and sex ratios of Brazilian free-tailed bats (Tadarida brasiliensis) at wind energy facilities in western California and Texas

Carcass collection

For all tissue samples, we requested information on the collection date, location, estimated time since death, species, and sex (male, female, unknown). Species and sex were determined in the field using morphological characteristics. While sex was provided for some bat samples based on morphological characteristics, all bats were sexed using molecular methods for our analysis. Molecular methods allowed us to obtain more data by sexing samples that could not be identified in the field and provided more accurate understanding of the sex ratios of fatalities (Chipps et al., 2020).

DNA extraction

Portions of this text were previously published as part of a thesis (LiCari, 2023). We extracted DNA from the wing tissue samples using 96-well PCR plates, following the methods detailed in Ivanova, Dewaard & Herbert (2006), with an extra wash step. We placed a small piece of each sample (approx. 2–3 mm²) in an individual well in a 96-well PCR plate. Each well contained 50 μl of Vertebrate Lysis Buffer (VLB) (100 mM NaCl, 50 mM Tris-HCL pH 8.0, 10 mM EDTA pH 8.0, and 0.5% SDS) and 5 μl of Proteinase K (20 mg/ml). We then covered the 96-well PCR plate with a PCR mat and incubated it at 56 °C overnight to allow digestion. After the wing tissue was digested, we centrifuged the 96-well PCR plate at 1,500 g for 10 s to remove condensation from the PCR mat. Next, we added 100 μl of Binding Mix (BM), made from equal volumes of EtOH 96% and Binding Buffer (BB) (6 M GuSCN, 20 mM EDTA pH 8.0, 10 mM Tris-HCL pH 6.4, 4% Triton X-100) to each sample. We then covered the 96-well PCR plate with a new PCR mat, vortexed the plate vigorously for 15 s, and centrifuged the 96-well PCR plate at 1,000 g for 20 s. We transferred the lysate (about 150 μl) from the wells of the 96-well PCR plate into the wells of a Glass Fiber (GF) Plate (PALL1) placed on top of a square-well block. We sealed the GF Plate with a self-adhering cover. The covered GF Plate and square-well block were centrifuged at 5,000 g for 5 min to bind the DNA to the GF membrane. Three wash steps followed. For the first wash step, we added 180 μl of Protein Wash Buffer (PWB) (46.8 μl BB, 126 μl EtOH 96%, and 7.2 μl autoclaved dH₂O) to each well of the GF Plate, sealed the GF Plate with a new self-adhering cover, and centrifuged the GF Plate and square-well block at 5,000 g for 2 min. For the second wash step, we added 750 μl of Wash Buffer (WB) (60% EtOH 96%, 50 mM NaCl, 10 mM Tris-HCL pH 7.4, and 0.5 mM EDTA pH 8.0) to each well of the GF Plate, sealed the plate with a new self-adhering cover, and centrifuged the GF Plate and square-well block at 5,000 g for 5 min. For the third wash step, we added 500 μl of 70% ethanol to each well of the GF Plate, sealed the plate with a new self-adhering cover, and centrifuged the GF Plate and square-well block at 5,000 g for 3 min. After the wash steps, we centrifuged the GF Plate and square-well block at 5,000 g for 5 min to remove excess moisture from the samples. We then placed the GF Plate on top of a microplate equipped with a PALL collar and dispensed 70 μl of 10 mM Tris-HCL, pH 8.5 (pre-warmed to 56 °C) directly onto the membrane in each well of the GF Plate. We then incubated the GF Plate on top of a microplate equipped with a PALL collar at room temperature for 1 min. After incubating at room temperature, we covered the GF Plate with a new self-adhering cover. We then centrifuged the assembled plate at 5,000 g for 5 min to collect the DNA eluate into the bottom microplate. We then sealed the DNA plate with a new self-adhering cover and stored the extracted samples at −20 °C.

Molecular sex determination

We determined sex by amplifying the zfx and zfy introns found on the sex chromosomes using the Brazilian free-tailed primers in Korstian et al. (2013). All PCRs (10 μl) contained between 2 to 50 ng/μl of template DNA, 0.5 μM of each X-primer, 0.35 μM of each Y-primer, and 1X AccuStart II™ PCR SuperMix. The cycling parameters for the PCRs were one cycle at 95 °C for 15 min, followed by 35 cycles of 30 s at 94 °C, 1 min at 57 °C, and 30 s at 72 °C. We stained PCR products using Gel Green (Biotium), electrophoresed at 200 V in 1% agarose gel for 35 min and visualized the bands using blue-green LED light.

We used the number of bands present to determine sex. Females produce a single band corresponding to the X-chromosome intron (245 bp). Males produce two bands, one corresponding to the X-chromosome intron (245 bp) and one corresponding to the Y-chromosome intron (80 bp). We disregarded bands larger than ~300 bp when determining sex.

Data analysis

The groupings we analyzed were location (GH; GHN; LM; LV; WFA), year by location, and month by location. We were interested in determining whether the proportion of female fatalities differed from 0.5 for these groupings and whether the proportion of female fatalities changed over time (year and month) by location. Due to variation in search efforts across locations and time, aspects of data analysis differed between California and Texas, and different constraints were applied to each dataset. We will discuss these differences and provide more detailed information on the statistical tests used in subsequent paragraphs.

To compare sex ratios, we used carcass count data rather than fatality estimates, such as those obtained from GenEst, that take into account searcher efficiency, carcass persistence times, and density-weighted proportion (DWP) of the area searched (Dalthorp et al., 2018a, 2018b). Estimators like GenEst adjust the observed carcass counts to generate more accurate estimates of fatality that account for these biases. Nonetheless, we did not expect these biases to impact estimated sex ratios based on carcass counts. Previous studies suggested that carcass removal by scavengers at wind energy facilities was not sex-specific, so we expected each sex to remain available to be detected for equal amounts of time (Villegas-Patraca et al., 2012; DeVault et al., 2017). Similarly, we also did not expect searcher efficiency to influence the sex ratios. The Brazilian free-tailed bat is not sexually dimorphic (Twente, 1956; Davis, Herreid & Short, 1962; Kunz & Robson, 1995), so we expected each sex to be equally detectable to both detection dog teams and technicians. Searcher efficiency could change between seasons due to a change in environmental conditions that influenced detectability. This could result in sex ratios that are biased towards seasons when searches were less difficult, and more carcasses were collected. For example, the extreme heat of Texas summers could have decreased searcher efficiency during this season, which in turn could have led to an underrepresentation of data from the summer months when calculating annual sex ratios. This bias could be an issue for the Texas population, where we expected a difference in the sex ratio of fatalities across seasons. However, in Texas, we only analyzed data from August to October, so we did not analyze all seasons. Lastly, we did not expect the DWP to bias the sex ratios. Since the sexes are similar in mass (Twente, 1956; Davis, Herreid & Short, 1962; Kunz & Robson, 1995), we would not expect the carcass fall distributions to differ between the sexes.

For California, the groups used for analysis were year by location (GH, GHN) and month by location (GH, GHN). For the by-year analysis, we analyzed data only for years, with searches conducted from January–December. With this restriction, the two locations had no data overlap for the by-year analysis, excluding GH 2017, GH 2019, and GHN 2018 due to partial search years. We separated the month analysis by location and only analyzed months with combined totals of ≥30 bats pooled across years at the location (GH: September, October; GHN: May, June, July, August, September, October). Thirty was chosen as the minimum sample size to align with the central limit theorem. We used a one-proportion z-test for all groupings to determine if the proportion of female fatalities differed from 0.5 for location, years by location, and months by location. We conducted pairwise comparisons using multiple two-proportion z-tests, with a Bonferroni correction when applicable, to determine if the proportion of female fatalities differed between locations (GH vs. GHN), years by location (GH year vs. year; GHN year vs. year), and months by location (GH month vs. month; GHN month vs. month). In the California data analysis, we applied the Bonferroni correction to the GHN (month vs. month) analysis because the analysis required multiple pairwise comparisons within the group.

In Texas, search effort varied across locations and years. Since most searches were conducted at facilities from August to October, our analysis focused exclusively on bat fatalities during that period. This time frame also coincided with peak bat fatalities. We did not investigate the overall proportion of female fatalities by month due to differences in survey efforts across sites and seasons. The groups analyzed were location (LM, LV, and WFA) and year by location (LM 2017; LV 2017–2018; WFA 2021–2022). We used a one-proportion z-test for all groupings to determine if the proportion of female fatalities differed from 0.5. We used pairwise comparisons using multiple two-proportion z-tests, with a Bonferroni correction when applicable, to determine if the proportion of female fatalities differed between locations (LM vs. LV vs. WFA). For locations with two years of data, we used a two-proportion z-test to determine if the sex ratio changed between years at the location. In the Texas data analysis, we applied the Bonferroni correction to the comparisons between locations because the analysis required multiple pairwise comparisons within the group.

California

Although month-specific sample sizes were fairly small for the presumed non-migratory California population, we did not detect a difference in the proportion of female fatalities from 0.5 for year by location (1.14 > |z| > 0.44, 0.25 < p < 0.66, α = 0.05) or month by location (1.18 > |z| > 0.00, 0.24 < p < 1.00, α = 0.05), except for June at GHN (p_female = 0.33, n = 55, z = −2.56, p = 0.01, α = 0.05). We did not detect differences among groupings: locations (GH vs. GHN, z = −0.65, p = 0.52, α = 0.05), years by location (GH year vs. year, z = 0.19, p = 0.85, α = 0.05; GHN year vs. year, z = 0.76, p = 0.45, α = 0.05), and months by location (GH month vs. month, z = −0.31, p = 0.76, α = 0.05; GHN month vs. month, 1.66 > |z| > 0.11, 0.01 < p < 0.92, α = 0.008).

One unexpected finding from the analysis was that the California population experienced a difference in the timing of peak fatalities across the facilities surveyed. The shift in the timing of peak fatalities we observed in the raw carcass counts was also evident in the GenEst corrected seasonal fatality estimates in the PCM reports for GH (H. T. Harvey & Associates (HTH), 2021) and GHN (Great Basin Bird Observatory (GBBO) & H. T. Harvey & Associates (HTH), 2022). For the first two complete years of surveys at GH (2017 and 2018), fatalities peaked between September and October, as expected (American Wind Wildlife Institute (AWWI), 2020; Figs. 1B, 1C). In 2019 and 2020, however, fatalities at GHN peaked between May and June (Figs. 1F, 1G).

Texas

For the migratory Texas population, the proportion of female fatalities varied among the three wind energy facilities. The proportion of female fatalities at Los Mirasoles did not differ from 0.5 (p_female = 0.56, n = 124, z = 1.44, p = 0.150, α = 0.05; Fig. 2A). In contrast, fatalities at Los Vientos were female-skewed (p_female = 0.55, n = 516, z = 2.25 p = 0.024; α = 0.05; Fig. 2A); and fatalities at WFA were male-skewed (p_female = 0.41, n = 341, z = −3.20, p = 0.001, α = 0.05; Fig. 2A).

When comparing the proportion of female fatalities between the three locations, the proportion of female fatalities at WFA was lower than the proportion of female fatalities at Los Mirasoles (z = −2.88, p = 0.004, α = 0.017) and Los Vientos (z = −3.90, p = 0.0001, α = 0.017; Fig. 2B).

The proportion of female fatalities varied between years at wind energy facilities in Texas (Fig. 3). At Los Vientos, fatalities were female-skewed in 2017 (p_female = 0.60, n = 285, z = 3.38, p = 0.0007, α = 0.05) but not in 2018 (p_female = 0.49, n = 230, z = −0.40, p = 0.69, α = 0.05). The proportion of female fatalities at Los Vientos decreased by 0.11 between 2017 and 2018 (z = −2.55, p = 0.01, α = 0.05). Lastly, the peak in Brazilian free-tailed bat fatalities consistently occurred in late summer and early autumn at the Texas wind energy facilities included in this study (Fig. 4).

California

The overall lack of bias in the sex ratio of Brazilian free-tailed bats at two wind energy facilities in California is consistent with what we know about this population’s presumed non-migratory behavior. The high proportion of male fatalities at one site (GHN) in June could be due to the local movement of male Brazilian free-tailed bats leaving roosting caves due to overcrowding as females move into caves and prepare for pupping in June (Davis, Herreid & Short, 1962). However, this pattern was not observed in June for GH, which is only 4 km from GHN, illustrated in Fig. 5. Under the assumption that the sex ratio of the non-migratory Brazilian free-tailed bat population in California is 50:50 (F:M), as seen in other non-migratory bat species (Davis, 1969; Cheng et al., 2019), it appeared that the proportion of female fatalities reflects the proportion of females in the population. Although our analysis does not specifically address this question, the absence of sex-specific bat-wind turbine behaviors in this non-migratory population could also explain these findings.

We considered two possible explanations for the shift in the timing of peak fatalities, which occurred in September and October at GH and May and June at GHN, as seen in both raw carcass counts and in the GenEst-corrected seasonal fatality monitoring estimates in the PCM reports for GH (H. T. Harvey & Associates (HTH), 2021) and GHN (Great Basin Bird Observatory (GBBO) & H. T. Harvey & Associates (HTH), 2022). One possibility proposed in the PCM report (Great Basin Bird Observatory (GBBO) & H. T. Harvey & Associates (HTH), 2022) is that GHN is near a reservoir used for foraging (Fig. 5), as Brazilian free-tailed bats are known to forage over large bodies of fresh lacustrine water during summer and fall (Johnston, 2013). Nevertheless, it was interesting that GH did not exhibit the same shift in 2017 and 2018, despite its close proximity to GHN (4 km apart), as depicted in Fig. 5, given that Brazilian free-tailed bats are known to fly over 56 km from their roosting sites when foraging at night (Best, 2009). Another possible factor that could have influenced this shift in timing is a change in moth abundance. Moths are a significant component of the Brazilian free-tailed bat’s diet (Kunz & Robson, 1995; Whitaker, Neefus & Kunz, 1996; Lee & McCracken, 2002), and bats have even been observed changing their foraging behavior in response to moth migrations (Krauel et al., 2018). Outbreaks of the California oak moth (Phryganidia californica) in May and June of 2019 and 2020 occurred in the same county in which the wind energy facilities in this study are located (Oboyski, 2020). The California oak moth experiences population outbreaks at irregular intervals that last one to two years (Burke & Herbert, 1920; Swiecki & Bernhardt, 2006). It is possible that GH also experienced a change in the timing of peak fatalities in 2019, but sample sizes were too small to detect a shift (n = 32). DNA barcoding of the stomach contents of the bat carcasses collected at GH and GHN could be used to determine if there were differences in their diets between locations (waterbody influence) or time periods (changes in moth abundance across years). Such data would help us better understand which external factors influence variation in bat fatalities at California’s wind energy facilities.

Texas

Estimates of the proportions of female Brazilian free-tailed bat fatalities varied across wind energy facilities and across time in Texas. The proportion of female fatalities differed between WFA and the two more southern locations, LM and LV, which are separated by only 2 km at the closest edge. Wind Farm A is approximately 200 km from LM and LV; it is possible that the sex ratio in the migratory Brazilian free-tailed bat population is not a uniform 9:1 (F:M; Davis, Herreid & Short, 1962) across the summer range. Differences in local sex ratios could have caused a higher proportion of male fatalities at WFA compared to the other wind energy facilities.

Previous studies regarding the Brazilian free-tailed bat sex ratio by roost type (bridge and cave) reported a difference in the proportion of females present in bridge roosts and cave roosts (Table S1). The average proportion of females at five surveyed bridge roosts was 0.44 (Sgro & Wilkins, 2003; Turmelle et al., 2010; Martinez, 2015). In contrast, the average proportion of females at 12 surveyed cave roosts was 0.73 (Cagle, 1950; Twente, 1956; Davis, Herreid & Short, 1962; Rogers, 1972; Thies, 1993; Turmelle et al., 2010; Danielson et al., 2022). Sex-specific roosting patterns may explain these differences (Hermanson & Wilkins, 1986; Sgro & Wilkins, 2003). For geographic reference, we have included known Brazilian free-tailed roost locations (bridges and caves), known bat roosts with assumed Brazilian free-tailed bats within the Texas range of the migratory population (S Weaver, 2023, personal communication), and roosts along the Mexico-United States border in the map provided in Fig. 6.

Because of the sex bias in roost type utilization in this population, roost proximity could influence the proportion of female fatalities at WFA. Wind Farm A was adjacent to four bridge roosts, all within 50 km (Fig. 6), which could result in an estimated proportion of male fatalities that is higher than expected based on the overall population sex ratio. Additionally, WFA is 112.7 km north of a bridge roost location that males are known to use during the winter (Mink, 2012). Proximity to this roost could be causing the male skew we observed in September and October, as males travel from their bridge roosts north of WFA to this winter roost in southern Texas (Fig. 6). However, because these sites were surveyed in different years, we cannot exclude the potential impact of temporal variation in the population sex ratio on the observed differences in the proportion of female fatalities at the southern wind energy facilities (LM and LV) and the northern facility (WFA).

Temporal effects could also influence the differences in the sex ratio of fatalities across locations, as LM and LV were surveyed from 2017 to 2018, and WFA was surveyed from 2021 to 2022. We observed a change in the sex ratio over time at locations with multiple years of surveys, with LV showing a decrease in the proportion of female fatalities between the two years. Further research is required to ascertain if the decrease in the proportion of female fatalities is consistent across facilities or the result of inherent between-year variation and sampling limitations.