Renewable energies and biodiversity: impact of ground-mounted solar photovoltaic sites on bat activity

Sampling design We implemented a paired study design across 19 solar PV sites to assess whether bat species richness and activity were higher in fields and along boundary habitats that contained PV panels, compared with ‘empty’, matched control sites. This resulted in 19 sampling points for solar boundary habitat, 19 for solar open habitat, 19 for control boundary habitat, and 19 for control open habitat. All sites were located in south-west England, where the highest concentration of solar PV sites and greatest bat species richness in the UK coincide (Mathews et al. 2018; Department for Business Energy and Industrial Strategy 2021). Where private land was entered, permissions were granted by land owners and the relevant solar farm companies. No ethical approval was required for this study as we passively monitored bats through acoustic recordings. The control sites were within the same land management boundary as the solar PV site, and matched as closely as possible in plot size, habitat type, land use and boundary habitats. There was no difference in the average size of solar PV and control fields (solar PV mean = 59.6 ha, SD = 32.0; control mean = 53.2 ha, SD = 28.4; paired t-test: t(18) = 1.3, P = 0.203) (see Appendix S3 in Supporting Information). All solar PV sites were on grassland that was either grazed or managed through mowing or were on cut arable crops. Field boundaries corresponded to hedgerows, treelines, woodland, or vegetated ditches and were exactly matched. The paired fields were a minimum of 500 m apart and not adjacent to each other to maximise the chances of obtaining independent data within comparable landscapes (Froidevaux, Louboutin & Jones 2017). Bat echolocation call recording and species identification Fieldwork was completed between July and October 2019 and the same period in 2020. Bat activity was monitored for seven consecutive nights at each site (30 minutes before sunset to 30 minutes after sunrise), simultaneously across the four locations (open and boundary habitats within the field with solar PV panels and paired control). Recordings were made using SM3 bat detectors (Wildlife Acoustics, Inc., Maynard, MA, USA). All detector microphones (SMM-U2f (frequency response +/- 6dB 20-100 kHz see https://www.wildcare.co.uk/amfile/file/download/file/56/product/94208/), Wildlife Acoustics)) were elevated to 1.27 m using identical tripods. Detectors were set to auto trigger between 8–120 kHz, 1–88 dB and recorded for a maximum of 10 s (384 kHz, sampling rate). A detector was placed within the centre of the control and solar fields, and along the associated boundary habitats of the control and solar fields. Detectors recording the open and boundary habitat within the solar and control field were a minimum of 50 m apart. Sampling took place during optimal weather condition for bats to forage (i.e., no rain, low wind speed and temperature >10°C). The mean (+ SD) temperature at dusk over the recording period was 16.2 ± 3.1 °C (https://www.timeanddate.com/weather/). Sound files were analysed using zero crossing software Kaleidoscope Pro (v. 5.4.1, Wildlife Acoustics, Inc.) with Bats Of Europe Classifiers (United Kingdom) (v. 5.4.0) selected. All 10 s recordings were automatically scanned and the call sequences were identified and then manually checked to confirm the species (Barbastella barbastellus, Eptesicus serotinus, Pipistrellus nathusii, P. pipistrellus, P. pygmaeus, Rhinolophus ferrumequinum and R. hipposideros) or species group (Nyctalus spp., Myotis spp., Plecotus spp.). The grouping of Myotis spp. is widely used due to the difficulty of separating the echolocation calls of the different species (Russ 2012). Similarly Nyctalus noctula and N. leisleri, as well as Plecotus auritus and P. austriacus could not always be separated so these calls were grouped as Nyctalus spp. and Plecotus spp., respectively. All files which Kaleidoscope Pro could not automatically assign a species to were identified manually (Russ 2012). All files which Kaleidoscope Pro classified as ‘Noise’ (195,375 files) were run through the full spectrum software Bat Classify (https://bitbucket.org/chrisscott/batclassify/src/master/). This was to ensure no call sequences within the large number of ‘Noise’ labelled files were missed. Following analysis, 0.5% of labelled files were randomly checked to ensure that the automated identification was reliable (Rowse, Harris & Jones 2018). For all call sequences with >80% certainty in the automated identification, the classification to species was accepted, except for Myotis species where >50% certainty was accepted to ensure call sequences were not excluded from the dataset. These parameters were designed to apply a precautionary approach based on the Precision-Recall metric of the Bat Classify software (https://bitbucket.org/chrisscott/batclassify/src/master/). Statistical analysis All analyses were performed in R statistical software v.4.1.1 (R Core Team 2021) and all statistical tests were considered significant at p < 0.05. We performed generalized linear mixed-effect models (GLMMs) with “glmmTMB” package (Brooks et al. 2017) to assess the effects of PV panels on species-specific bat activity and bat species richness in agricultural landscapes. Echolocation call sequence data were pooled by site and location over the seven-night period, and we defined bat activity as total number of bat call sequences for species or species groups. Due to their low occurrences (<40% of the sites), R. hipposideros and P. nathusii were disregarded for the analysis on species-specific activity. GLMMs on bat species were fitted with a Gaussian distribution (since diagnostic plots were largely unsatisfactory with Poisson or negative binomial distributions) and we applied a squared transformation to the response variable to meet the normality assumption. GLMMs on bat activity were fitted with a negative binomial distribution and we employed zero-inflated models when necessary. We included the presence/absence of PV panels (treatment: solar vs. control site) in interaction with the habitat type surveyed (boundary vs. open field) as explanatory variables while pair IDs were considered as random factors to account for the paired-sampling design. We also included in the models, landscape variables that could potentially affect bat activity in agricultural landscapes, including the proportion of urban, arable land, grassland and broadleaf woodland, and the Euclidean distance to the nearest watercourse. For area-based landscape variables, we considered eight spatial scales (buffers ranging from 250 m to 10 km radii) to qualify local habitats around each site, and to encompass the wide foraging ranges of the bat species studied (Laforge et al. 2021). Landscape variables were derived in QGIS using the Land Cover Map (Environmental Information Data Centre 2019) (20 m resolution) supplied by the Centre of Ecology and Hydrology. When comparing solar PV sites with control sites no statistical differences occurred in the distance to the nearest water source, or in cover of arable land, grassland, broadleaved woodland or urban areas at the different spatial scales with the exception of cover of grassland and arable habitat surrounding the control and solar PV site at the 250m and 500m scales (Appendix S2). To reduce the number of landscape variables and avoid model overparameterisation, we assessed independently the relationships between the response variables and each landscape variable using GLMMs with the same model structure as described above (i.e., including the same random effect and the interaction and using the same distribution family). We compared the second-order Akaike information criterion (AICc) of each model with the model that included the interaction only and retained in the final models only landscape variables at their best scale of effect (Martin 2018) that led to lower AICc (i.e. ΔAICc ≥ 2) (Burnham & Anderson 2002). For highly correlated variables (Spearman coefficient correlation |r| > 0.7), we retained the one leading to lower AICc. From the final full models, we finally ran post hoc pairwise comparisons corrected for multiple testing using the Tukey method in the “lsmeans” package (Lenth 2014). Residual diagnostics were checked with the “DHARMa” package (Hartig 2022). We also checked for multicollinearity, overdispersion, influential outlier, and zero inflation with the “performance” package (Lüdecke et al. 2023).

采样设计
我们在19个太阳能光伏（solar PV）站点采用配对研究设计，以评估田野内及毗邻生境中安装光伏面板的区域，其蝙蝠物种丰富度与活动水平是否高于“空置”匹配对照样地。最终共设置19个太阳能边界生境采样点、19个太阳能开阔生境采样点、19个对照边界生境采样点及19个对照开阔生境采样点。
所有站点均位于英格兰西南部，该区域是英国太阳能光伏站点最集中、蝙蝠物种丰富度最高的区域（Mathews等，2018；英国商业、能源与工业战略部，2021）。涉及私人土地的采样均已获得土地所有者及相关太阳能农场公司的许可。本研究通过声学记录被动监测蝙蝠，无需伦理审批。
对照样地与对应光伏站点处于同一土地管理边界，在样地面积、生境类型、土地利用方式及毗邻生境方面均尽可能匹配。光伏样地与对照样地的平均面积无显著差异（光伏样地均值=59.6公顷，标准差=32.0；对照样地均值=53.2公顷，标准差=28.4；配对t检验：t(18)=1.3，P=0.203）（详见支持信息附录S3）。所有光伏站点均位于放牧或刈割管理的草地，或刈割后的可耕作物田。田间边界对应树篱、林线、林地或植被沟渠，且严格匹配。配对样地间距至少500米且不相邻，以确保在可比景观中获取独立数据（Froidevaux、Louboutin与Jones，2017）。

蝙蝠回声定位录音与物种识别
野外工作于2019年7-10月及2020年同期开展。每个站点连续7晚监测蝙蝠活动（日落前30分钟至日出后30分钟），同时在4个位置开展监测（光伏面板田内的开阔生境与边界生境，以及配对对照样地的对应生境）。使用SM3蝙蝠探测器（SM3 bat detectors，Wildlife Acoustics, Inc.，美国马萨诸塞州梅纳德市）进行录音。所有探测器麦克风（SMM-U2f，频响±6dB 20-100kHz，详见https://www.wildcare.co.uk/amfile/file/download/file/56/product/94208/，Wildlife Acoustics）均通过统一三脚架架设至1.27米高度。探测器设置为8-120kHz自动触发，阈值1-88dB，单段最长录音时长10秒（采样率384kHz）。在对照样地与光伏样地的中心各放置一台探测器，同时在对照与光伏样地的毗邻生境沿线布设探测器。光伏与对照样地内的开阔生境与边界生境对应的探测器间距至少50米。采样在蝙蝠觅食的最优天气条件下进行（无降雨、风速低且气温高于10℃）。录音期间黄昏时的平均气温（±标准差）为16.2±3.1℃（数据来源：https://www.timeanddate.com/weather/）。
声音文件使用零交叉分析软件Kaleidoscope Pro（v.5.4.1，Wildlife Acoustics, Inc.）结合《欧洲蝙蝠分类器（英国版）》（v.5.4.0）进行分析。自动扫描所有10秒录音，识别叫声序列后经人工复核确认物种（Barbastella barbastellus、Eptesicus serotinus、Pipistrellus nathusii、P. pipistrellus、P. pygmaeus、Rhinolophus ferrumequinum及R. hipposideros）或物种类群（Nyctalus spp.、Myotis spp.、Plecotus spp.）。由于难以区分不同物种的回声定位叫声，Myotis属常被归为一类（Russ，2012）。同理，Nyctalus noctula与N. leisleri，以及Plecotus auritus与P. austriacus亦无法始终区分，因此分别归为Nyctalus spp.与Plecotus spp.类群。所有Kaleidoscope Pro无法自动分配物种的文件均经人工识别（Russ，2012）。所有被Kaleidoscope Pro归类为“噪音”的文件（共195375个），均使用全频谱软件Bat Classify（https://bitbucket.org/chrisscott/batclassify/src/master/）重新分析，以确保未遗漏大量“噪音”标签文件中的叫声序列。分析完成后，随机抽取0.5%的已标记文件进行复核，以验证自动识别的可靠性（Rowse、Harris与Jones，2018）。对于自动识别置信度>80%的叫声序列，接受其物种分类结果；但Myotis属物种除外，仅需置信度>50%即可接受，以避免将叫声序列排除出数据集。该参数设置基于Bat Classify软件的精确召回率（Precision-Recall）指标，采用预防性分析原则（https://bitbucket.org/chrisscott/batclassify/src/master/）。

统计分析
所有分析均使用R统计软件v.4.1.1（R核心团队，2021）完成，所有统计检验的显著性阈值设为p<0.05。我们使用“glmmTMB”包（Brooks等，2017）构建广义线性混合模型（generalized linear mixed-effect models, GLMMs），以评估光伏面板对农业景观中物种特异性蝙蝠活动及蝙蝠物种丰富度的影响。将7晚监测期内每个站点及生境类型的回声定位叫声序列数据合并，将蝙蝠活动定义为某一物种或类群的总叫声序列数。由于R. hipposideros与P. nathusii的出现频率较低（<40%的站点），因此在物种特异性活动分析中未纳入这两个物种。针对蝙蝠物种的模型采用高斯分布（因诊断图显示泊松分布或负二项分布拟合效果不佳），并对响应变量进行平方变换以满足正态性假设。针对蝙蝠活动的模型采用负二项分布，必要时使用零膨胀模型。解释变量包括光伏面板的有无（处理组：光伏样地vs对照样地）与调查生境类型（边界生境vs开阔田野）的交互项，同时将配对ID作为随机因子以匹配配对采样设计。模型中还纳入了可能影响农业景观中蝙蝠活动的景观变量，包括城镇、可耕土地、草地及阔叶林地的占比，以及到最近水道的欧氏距离。针对面积类景观变量，我们设置了8个空间尺度（缓冲区半径从250米至10千米），以表征每个站点周边的局部生境，覆盖研究中蝙蝠物种的宽泛觅食范围（Laforge等，2021）。景观变量通过QGIS软件，基于生态与水文学中心提供的《土地覆盖图》（环境信息数据中心，2019，分辨率20米）提取。对比光伏样地与对照样地时，在不同空间尺度下，到最近水源的距离、可耕土地、草地、阔叶林地或城镇用地的占比均无显著差异，仅在250米与500米尺度下，对照与光伏样地周边的草地及可耕生境占比存在差异（附录S2）。为减少景观变量数量并避免模型过度参数化，我们使用与前述相同的模型结构（即包含相同的随机效应、交互项及分布族），通过广义线性混合模型分别评估响应变量与每个景观变量的关系。将每个模型的二阶赤池信息准则（AICc）与仅包含交互项的模型进行比较，最终模型仅保留能降低AICc（即ΔAICc≥2）的最佳尺度下的景观变量（Martin，2018；Burnham与Anderson，2002）。对于高度相关的变量（斯皮尔曼相关系数|r|>0.7），保留能使AICc更低的变量。从最终的全模型出发，我们使用“lsmeans”包（Lenth，2014）中的Tukey法进行经多重检验校正的事后两两比较。使用“DHARMa”包（Hartig，2022）进行残差诊断。同时使用“performance”包（Lüdecke等，2023）检验多重共线性、过度离散、有影响的异常值及零膨胀情况。