Chironomid-climate continentality conundrum

PLOS One | https://doi.org/10.1371/journal.pone.0327780 August 5, 2025 6 / 18 (analysed by Varvara Bakumenko), and Swedish (analysed by Simon Belle) parts were identified using keys by Klink and Pillot [54], Brooks et al. [53], Larocque-Tobler [55], and Andersen et al. [56]. On average, 69 chironomid head cap- sules were identified per sample, with a range of 47–139 (S1 File).Final taxonomic harmonisation was done after merg- ing all the above-described datasets. All identifications at genus or subfamily taxonomic level (Tanytarsini, Tanytarsus spp., Paratanytarsus spp., Tanypodinae, Chironomini, Orthocladinae) were excluded from the merged dataset to avoid including broad groups of Chironomidae species with a wide range of ecological preferences. The excluded taxa made up less than 10% of the total head capsules count in the dataset and no more than 8% in any individual sample. Cricoto- pus intersectus- type was merged with Cricotopus laricomalis-type into one type due to the likelihood of misidentification of these morphotypes. Morphotype-level identifications from Einfeldia (Einfeldia dissidens-type), Zalutschia (Zalutschia zalutschicola-type), Eukieferiella (Eukiefferiella coerulescens-type), and Dicrotendipes (Dicrotendipes nervosus-type and Dicrotendipes notatus-type) were merged into corresponding genera level groups due to the differences in identification resolution in parts of the dataset. Harmonisation was done before transforming the data into relative abundances. Numerical analysis

The harmonized chironomid count data was transformed into relative abundances, and thereafter square-root trans- formed. To remove rare taxa, only morphotypes with an abundance higher than 2% in at least one sample were included for numerical analysis to improve performance of the inference model [57]. Based on KOI, the dataset was divided into 3 parts: continental with KOI −10–0, transitional with KOI 0–10, and oceanic with KOI 10–20 [2]. Analysis of similarities (ANOSIM; [58]) was applied to justify the KOI-based division of the chironomid assemblages. Bedrock data were grouped in 4 groups (sand-containing, clay-containing, carbonates-containing, granits/granitoids) and coded as numbers (1–4) for the numerical analyses. Principal component analysis (PCA) was applied to the environmental and climate data of the dataset to investigate their gradients length. The Shapiro-Wilcox test and Spearman correlations were used to test for collinearity between envi- ronmental and climatic variables. Variables with correlation coefficient > ±0.7 were considered highly correlated and their effect on chironomid assemblages could not be distinguished from each other. Detrended correspondence analysis (DCA; [59]) was applied to the chironomid assemblages data to examine the distribution of taxa and the compositional gradient lengths along the first two DCA axes. Redundancy analysis (RDA) was chosen based on the length of DCA axis-1 and axis-2 of the dataset (2.9 and 2 SD units, respectively; [60,61]). RDA was

applied to determine which environmental variables explain significant compositional variation in the chironomid data. Weighted averaging-partial least squares (WA-PLS; [62]) was performed to evaluate the idea of developing a chironomid- based training set applicable to the reconstructions of continentality. The continentality related variable which showed the strongest relationship to the chironomid assemblages and had a λ1:λ2 ratio of more than 1 in RDA was used. The stron- gest transfer function was determined as the one producing the lowest cross-validated root mean square error of pre- diction (RMSEP). The relevant components were accepted as statistically significant at the p ≤ 0.05 level. Bootstrapping techniques (9999 permutations; [60,63]) were used to estimate cross-validated error and performance statistics for the WA-PLS inference model, such as RMSEP, maximum and mean bias, and the coefficient of determination (R2) between inferred and predicted values. Indicator species analysis (INDVAL; [64]) was applied to reveal characteristic morphotypes for the best performing con- tinentality related variable according to the RDA results. Weighted-average regression with inverse deshrinking ([62]) was applied to the taxa that were determined as potential indicators to estimate taxon-specific continentality optima. The software program R version 4.1.1. (R Core Team, 2021) was used to perform numerical analyses and create plots. The following packages were used: ‘tidyverse’ for data visualisation [65], ‘dplyr’ for data restructuring and basic calcula-

tions [66], ‘vegan’ for ordination and ANOSIM [67], ‘rioja’ for WA-PLS and WA and plotting the stratigraphic diagram [68], and ‘indicspecies’ for performing the INDVAL [69].