Fluttering in a changing world: Effects of urbanization and nectar plants on butterfly movement patterns

We aimed to answer the general question whether urbanization affects butterfly movement patterns and more specific, whether (1) the mobility of the investigated small white (Pieris rapae) and the small heath butterfly (Coenonympha pamphilus) is affected differently and (2) these butterflies show altered tortuosity-patterns along a rural-urban gradient. The study sites were situated along a rural-urban gradient in the Berlin-Brandenburg metropolitan region (Germany). We recorded GPS-movement trajectories of two common butterfly species differing in territoriality, agility and habitat requirements. Movement trajectories were analyzed in terms of mobility (flight speed and time investment in stopping, resting and nectaring) and tortuosity measures and the effect of urbanization on the derived variables was investigated.

Methods
The fieldwork was carried out from June to September 2020, between 9:00 and 15:00 h on sunny and calm days. To gain spatial data of butterfly movements, we tracked two to five butterflies per species and study site if the species was present (P. rapae: 3,7±0,9 and C. pamphilus: 3,4±0,8 mean±SD of butterflies per study site), by constantly keeping a distance of about 2.5±0.5 m between observer and animal. Movements of the observer were logged by a DGPS-Receiver with real-time differential correction (“Trimble R10”), mounted on a GPS-backpack and operated via handheld control. The GPS-position was recorded with a temporal resolution of one second and an average horizontal accuracy of 0.02 m.

The high accuracy of the DGPS-Receiver allowed for the tracking of both, large scale fast flight, as well as small scale fluttering flight. In order to be able to identify stops, we noted the respective second on a prepared data sheet every time the observed butterfly was landing. Additionally, we noted the respective action of the butterfly during a stop, being either nectaring, resting, basking or oviposition. The recording stopped, when the observer lost sight of the butterfly or after a maximum of 12 minutes (00:05:36 ±00:02:38; mean ±SD). Since the high volatility of the butterfly movements often made it impossible to write all information down immediately, comments regarding the tracks were recorded as audio files via “AGPTEK” Lavalier microphone on a smartphone. The audio-recording started simultaneously with the GPS-recording in order to achieve initial temporal synchronization. This allowed for a subsequent completion of the data sheets and facilitated the identification of stops or data gaps (the latter sometimes occurred when butterflies were flying through groves and the GPS-signal was insufficient to locate the observers position).

The meteorological parameters wind speed (m/s), air temperature (°C) and relative humidity (%) were measured with a Thermo-Hygro-Anemometer („PCE-THA 10“) approximately 2 m above ground subsequently to each GPS-tracking. 

To attain an approximation for food availability, a transect was drawn through the respective movement area and the coverage of flowering nectar plants was estimated for each 1-m segment of the transect within a margin of 0.5 m to the left and right side respectively. The mean coverage of nectar plants was then calculated for each transect. 

The recorded movement patterns were exported as point features in an ESRI shapefile and imported into ArcGIS 10.3.1 (Esri Inc. 2018) for processing. Line geometries were drawn along the subsequent GPS-positions, using the time information (seconds since start of the track) as guidance. Due to small movements of the observer, even when trying to stand perfectly still, the butterfly landings were not recorded as perfect stops but rather as point clouds. For each point cloud, previously confirmed as a stop via comparison with the data sheet information, the median center was calculated with the ArcGIS tool “Median Center”. The computed coordinate was then attributed to all affected GPS-points and used for further analyses. For the calculation of the tortuosity however, all but the first (adjusted) point of each stop were excluded to avoid problems with the computation of the mean cosine. Stops during initiation and termination of the GPS-logging were excluded from all analyses.

To describe the butterfly movements, we derived several parameters from the spatial GPS-data and assigned them to one of the categories (1) mobility or (2) tortuosity.

As measures for butterfly mobility, the mean flight speed (“flight speed” in m/s) and the share of time spent stopping (“stopping time”), nectaring (“nectaring time”) and resting (“resting time”) were calculated. Due to low sample sizes, the time spend ovipositing and the time spent basking were not analyzed.

As tortuosity measure, the sinuosity of the flight path, as defined by Benhamou (2004) was calculated with the R package trajr and the function TrajSinuosity2 (McLean & Skowron Volponi 2018). High values of sinuosity indicate a high tortuosity, while low values indicate a virtually straight flight path. In cases of split tracks (due to data gaps), the sinuosity was computed for each segment individually and the mean sinuosity, weighted by relative track segment length was then calculated for the entire butterfly track. 

As a measure for urbanization, the percentage of sealed soil surface around each study site was calculated. Within the boundaries of Berlin, this calculation was based on the Berlin impervious soil coverage map (BISCM, SenUDH 2016). Since such data was not available for Brandenburg, approximations based on the BISCM (SenUDH 2016) and the Berlin biotope type map (SenUDH 2014) were computed. For this purpose, the average percentage of surface sealing was calculated for each biotope type via the Zonal Statistics tool in QGIS 2.18.11 (QGIS-DT 2018) based on the available rasterized data for Berlin. The calculated values were then attributed to the biotope patches of the Brandenburg biotope type map (LfU 2009). The vector data gained was merged with the BISCM and then transformed into raster data (resolution: 2 x 2 m). This map was used to calculate the mean percentage of sealed surface within a 500, 1000 and 2000 m radius around the center point of each study site via the Zonal Statistics tool in QGIS 3.4.12 (QGIS-DT 2018).

To attain an approximation for the size of the habitat patch, the area (ha) of the habitat site was calculated using the “calculate geometry” tool of ArcGIS Version 10.3.1 (Esri Inc. 2018). To do so, we first digitalized the continuous grassland site (not interrupted by urban matrix, forest, agricultural crops, waterbodies, ornamental lawns or other flower-free biotopes as reed beds) based on orthophotos of the Berlin-Brandenburg area in summer 2020 and winter 2021 (Geoportal Berlin 2020, 2021; LGB 2020) and supported by the Berlin and Brandenburg biotope type map (LfU 2009, SenUDH 2014).

References: 

Benhamou, S. (2004). How to reliably estimate the tortuosity of an animal’s path: straightness, sinuosity, or fractal dimension? Journal of Theoretical Biology, 229, 209–220. https://doi.org/10.1016/j.jtbi.2004.03.016

Esri Inc. ArcMap™. Version 10.3.1.4959. 18 June 2015 [27.04.2018]. Available online: https://www.esri.de/de-de/home.

Geoportal Berlin (2020). Digitale farbige Orthophotos 2020 (DOP20RGB). Datenlizenz Deutschland – Namensnennung – Version 2.0. URL: https://fbinter.stadt-berlin.de/fb/wms/senstadt/k_luftbild2020_rgb [19.04.2021].

Geoportal Berlin (2021). Digitale farbige Orthophotos 2021 (DOP20RGBI). Datenlizenz Deutschland – Namensnennung – Version 2.0. URL: https://fbinter.stadt-berlin.de/fb/wms/senstadt/k_luftbild2021_rgb [19.04.2021].

LfU (Landesamt für Umwelt Brandenburg) (2009). CIR-Biotoptypen 2009 – Flächendeckende Biotop- und Landnutzungskartierung im Land Brandenburg (BTLN). dl-de/by-2-0, URL: https://geoportal.brandenburg.de/detailansichtdienst/render?view=gdibb&url=https://geoportal.brandenburg.de/gs-json/xml?fileid=B57B9F35-AFFF-49F2-BA32-618D1A1CD412 [19.04.2021].

LGB (Landesvermessung und Geobasisinformation Brandenburg) (2020). Digitale Orthophotos 20cm Bodenauflösung Farbe Brandenburg mit Berlin (WMS). GeoBasis-DE/LGB, dl-de/by-2-0, URL: https://geobroker.geobasis-bb.de/gbss.php?MODE=GetProductInformation&PRODUCTID=253b7d3d-6b42-47dc-b127-682de078b7ae [19.04.2021].

McLean, D.J., Skowron Volponi, M.A. (2018). trajr: An R package for characterisation of animal trajectories. Ethology, 124, 440–448. https://doi.org/10.1111/eth.12739

QGIS-DT (QGIS Development Team) (2018). “QGIS Geographic Information System: Open Source Geospatial Foundation.” http://qgis.osgeo.org.

SenUDH (Senate Department for Urban Development and Housing) (2014). Geoportal Berlin. 05.08 Biotope Types. dl-de/by-2-0, URL: https://www.berlin.de/umweltatlas/en/biotopes/biotope-types/2013/maps/artikel.970219.en.php [02.03.2020].

本研究旨在解答两个核心科学问题：其一，城市化是否会对蝶类的运动模式产生影响；其二，更为具体地，（1）供试的菜粉蝶（Pieris rapae）与小赭弄蝶（Coenonympha pamphilus）的运动活动性是否会受到城市化的差异化影响，（2）这两种蝶类的运动曲折模式是否会沿城乡梯度发生改变。本研究的样地设置于德国柏林-勃兰登堡大都市区的城乡梯度带上。本研究记录了两种在领地性、活动敏捷性与栖息地需求上存在显著差异的常见蝶类的GPS运动轨迹，并从运动活动性（飞行速度、停驻、休憩与取食花蜜的时间投入）以及轨迹曲折度指标两个维度对轨迹展开分析，同时探究城市化对衍生变量的影响。

研究方法
野外调查于2020年6月至9月期间开展，观测时段为每日9:00至15:00，且仅选择晴朗无风的天气。为获取蝶类运动的空间数据，本研究在样地中出现目标蝶类时，对每个物种追踪2~5只个体（每个样地的菜粉蝶追踪数量为3.7±0.9只，小赭弄蝶为3.4±0.8只，均为平均值±标准差），观测者与蝶类之间始终保持2.5±0.5米的距离。观测者的运动轨迹由实时差分GPS接收机（DGPS-Receiver，型号"Trimble R10"）记录，该设备搭载于GPS背包中，通过手持控制器进行操作。GPS定位的时间分辨率为1秒，平均水平定位精度达0.02米。
该DGPS接收机的高精度特性，既可以记录大范围的快速飞行轨迹，也能捕捉小范围的振翅飞行细节。为准确识别蝶类的停驻行为，每当观测到蝶类降落时，研究人员会在预先准备的数据记录表上标注对应的时刻。同时记录停驻期间蝶类的具体行为，包括取食花蜜、休憩、晒翅调温与产卵。当观测者无法再看到目标蝶类，或单次追踪时长达到12分钟上限时，本次追踪结束（单次追踪时长为5分36秒±2分38秒，平均值±标准差）。
由于蝶类运动极具机动性，研究人员往往无法立刻记录全部观测信息，因此通过智能手机搭配"AGPTEK"领夹式麦克风，将轨迹相关的观测记录存储为音频文件。音频录制与GPS记录同步启动，以实现初始的时间同步。这一设置便于后续补全数据记录表，同时也能更便捷地识别停驻行为或数据缺失段（数据缺失常发生在蝶类穿越树林时，此时GPS信号不足以定位观测者位置）。
在每次GPS追踪结束后，于距地面约2米高度处，使用温湿风速仪（Thermo-Hygro-Anemometer，型号"PCE-THA 10"）测量风速（单位：m/s）、气温（单位：℃）与相对湿度（单位：%）等气象参数。
为近似评估食物可获得性，研究人员在每个样地的蝶类运动区域内设置样线，对样线每1米区段左右各0.5米范围内的开花蜜源植物覆盖度进行估算，最终计算每条样线的蜜源植物平均覆盖度。

记录的运动轨迹数据以点要素形式导出为ESRI形状文件（ESRI shapefile），并导入ArcGIS 10.3.1（Esri Inc. 2018）进行预处理。根据GPS记录的时间信息（追踪启动后的秒数），将连续的GPS点位连接为线要素几何图形。由于观测者即便试图保持静止也会存在微小位移，蝶类的降落停驻并未被记录为精确的单点，而是形成点云。对于每一处经数据记录表确认的停驻点云，使用ArcGIS的"中位数中心（Median Center）"工具计算其中位中心坐标，将该坐标赋予该停驻对应的所有GPS点位，用于后续分析。但在计算轨迹曲折度时，为避免平均余弦值计算出现问题，需剔除每个停驻段中除首个（已校正）点位之外的其余点位。此外，GPS记录启动与结束阶段的停驻行为，均被排除在所有分析之外。

为描述蝶类的运动特征，研究人员从空间GPS数据中提取了多项参数，并将其划分为（1）运动活动性与（2）轨迹曲折度两大类。
作为蝶类运动活动性的衡量指标，本研究计算了平均飞行速度（单位：m/s，记为"flight speed"）、停驻时间占比（记为"stopping time"）、取食花蜜时间占比（记为"nectaring time"）与休憩时间占比（记为"resting time"）。由于产卵与晒翅调温的行为样本量过小，未对这两种行为的时间投入展开分析。
作为轨迹曲折度的衡量指标，本研究采用Benhamou（2004）定义的飞行路径弯曲度（sinuosity），通过R语言包trajr的TrajSinuosity2函数进行计算（McLean & Skowron Volponi 2018）。弯曲度数值越高，代表轨迹曲折度越大；数值越低，则代表飞行路径越趋近于直线。当轨迹因数据缺失出现分段时，先分别计算各分段的弯曲度，再以各分段的相对长度为权重，计算整条蝶类运动轨迹的平均弯曲度。
本研究以每个样地周边的地表硬化率作为城市化程度的衡量指标。在柏林市域范围内，该指标基于柏林不透水地表覆盖图（Berlin impervious soil coverage map，简称BISCM，SenUDH 2016）计算。由于勃兰登堡州暂无此类数据，研究人员基于柏林不透水地表覆盖图（SenUDH 2016）与柏林生境类型图（SenUDH 2014）进行近似计算。具体步骤为：基于柏林已有的栅格数据，通过QGIS 2.18.11的分区统计工具（Zonal Statistics）计算每种生境类型的平均地表硬化率，并将该数值赋予勃兰登堡生境类型图（LfU 2009）中的对应生境斑块。将得到的矢量数据与柏林不透水地表覆盖图合并后，转换为分辨率为2×2米的栅格数据。随后，通过QGIS 3.4.12的分区统计工具，计算每个样地中心点周边500米、1000米与2000米半径范围内的平均地表硬化率。

为近似估算栖息地斑块的面积，研究人员使用ArcGIS 10.3.1的"计算几何"工具计算样地栖息地的面积（单位：公顷）。具体步骤为：基于2020年夏季与2021年冬季的柏林-勃兰登堡地区正射影像（Geoportal Berlin 2020, 2021; LGB 2020），并结合柏林与勃兰登堡生境类型图（LfU 2009, SenUDH 2014），对连续的草地生境进行数字化描边——该类生境未被城市建成区、林地、农作物、水体、装饰性草坪或芦苇床等无花植物生境所割裂。

参考文献：
Benhamou, S. (2004). How to reliably estimate the tortuosity of an animal’s path: straightness, sinuosity, or fractal dimension? Journal of Theoretical Biology, 229, 209–220. https://doi.org/10.1016/j.jtbi.2004.03.016
Esri Inc. ArcMap™. Version 10.3.1.4959. 18 June 2015 [27.04.2018]. Available online: https://www.esri.de/de-de/home.
Geoportal Berlin (2020). Digitale farbige Orthophotos 2020 (DOP20RGB). Datenlizenz Deutschland – Namensnennung – Version 2.0. URL: https://fbinter.stadt-berlin.de/fb/wms/senstadt/k_luftbild2020_rgb [19.04.2021].
Geoportal Berlin (2021). Digitale farbige Orthophotos 2021 (DOP20RGBI). Datenlizenz Deutschland – Namensnennung – Version 2.0. URL: https://fbinter.stadt-berlin.de/fb/wms/senstadt/k_luftbild2021_rgb [19.04.2021].
LfU (Landesamt für Umwelt Brandenburg) (2009). CIR-Biotoptypen 2009 – Flächendeckende Biotop- und Landnutzungskartierung im Land Brandenburg (BTLN). dl-de/by-2-0, URL: https://geoportal.brandenburg.de/detailansichtdienst/render?view=gdibb&url=https://geoportal.brandenburg.de/gs-json/xml?fileid=B57B9F35-AFFF-49F2-BA32-618D1A1CD412 [19.04.2021].
LGB (Landesvermessung und Geobasisinformation Brandenburg) (2020). Digitale Orthophotos 20cm Bodenauflösung Farbe Brandenburg mit Berlin (WMS). GeoBasis-DE/LGB, dl-de/by-2-0, URL: https://geobroker.geobasis-bb.de/gbss.php?MODE=GetProductInformation&PRODUCTID=253b7d3d-6b42-47dc-b127-682de078b7ae [19.04.2021].
McLean, D.J., Skowron Volponi, M.A. (2018). trajr: An R package for characterisation of animal trajectories. Ethology, 124, 440–448. https://doi.org/10.1111/eth.12739
QGIS-DT (QGIS Development Team) (2018). "QGIS Geographic Information System: Open Source Geospatial Foundation." http://qgis.osgeo.org.
SenUDH (Senate Department for Urban Development and Housing) (2014). Geoportal Berlin. 05.08 Biotope Types. dl-de/by-2-0, URL: https://www.berlin.de/umweltatlas/en/biotopes/biotope-types/2013/maps/artikel.970219.en.php [02.03.2020].