Contrasting response of microeukaryotic and bacterial communities to the interplay of seasonality and local stressors in shallow soda lakes

Abstract Seasonal environmental variation is a leading driver of microbial planktonic community assembly and interactions. However, departures from usual seasonal trends are often reported. To understand the role of local stressors in modifying seasonal succession, we sampled fortnightly, throughout three seasons, five nearby shallow soda lakes exposed to identical seasonal and meteorological changes. We characterised their microeukaryotic and bacterial communities by amplicon sequencing of the 16S and 18S rRNA gene, respectively. Biological interactions were inferred by analyses of synchronous and time-shifted interaction networks, and the keystone taxa of the communities were topologically identified. The lakes showed similar succession patterns during the study period with spring being characterised by the relevance of trophic interactions and a certain level of community stability followed by a more dynamic and variable summer-autumn period. Adaptation to general seasonal changes happened through shared core microbiome of the lakes. Stochastic events such as desiccation disrupted common network attributes and introduced shifts from the prevalent seasonal trajectory. Our results demonstrated that, despite being extreme and highly variable habitats, shallow soda lakes exhibit certain similarities in the seasonality of their planktonic communities, yet local stressors such as droughts instigate deviations from prevalent trends to a greater extent for microeukaryotic than for bacterial communities.


Introduction
Seasonal changes of environmental parameters are the primary dri vers of ann ual succession d ynamics of planktonic communities instigating c har acteristic ecological pr ocesses r epeated eac h year (Sommer et al. 1986(Sommer et al. , 2012, making seasonality the focal point of numerous aquatic microbial ecology studies (Bista et al. 2017, Lambert et al. 2018, Reji et al. 2020. Departures from usual seasonal patterns are typically explained by interannual climatic variations or long-term trends (Fuhrman et al. 2015 ). Seasonal changes affect shallow lakes because they have a high surfaceto-volume ratio (Jeppesen et al. 2009, Cobbaert et al. 2014, Li et al. 2021. Endorheic shallow lakes are even more impacted by external envir onmental c hanges due to their intrinsic lack of outflow and consequent importance of e v a por ation (Fr ondini et al. 2019, Boros et al. 2020, Wang et al. 2022, making them ideal systems for the study of the seasonality of pela gic comm unities. Indications of interannual variation in the plankton dynamics of endorheic aquatic systems have been shown before (García et al. 1997, Márton et al. 2023. Ho w e v er, it is not known whether such variations arise because of changes of regional magnitude such as y ear-to-y ear w eather and/or climate variability, or if they are the consequence of stress events of only local influence such as algae blooms or complete dry ups. Soda lakes are a form of saline endorheic lakes c har acterised by carbonate, bicarbonate and sodium as dominant ions (Boros and K olpako va 2018 ). T hese systems represent the most alkaline aquatic habitats on the surface of our planet (Jones and Grant 2000, Sorokin et al. 2011. The Kiskunság National Park (Hungary), in the central part of the Carpathian Basin (Central Eur ope), is a r egion with a particularl y high density of soda pans (i.e. shallow soda lakes) (Boros et al. 2014 ). The climate of this region is temperate continental (K o vács and Jakab 2021 ), which cr eates consider able seasonal temper atur e and water le v el fluctuations . T he drought prolonging and precipitation pattern modifying the effect of climate change (Konapala et al. 2020 ) is also substantial, resulting, among others, in an incr easing fr equency of desiccation e v ents (Bor os et al. 2013 ). It has been implicated that climate change-induced weather anomalies also exacerbate the inter annual v ariation in the seasonality of the biogeochemistry of individual pans (Boros et al. 2020 ). Nevertheless, while the seasonal envir onmental c hanges ar e becoming less pr edictable, the close proximity of these pans means that they are still exposed to pr acticall y identical climatic and weather conditions . T his makes the soda pans of Kiskunság a perfect system with which to disentangle the impact of local stress events on planktonic communities from the effect of common seasonal environmental changes.
Micr obial inter actions play a crucial r ole in micr obial functioning (Tshikantwa et al. 2018 ) and comm unity assembl y (Leibold et al. 2022 ), highlighting the need to consider species interactions alongside compositional changes when studying the seasonality of aquatic microbial communities (Gilbert et al. 2011 ). Networkbased a ppr oac hes allow insight into inter actions of e v en highl y complex communities such as those of microbes (Leibold et al. 2022 ). Sync hr onous networks can assess the simultaneous occurrence and abundance of organisms, while time-shifted association networks can elucidate time-delayed impacts of species interactions (Faust et al. 2015, Fuhrman et al. 2015, Röttjers and Faust 2018. Network analysis can also identify k e ystone species (Fuhrman et al. 2015, Banerjee et al. 2018, Röttjers and Faust 2018, whic h ar e usuall y defined as hubs (i.e. highl y associated species in the network) that hav e dispr oportionatel y high importance in the comm unity r elativ e to their abundance (Jordán 2009, Berry and Widder 2014, Banerjee et al. 2018. Although micr obial comm unities typicall y consist of hundr eds or thousands of species (Bengtsson-P alme 2020 ), onl y a smaller fraction of the microbial species are shared among communities of a specific habitat type (e.g. soda pans) and can be defined as its cor e micr obiome (Shu et al. 2020 ). The cor e micr obiome is hypothesised to r epr esent the functionall y or ecologicall y most important taxa of a habitat type (Degenhardt et al. 2020, Neu et al. 2021. It has also been suggested that microbial adaptation to envir onmental c hanges thr ough dispersal mediated species sorting also oper ates primaril y on the cor e micr obiome instead of the r ar e or sporadic (i.e. non-core) community members (Niño-García et al. 2016 ).
Soda lakes and pans ar e c har acterised by abundant phytoplankton (Afonina and Tashl yk ov a 2020 , Somogyi et al. 2022 ) and zooplankton (Burian et al. 2013, Horváth et al. 2013, which feed large amounts of migratory birds (Wagaw et al. 2019, Boros et al. 2023, while fish are typically absent (Felföldi 2020 ). Studies describing the planktonic microbiome of soda pans in the Carpathian Basin (Felföldi 2020, Szabó et al. 2020, Márton et al. 2023, Canada (Zorz et al. 2019 ), Brazil (Cotta et al. 2022 ), China (Zhao et al. 2020 ) and East Africa (Sc ha gerl 2016 ) have identified unique bacterial and eukaryotic communities . T he seasonal changes of the composition (Szabó et al. 2020 ) and function (Eiler et al. 2003 ) of bacterial communities of soda pans have also been described before, but the impact of seasonality on microbial interactions remains unexplored. Moreover, the seasonal succession mechanisms of bacterial and microeukaryotic communities of soda pans have never been compared. Microeukaryotes and bacteria have various inherent biological differences that might influence their seasonal dynamics, such as their different cell size and structur e, gener ation times and life history tr aits. Accordingl y, pr e vious studies hav e demonstr ated higher turnover for microeukaryotes than bacteria (Zhang et al. 2017 ) and that drift and dispersal play a more important role in the structuring of their communities (Logares et al. 2018, Vass et al. 2020. Meanwhile, bacterial communities tend to have higher ada ptability to envir onmental fluctuations than eukaryotes  ).
We carried out an extensive field study e v aluating the seasonal dynamics of planktonic microbial communities of five soda pans located in the same region (i.e. Kiskunság NP) throughout three seasons (spring, summer, autumn) by fortnightly sampling. We aimed to understand the impact of seasonality on the structure and interactions of bacterial and microeukaryotic communities and the contribution of core and non-core microbiome to the adaptation to en vironmental changes . Our primary hypothesis was that common habitat type (i.e. shallow soda lake) and identical weather and climatic exposure due to the proximity and sync hr onous sampling of the five sites results in analogous seasonal succession patterns, both in respect of community composition and interactions. We further hypothesised that adaptation to seasonal c hanges ha ppens primaril y thr ough species r ecruitment fr om the cor e r ather than the non-cor e comm unity, while noncor e taxa ar e to a lar ge extent involv ed in the r esponse to local str essors. Finall y, we hypothesise that bacterial communities are less affected by local stress events than eukaryotic communities.

Sample collection and environmental parameters
We selected five soda pans located within an area of 14 km 2 of Kiskunság National Park, Hungary: Böddi-szék, Kelemen-szék, Sós-ér, Zab-szék and an unnamed pan (Pan no. 60) (Fig. 1 ). Two different subtypes of soda pan are distinguished in this region: the pr e v ailing turbid-white subtype (r eferr ed to as turbid from here on), whic h is c har acterised by high amounts of suspended inorganic cla y particles , and the less common non-turbid, brown subtype (r eferr ed to her e as br own), whic h has a c har acteristic br own colour due to high organic material content and a low amount of suspended inorganic material (Boros et al. 2014 ). Among the sampled pans Sós-ér was brown, while the other pans were turbid.
The pans were sampled fortnightly through three seasons: spring (sampling time 1-4), summer (sampling time 5-10) and autumn (sampling time 11-14) from 12 April to 14 November in 2017 to cover the main productivity period of the year and to ensure icefree conditions (Fig. 2 A). All samples were pooled by proportionally mixing water collected from 1 cm below the surface, from at least five different points near the deepest part of the pans. For micr obial comm unity anal yses, w ater samples w ere collected into a 1-litre sterile bottle after filtering through a 40 μm mesh size plankton net to r emov e lar ge zooplankton. Lar ger algae ar e pr acticall y absent fr om these pans; the mean contribution of picosized (cell diameter < 2 μm) algae to total phytoplankton biomass is ∼85% in the turbid pans (Somogyi et al. 2022 ), while in the brown soda pans nanoplanktonic (2-20 μm) algae have a similar contribution ( ∼80%), with a moderate amount of picophytoplankton ( ∼8%) (Szabó et al. 2020, Somogyi et al. 2022. The water samples wer e tr ansferr ed to the labor atory in a cooling box, wher e planktonic microbial cells were collected onto 0.1 μm pore size filters (MF-Millipor e membr ane filter) b y filtering 30 ml w ater for turbid pans and 50 ml for the brown pan. The filters were stored at -80 • C until further processing.
Temper atur e, pH (SenTix 41 electr ode), conductivity (Tetr aCon 325 cell) and dissolved O 2 (CellOx-325 electrode) were measured on site with a MultiLine Handheld Meter model 340i (WTW, Weilheim in Oberbayern, German y). Meteor ological data fr om the Soltszentimre station (within 15 km distance from all sampling sites) wer e pr ovided by the Hungarian Meteorological Service. Further envir onmental par ameters wer e determined in the labor atory. Chlorophyll a was measured based on a pr e viousl y published   Numbers r epr esent the sampling times, differ ent symbol sha pes the fiv e pans, while differ ent colours r epr esent the thr ee studied seasons (gr een for spring, y ello w for summer and red for autumn). study (Mentes et al. 2018 ), while total phosphorus (TP), total nitrogen (TN) and dissolved organic carbon (DOC) concentrations were determined as described in Nydahl et al. ( 2019 ). Electrical conductivity data (EC) were converted to salinity based on a previously published equation estimated specifically for soda lakes of the same region (Total ions (g/L) = 0.792 * EC (mS/cm) + 179) (Boros et al. 2014 ). A total of 20 litres of water collected at different points around the deepest part of the open water area was sie v ed thr ough a 40 μm mesh size plankton net and pr eserv ed in 70% ethanol until the abundance and composition of zooplankton was determined . Subsampling was used to reduce the counting effort in the case of samples containing more than 300 specimens, otherwise all the specimens were counted. The abundance of cladocera and copepoda was calculated as individuals per litre.

Community analysis
DN A w as extracted using DNeasy Po w erSoil Kit (QIAGEN) accor ding to the manufactur er's instructions. Micr oeukaryotic and bacterial community composition was determined by amplification and sequencing of the V4-V5 region of the 18S rRNA gene and V3-V4 region of the 16S rRNA gene, respectively. No archaea-specific PCRs were carried out, because they r epr esent onl y a minor fraction of the pr okaryotic comm unity in these sites (Korponai et al. 2019, Szabó et al. 2020. Amplification and libr ary pr epar ation for amplicon sequencing were carried out based on pr e viousl y published studies (Székely et al. 2019, Vass et al. 2020 ) ( Text S1 .), while sequencing was performed at the Swedish National Genomics Infr astructur e (Uppsala, Sweden) on Illumina MiSeq platform (Illumina Inc, San Diego, CA, USA) in a 2 × 300 bp paired-end format and with v3 chemistry.
Sequence r ead pr ocessing, alignment and taxonomic assignments were carried out using mothur v. 1.41.1 (Schloss et al. 2009 ). Operational taxonomic units (OTUs) were assigned at 99% similarity cutoff and r ar efied OTU sets wer e cr eated as a basis for subsequent analyses. For the amplified region and based on previous studies of the surveyed habitat, this cutoff value was found to be suitable to avoid most diversity estimation biases caused by clustering different species to the same OTU or splitting organisms having multiple rRNA copies to separate clusters (Johnson et al. 2019, Schloss 2021. The seventh sampling time of Pan no. 60 was discarded from the 18S rRNA gene amplicon dataset due to low sequencing quality. Reads were subsampled to the read number of the samples with the lo w est sequence counts (62 samples in the 18S rRNA amplicon set, n = 2407; and 63 samples in the 16S rRNA amplicon set, n = 3188). A detailed description of the sequence analysis is provided in Supplementary Text S2 .
OTUs present in all five studied soda pans were defined as core5 and OTUs shared between the four turbid soda pans as core4.
OTUs not shared between the pans were defined as non-core5 and non-cor e4, r espectiv el y.

Sta tistical anal yses
T he measured en vironmental variables were scaled to unit variance and compared using principal component analysis (PCA). The communities of the turbid and brown pans were compared b y one-w ay perm utational m ultiv ariate anal ysis of v ariance (PER-MANOVA, 'adonis' function, permutations = 999) based on Bray-Curtis (BC) dissimilarity of the microeukaryotic O TUs (eO TUs) and bacterial OTUs (bOTUs), while differences in community compositions among seasons and pan identity were tested by two-way PERMANOVA (permutations = 999). Non-metric multidimensional scaling (NMDS) based on BC distance was used to visualise the microeukaryotic and bacterial communities, while the 'envfit' function was applied to plot significantly fitted ( P < 0.05) environmental vectors onto the NMDS or dinations. Mantel-tests w ere implemented to verify the results of 'envfit' by identifying significant correlation between microeukaryotic, and bacterial communities and envir onmental v ariables, and zooplankton abundance. Highl y collinear v ariables wer e identified based on Pearson correlation ( | r | > 0.7) (Dormann et al. 2013 ). To assess the structural temporal dynamics of each pan, BC dissimilarity between sampling occasions (i.e. time distance decay curve) was calculated for both 18S and 16S rRNA gene O TUs (i.e . eO TUs and bO TUs , r espectiv el y) and the BC dissimilarity between consecutive sampling occasions was considered as a proxy of community turno ver. T he effect of dr ought on comm unity turnov er of micr oeukaryotes and bacteria w as tested b y comparing the BC dissimilarities between consecutive sampling times during the drought period in late summer and autumn (i.e. between sampling times 7 and 14) for the three nondrying pans (Böddi-szék, Pan no. 60 and Sós-ér) and the two drying pans that underwent various drought events during this period (Zab-szék and K elemen-szék). Anal ysis of variance (ANOVA) with Tuk e y's post-hoc test was used to test the differences of environmental v ariables, cor e O TU contribution and turno ver between soda pans and seasons.

Netw ork anal yses
The extended Local Similarity Analysis (eLSA) is a robust time series analysis tool that not only detects microbial associations that ar e pr esent thr oughout the study period (i.e. global corr elations), but also those that only occur in a subinterval of the time series (i.e. local associations). Furthermore, eLSA considers both associations between taxa that coexist in time (i.e. co-occurence) and time-shifted or time-lagged correlations (Ruan et al. 2006, Xia et al. 2011, Fuhrman et al. 2015. To better understand sync hr onous and async hr onous inter actions, we gener ated tw o eLSA netw orks of the microeukaryotic and bacterial communities of each soda pan: one using the sync hr onous corr elations (i.e. co-occurr ence networks , dela y 0) and another using only time-shifted correlations (delay 1 or -1). The eLSA (v. 1.02) was carried out for each pan using the default settings, except for adjusting the delay limit to 1 and data normalisation with the percentileZ function. To reduce the complexity, only OTUs with > 1% relative abundance in at least one sample and present with more than 10 reads in at least thr ee differ ent subsampled samples wer e included in the network analysis. For both global (Spearman's rank correlation coefficients [SSCC]) and local associations (local similarity scores [LS]), only str ongl y significant ( P < 0.01 and q < 0.01) correlations were included. Network visualisation was carried out with Cytoscape v. 3.8.2 using the edge-w eighted, spring-embedded lay out (Shannon et al. 2003 ).
To identify OTUs in important network positions and to quantify their centrality, we used the weighted topological importance (WI) measur e gener alised b y Jor dán et al. ( 2006 ). This index calculates the number of neighbours and the number of their neighbours, in an ad diti v e and m ultiplicativ e way, while considering the strength of the interactions. To differentiate k e ystone OTUs (i.e . those O TUs that play a k e y role in the network and their removal would drastically impact the structure of the netw ork), w e consider ed indir ect inter actions up to thr ee steps (WI 3 ) and selected OTUs with WI i 3 > 1 (Jordán et al. 2006 , Berry andWidder 2014 ). If less than six OTUs fulfilled this r equir ement in a netw ork, the selection w as expanded to include OTUs with WI i

En vironmental par ameters
Meteor ological data r e v ealed typical seasonal air temper atur e dynamics with an incr easing tr end in spring (rate: + 0.16 • C/day, mean: 14.1 • C), no trend in summer (mean: 22.3 • C) and a decreasing trend in autumn (rate: -0.17 • C/day, mean: 11.3 • C) (Fig. 2 A). T he measured en vironmental parameters had values and follo w ed tr ends pr e viousl y described for the soda pans of this region ( Table S1 ) (Boros et al. 2014, Felföldi 2020, Szabó et al. 2020. In general, for samples collected in spring, similar environmental parameters have been measured, while in summer and autumn v ariation incr eased both between sampling times and sites ( Fig. 2 B , Figure S1 ). W ater depth v aried gr eatl y (1.5-46.0 cm) during the sampling period with the deepest water le v els measur ed during spring. Mor eov er, Zab-szék and K elemen-szék wer e completely dry on some occasions (sampling time 7, 8, 10, 11 and 13 for Kelemen-szék; and 11 and 13 for Zab-szék), making water sampling impossible. According to the PCA biplot, water depth was negativ el y r elated to the concentr ation of soluble compounds, pH, as well as copepoda and cladocera abundance, which had the highest values in mid-summer, while pH, DO and TP were elevated in summer and early autumn ( Fig. 2 B, Figure S1 ). Despite the common tr ends, onl y salinity and TN, and salinity and DOC, wer e highl y collinear. Although the PCA did not markedly distinguish the environmental parameters of the brown from the turbid pans ( Fig. 2 B), Sós-ér had, on av er a ge, significantl y deeper waters and higher TN concentrations as well as the highest median DOC, and lo w est TP and pH ( Figure S1 ). For a detailed description of the trends in environmental parameters, check Text S3 .
The seasonal community dynamics of the brown Sós-ér differ ed fr om those of the turbid pans . Here , in spring, a bOTU belonging to the famil y Erysipelotric haceae (Er) was dominant, while in the second half of summer there was a c y anobacterial bloom by a filamentous nitrogen-fixing Nodularia _PCC-9350 (No) bOTU. Some eO TUs , such as the members of Chrysophyceae Clades D (Cd) and F (Cf) or the parasitic fungus genus Pythium (Py), were also only dominant ( > 1%) in Sós-ér (Fig. 3 ). Meanwhile, the turbid pans had r elativ el y similar comm unity dynamics; the onl y exception were the drastic shifts in the microeukaryotic community composition observed following the desiccation-refillment events in the drying pans (i.e. Kelemen-szék and Zab-szék) when specific eOTUs became dominant (e.g. after the first desiccation in Kelemen

Core microbial community
From the almost 9000 identified O TUs , only 97 eO TUs and 191 bO-TUs were detected in all five pans; ho w ever, these core5 OTUs r epr esented 62% and 67% of the 18S and 16S rRNA gene reads, r espectiv el y ( Figur e S2 ). The cor e5 OTUs wer e pr edominant in turbid pans but not in Sós-ér, where core5 eOTUs and bOTUs constituted only 30% and 51% relative abundance, respectively (Fig. 4 A). When considering all the time points, the contribution of core5 eOTUs to the communities of Sós-ér was substantially lo w er than to the turbid pans ( P < 0.001), while, despite the lo w er contribution of core5 bOTUs to Sós-ér during spring, there was no significant difference between the brown and the turbid pans regarding the r elativ e abundance of core5 bOTUs ( P > 0.05).
The core4 OTUs shared only by the turbid pans consisted of 5-10 times more OTUs than core5 (952 eOTUs and 988 bOTUs) and r epr esented 80% and 84% of the r espectiv e r eads (Fig. 4 B). Differences among the turbid pans in respect of core OTUs were also detected. Mor e pr ecisel y, the r elativ e abundance of non-cor e OTUs was higher in the drying pans (Kelemen-szék and Zab-szék) than in the non-drying pans (Böddi-szék and Pan no. 60), and this difference was more notable for eOTUs than for bOTUs ( Fig. 4 A and  B). Furthermore, the contribution of both core5 and core4 eOTU r eads went markedl y down after dr ought e v ents, while cor e5 and core4 bOTUs sho w ed a r elativ el y stable contribution ov er time, irr espectiv e of desiccation ( Figure S3 ).

Dri v ers of community changes
T he 'en vfit' anal ysis significantl y ( P < 0.05) fitted salinity, pH, DOC, TN, TP and DO on the NMDS plots of both the microeukaryotic and bacterial communities. Ho w ever, w ater depth and Daphnia magna abundance wer e significantl y fitted onl y on the micr oeukaryotic NMDS, while water temper atur e and c hlor ophyll wer e onl y significant for bacterial communities (Fig. 4 ). The importance of DOC, TN and TP for both microeukaryotic and bacterial communities, and water temper atur e for bO TUs , w as enfor ced b y significant Mantel tests, although water temper atur e and DOC for bOTUs were onl y mar ginall y significant ( P = 0.046). Meanwhile, salinity and DO wer e onl y significant for eOTUs ( Table S3 ). The season of sampling had a significant effect on the communities (PERMANOVA Microeukaryotes: R 2 = 0.144, P = 0.001; Bacteria: R 2 = 0.115, P = 0.001) and the differences between seasons were also significant, with the strongest differentiation of spring communities from summer and autumn ( Table S4 ). Interestingly, the last autumn samples (sampling 14) appear close to the spring samples on both microeukaryotic and bacterial NMDS plots, suggesting similarity between the spring and late autumn samples (Fig. 5 ). This was also supported by the time distance decay curv e, whic h decr eased a gain for the sample pairs with a high time difference in almost all cases, except for the microeukaryotic communities of the two drying pans (Fig. 6 A). The NMDS plots and the PER-MANOVA analysis testing the impact of soda pan subtype (i.e. br own or turbid), clearl y separ ated the communities of the brown Sós-ér from those of the turbid pans (Microeukaryotes: R 2 = 0.086, P = 0.001; Bacteria: R 2 = 0.145, P = 0.001).
For cor e comm unities, the tw o-w ay PERMANOVAs assessing the effect of pan identity and sampling season, irr espectiv e of the analysed domain (microeukaryotes or bacteria) and the number of lakes included (all five or only four turbid), always explained  oeukaryotic and (B) bacterial OTUs with > 1% r elativ e abundance in at least one sample of the given pan coloured according to the corresponding phyla. OTUs with > 5% relative abundance in at least one sample are highlighted and indicated by the abbr e viated name of the corr esponding micr oeukaryotic gener a or bacterial clade, r espectiv el y. A k e y to the abbr e viations can be found in Supplementary Table S2 . mor e v ariance than for the non-cor e comm unities (Fig. 4 C and D ,  Table S5 ). Furthermore, seasonality explained more variance for cor e comm unities than for non-cor e comm unities, particularl y in the case of microeukaryotes (Fig. 4 C

Comm unity turnov er
Both seasonality and lake identity had a significant effect on turno ver (i.e . BC dissimilarity between sampling times) ( Table S6 ). Spring was c har acterised by significantl y lo w er values than summer and autumn, indicating a period of stability in spring (Fig. 6 B). Mor e pr ecisel y, while in Sós-ér the micr oeukaryotic turnov er was high throughout the study period, in the turbid pans in spring the microeukaryotic BC dissimilarity between sampling times was mostly low ( < 0.5). In late spring and early summer, eO TU turno ver substantiall y incr eased for P an no. 60, Zab-szék and K elemen-szék and gr aduall y for Böddi-szék, indicating major shifts in the micr oeukaryotic comm unity structur e of the turbid pans between the two seasons. While in June and July there were single consecutive sampling times with low BC dissimilarity in Zab-szék and Pan no. 60, starting from mid-July to the end of the study period, micr oeukaryotic turnov er r emained high ( > 0.5) in all pans. In the case of the drying pans (i.e. Zab-szék and Kelemen-szék), the eOTU BC dissimilarity was e v en higher, indicating that the micr oeukaryotic comm unity composition after the refillment was dr asticall y differ ent fr om the comm unities befor e the dr ought, which was further corroborated by the significant difference  ( P = 0.001) in the microeukaryotic turnover in this period between the non-drying and drying soda pans.
The bacterial BC dissimilarity between consecutive samplings was ov er all lo w er (mean 0.5) than for eOTUs (mean 0.6) (Fig. 6 B). Bacterial turnov er incr ease between spring and summer occurred only for the brown Sós-ér, and Kelemen-szék, the turbid pan with the most desiccation e v ents . After this , Sós-ér maintained high bacterial turnover until mid-August, while the turnover of the turbid pans r emained r elativ el y low throughout the summer. In the drought period of late summer and autumn, opposite to the micr oeukaryotic comm unities , bacterial turno ver of the drying pans was not higher ( P = 0.162) than those of the non-drying pans . T he onl y desiccation-r efillment e v ent with notable bacterial turnov er increase was the first drought of Kelemen-szék.

Microbial interactions
All networks had mor e positiv e corr elations than negativ e, irr espective of the pan or type of network ( Table S7 ). The networks of the non-drying turbid pans (i.e. Böddi-szék and Pan no. 60) had various properties not shared with the networks of the other pans: (1) both their sync hr onous and time-shifted networks had two distinct clusters of hubs densely connected with mostly SSCC Figur e 6. (A) Bra y-Curtis (BC) dissimilarity index between consecutive samplings as proxy for community turno ver. T he x-axis indicates the date of the latter sampling date within consecutive sampling pairs. (B) Time distance decay curve of community dissimilarity during the sampling period. Bo xplots re present the pairwise BC dissimilarity index between samples of the same pan taken at the given sampling interval. The x-axis shows the sampling interval between sampling events. edges (i.e. global correlations) (Fig. 7 ), and (2) their synchronous and time-shifted networks had a similar number of nodes, but the time-shifted had more edges and higher density ( Table S7 ).
Meanwhile, the networks of the two drying turbid pans (i.e. Kelemen-szék and Zab-szék) also shared various similarities, partly in contrast to the networks of the other pans: (1) their sync hr onous networks had more edges and nodes, and were denser than their time-shifted networks; (2) the nodes of their sync hr onous netw orks w er e densel y connected with mostl y LS edges (i.e. local correlations); and (3) the topology of their time-shifted networks, especially in the case of Zab-szék, was more fragmented than their sync hr onous networks (Fig. 7 , Table S7 ).
The networks of the brown Sós-ér shared properties with the non-drying turbid pans such as the similarity of the sync hr onous and time-shifted network, and higher density of the time-shifted network. They also had similarities to the networks of the drying turbid pans, such as the higher number of nodes and edges in the sync hr onous netw orks. Ho w e v er, the networks of Sós-ér were also distinguished from all turbid networks because they had the lo w est numbers of edges, nodes and neighbours , and o v er all low density ( Table S7 ). The topology of both of the networks of Sós-ér was highly fragmented, including an additional hub corresponding to the Nodularia bloom period (Fig. 7 ). The pr eferr ed season of eac h OTU was defined as the season when the differ ence between their mean r elativ e abundance in the giv en season and their mean r elativ e abundance in the entir e sampling period was larger than the standard deviation of their r elativ e abundance in the entire sampling period ((MEAN season -MEAN study period ) > SD study period ). Green indicates spring as preferred season, y ello w summer and red autumn. Grey denoted OTUs that lacked a distinct preferred season. The k e ystone OTUs are distinguished by brighter colours. The colouring of the edges is based on the correlation type (grey edges = SSCC, pink edges = LS), while line type of the correlations corresponds to their direction (solid edges = positive correlation, dashed edges = negative correlations).
Although microeukaryotic and bacterial OTUs were both identified as k e ystone OTUs in all networks, the majority of k e ystone OTUs were bacteria and there were more bacterial k e ystone OTUs in the turbid pans than in Sós-ér (Fig. 8 ). Only onethird (36%) of eukaryotic k e ystones was core5, while for the turbid pans 58% were core4. Among the k e ystone bO TUs , the majority wer e cor e5 (79%) and almost all wer e cor e4 (91%) in the turbid pans. Many k e ystone OTUs were assigned to phyla with high r elativ e abundance suc h as Chlor ophyta, Oc hr ophyta and Fungi for k e ystone eOTUs and Actinobacteria for bO TUs . Ho w e v er, Cy anobacteria w er e underr epr esented among k e ystones with e v en v ery abundant taxa suc h as Nodularia _PCC-9350 and Synechococcus _MBIC10613 was not identified as a k e ystone. Furthermor e, the highl y abundant Erysipelotric haceae OTU in Sósér was also not a k e ystone . Meanwhile , among the less abundant ( < 0.1% of reads) k e ystones, there were two eOTUs assigned to the parasitic Cryptomycotina order.
Ther e wer e v arious differ ences among positiv e and negativ e k e ystone OTUs (Fig. 8 ). First, positive keystones were more common among those with highest abundance in spring (74% of spring k e ystones), in summer the n umber of positi ves and negatives was similar, while among those abundant in autumn, negativ es wer e mor e common (80%). For bO TUs , taxonomic differences between positives and negatives were also noted. On the phylum le v el, Bacter oidetes (10) was the most common among positi ve k e ystones, follo w ed b y Pr oteobacteria (se v en) and then Actinobacteria (six), while for the negative OTUs Actinobacteria was the most common (17), follo w ed b y Proteobacteria (13) and then Bacter oidetes (six). Differ ences wer e ob vious also at lo w er taxonomic le v els as negativ e actinobacterial k e ystones wer e pr esent in high abundance throughout the study period and primarily belonged to the Nitriliruptoraceae (10) and acIII-A1 (three) lineages, while the positives were only abundant in spring and belonged to Luna1-A (three) and acIV-C (one). Similarly, the most common famil y for negativ e pr oteobacterial k e ystones was summer and autumn abundant Rhodobacteraceae, while for positives it was Burkholderiaceae (five) in all three seasons.
No substantial taxonomic difference was notable between k e ystone eOTUs of the sync hr onous and time-shifted netw orks. Ho we v er, for bacterial k e ystones, Actinobacteria were more common in the sync hr onous networks than in the time-shifted ones (16 vs. 8), while Bacteroidetes k e ystones wer e mor e common in timeshifted networks (10 vs. 6).

Discussion
Planktonic micr obial comm unities of fiv e shallow soda lakes exposed to identical climatic and meteorological conditions were e v aluated by a sync hr onous time series analysis . T he results demonstrated that the pans shared a core microbial community and similar seasonal dynamics, both in respect of community composition and interactions. Ho w ever, the extent of shared micr obiome and tr ends was not uniform acr oss the pans. Substantial differ ences wer e identified based on habitat subtype (i.e. br own or turbid). Common seasonal succession tr ajectories wer e also disrupted by local stressors (e.g. desiccation). Such events prompted a stronger response from microeukaryotic than from bacterial communities.

Differences of the brown soda pan
Although both the brown and turbid pans investigated in this stud y re present characteristic soda pan habitats (i.e. high pH, spe-cific ion composition and shallowness) (Boros et al. 2014 ), the brown Sós-ér is also c har acterised by extr emel y high colour ed DOC matter and DOC concentrations, and abundant shoreline v egetation (Bor os et al. 2020 ). Furthermor e, it has been pr e viousl y demonstrated that Sós-ér harbours markedly different microbial communities than the turbid Zab-szék (Szabó et al. 2017(Szabó et al. , 2020. In our study, Sós-ér also had significantly different microeukaryotic and bacterial communities as well as significantly deeper waters, and higher TN concentrations combined with low TP and pH and high DOC le v els. Our anal yses (i.e . en vfit and Mantel tests) suggested that the primary drivers of the distinctive microbial communities of the brown pan in this period were the different nutrient concentrations (i.e. TP, TN, DOC). Sós-ér also had a substantially lo w er number of shar ed OTUs (cor e5) than the turbid pans. The contribution of core5 eOTUs to the Sós-ér community was low throughout the study period and various unique microeukaryotic taxa were detected with high r elativ e abundance only in this pan, indicating that different soda pan subtypes exert strong selective forces on microeukaryotes . Meanwhile , from mid-summer, the contribution of core5 bOTUs to the Sós-ér community was similarly high as in the turbid pans, suggesting for bacteria more similar selection processes and no dispersal limitation between the five pans. Sós-ér also had v arious strikingl y differ ent intrinsic c har acteristics suc h as fr a gmented, low density and connectivity netw orks, as w ell as a high turno ver of eO TUs that indicate highly dynamic communities. Interestingly, the bOTUs of Sós-ér had a high turnov er onl y during summer and autumn, while in spring the bacterial turnover was low, together with the low core5 bOTU contribution in this period, suggesting a distinctive but stable spring bacterial community. This spring community was characterised by the high abundance of Erysipelotrichaceae UCG-004 bOTU, a taxon that was pr e viousl y r eported fr om soda pans of this region (Szabó et al. 2020 ), but has been otherwise described primaril y fr om the digestiv e systems of mammals and insects (Tegtmeier et al. 2016, Cox et al. 2017, Wu et al. 2021 ).

Common seasonal trends
Seasonality significantly affected the microeukaryotic and bacterial communities of all five pans and all studied seasons had distinctiv e comm unities. Furthermor e, despite the differences of Sósér, various common attributes of seasonal microbial succession could be identified. First, a circular trajectory of seasonality indicated by the resemblance of early spring and late autumn communities was detected for all pans, which was driven b y lo w w ater temper atur e and deeper waters and consequently low pH and low concentration of solutes . T he seasonal variation explained by the cor e vs. non-cor e comm unities also sho w ed similarities, irrespective of including Sós-ér in the analyses (i.e. core5 or core4). Mor e specificall y, the higher seasonal v ariance explained by the core than by the respective non-core communities supported our hypothesis in respect of seasonal adaptation happening primarily through species recruitment from the core community.
Further similarities of the seasonality of the five pans were the analogous attributes of spring. Spring was the season that differentiated the most strongly from the others (i.e. summer and autumn), according to both community structure and interactions. In all five pans, spring was c har acterised by positiv e keystones . P ositiv e sync hr onous associations can reflect mutualistic and facilitative interactions, but also parasitism, predation or similar nic he pr efer ence. Ho w e v er, the dir ection of species corr elations is not always obvious, for example, similar niche preference can manifest both as positive synchronous correlation due to . The pr eferr ed season was determined as the season when the mean z-scor e transformed abundance of the OTU was positive. If more than one season had positive mean z-score transformed abundance, no preferred season was specified (unspecified). coexistence or as negative due to competitive exclusion. Similarly, predation can display as positive synchronous associations when a pr ey attr acts its pr edators to a habitat patch or as negative (both sync hr onous and time-shifted) when pr edators eliminate their pr eys (Barber án et al. 2011. Various positive k e ystone eOTUs abundant in spring were assigned to predatory fla gellates suc h as Colpodellida (Mylnik ov 2009 ) or intr acellular par asitic taxa suc h as Cryptomycotina (Letc her et al. 2017 ) and Perkinsozoa (Mangot et al. 2011 ). T hese , together with the significant effect of Daphnia magna abundance on microeukaryotic comm unity structur e, suggest an important r ole of top-down contr ols in spring. Meanwhile, the spring-abundant positive actinobacterial k e ystones belonged to lineages (i.e. Luna1-A and acIV-C) that ar e c har acterised by v ery small cell sizes ( < 0.1 μm 3 ) (Duda et al. 2012 ) and have been previously suggested to be grazing resistant (Tar ao et al. 2009, Ec kert et al. 2013, whic h might explain their coexistence with predatory eukaryotic taxa. Overall, the seasonal dynamics of actinobacterial keystones were similar to the seasonal dynamics described for this phylum in various other limnic systems (Eiler et al. 2012, Mikhailov et al. 2022. Except for the eOTUs of Sós-ér, spring was also c har acterised b y lo w species turnov er, whic h could be the r esult of r elativ e stability provided by deeper waters . T he spring samples had a high abundance of a single Choricystis eOTU, a widely distributed freshw ater picoeukary otic green algal genus (Kulak ov a et al. 2020 ). These r esults a gr ee with pr e vious studies of soda pans of the region that sho w ed that picoeukary otic gr een algae ar e the most abundant members of the phytoplankton and have a specific seasonal trend with the highest abundances in winter-spring and the lo w est in summer (Somogyi et al. 2011, 2016, Felföldi 2020, Szabó et al. 2020. The dominance of positive keystone taxa and the transitivity of positiv e corr elations combined with low turnover of both eOTUs and bOTUs in the turbid pans suggests that spring was a period when community assembly was primarily ruled by trophic interactions between primary producer eukary otic picoalgae (e.g. Choric ystis ), their parasites (e.g. Cryptomycotina), heter otr ophic bacteria (e.g. Luna1 lineages, Burkholderiaceae, Balneolaceae) consuming algal exudates and debris, and flagellates (e.g. Colpodellida) predating on bacteria.

The impact of local stressors
After the r elativ el y stable and sync hr onous period in spring, the variation between microbial communities increased even for the turbid pans . T his w as driven b y the higher and more variable concentr ations of dissolv ed substances r esulting fr om the lo w er w ater le v els, whic h a gr ees with pr e vious studies showing that environmental fluctuations induced by shrinking ecosystem size modulate comm unity assembl y pr ocesses and r educe stability (Bier et al. 2022 ). In the case of Kelemen-szék and Zab-szék, the decreasing water le v els r esulted in v arious desiccation e v ents follo w ed by r efillments. Drying-r e wetting cycles exert se v er e str ess on micr oor ganisms due to drastic changes in salt and nutrient content (Székely and Langenheder 2017, Schimel 2018, Truchy et al. 2020 ). While desiccation is common in soda pans of this region, not every pan dries out e v ery year and it is not always the same pans that dry out (Boros et al. 2020, Szabó et al. 2020 ), making desiccation not a part of the regular seasonality, but rather a local stressor.

Comparison of microeukaryotic and bacterial trends and stress response
Intense environmental fluctuation combined with increased gr owth r ates due to summer warming and decr easing habitat size (i.e. shrinking water le v els) wer e expected to stimulate micr obial turnov er (Vass et al. 2021, Bier et al. 2022. Accordingly, the turnover of microeukaryotic communities in the turbid pans incr eased r elativ el y uniforml y. Micr oeukaryotic turnov er also incr eased substantiall y as a consequence of eac h drying-r e wetting cycle, suggesting limited resilience of the microeukaryotic communities to such stressors. Simultaneously, the contribution of non-core eOTUs to the microeukaryotic communities of the drying pans was m uc h higher than to those of the non-drying turbid pans and their r elativ e abundance incr eased, particularl y following drought events, supporting our hypothesis about non-core OTUs becoming more abundant in response to local stressors. A potential explanation for this phenomenon is that, despite the extensiv e soda pan-ada pted micr oeukaryotic cor e comm unity, dr astic str ess e v ents like desiccation disrupt the species-sorting processes from the core community due to dispersal limitation or a lack of internal drought-resistant seed banks. Meanwhile , the turno ver of the bacterial communities of the turbid pans r emained r elativ el y similar through the study irrespective of desiccation, implying that bacteria were more resistant to the desiccation stress and the overall more extreme conditions of summer-autumn than microeukaryotes. Bacterial k e ystones of this period also suggest special adaptations to extreme conditions . For example , it has been shown for close r elativ es of the summer-autumn abundant Nitriliruptoraceae k e ystones that they can scav enge or ganic nitr ogen e v en fr om str ong nitrile bonds (Sorokin et al. 2009 ), allowing them versatility to overcome nitrogen limitation. The high number of negative bacterial keystones in summer-autumn also suggests that bacterial groups with different optima were dynamically outcompeting each other under the quic kl y c hanging conditions . T he contribution of core4 bOTUs was also very high at every sampling time and site, suggesting that the core bacterial community of turbid pans not only contributed to the adaptation to seasonal changes, but was also highly resistant to extreme conditions. As it has been shown that dryingr e w etting c ycles hav e str ong filtering effects on bacterial comm unities and dispersal is r equir ed for full r ecov ery (Fazi et al. 2008, Székely and Langenheder 2017, the similar contribution of core4 bOTUs to the drying and non-drying turbid pans, as well as the uninterrupted high contribution of core4 bOTUs irrespectiv e of dr ought e v ents, suggests no dispersal limitation for bacteria between the studied pans. Recently, it has been suggested that w aterbir ds play an important role in dispersing both prokaryotes and microeukary otes betw een soda pans (i.e. endozoochory) (Szabó et al. 2022 ), which implies that dispersal intensity would vary with bird visitation frequency. Although our study cov er ed periods with different bird abundances (Boros et al. 2023 ), the absence of indications of bacterial dispersal limitation suggests that the bacterial cor e micr obiome is highl y dispersed between the pans by non-endozoochory dispersal such as wind or precipitation (Langenheder and Székely 2011 ).
T he o v er all higher contribution of cor e bOTUs to k e ystone taxa than that of core eOTUs indicates more synchronised community trends for bacteria than for microeukaryotes. In general, micr oeukaryotic comm unities wer e mor e sensitiv e to the local str essors, pr obabl y both due to less physiological resistance and dispersal limitation. All in all, this supports our hypothesis regarding the lesser impact of local str ess e v ents on bacterial compared with micr oeukaryotic comm unities. Our r esults ar e also in a gr eement with studies suggesting that dispersal limitation and stochasticity ar e mor e important in sha ping micr oeukaryotic comm unities, while selection pr ocesses ar e mor e pr ominent in bacterial comm unity assembl y and with studies that identified incr eased dispersal limitation during dry periods for microeukaryotic communities (Beisner et al. 2006, Wang et al. 2015, Logares et al. 2018, Chen et al. 2019, Mikhailov et al. 2022 ).

Effect of seasonal trends and local stressors on microbial interactions
The networks generated in this study allo w ed for joint analyses of the interactions of microeukaryotes and bacteria. All networks had mor e positiv e corr elations than negativ e, suggesting pr edominance of positive interactions in the comm unities. Pr e vious studies demonstrated that positive associations are more common in ecosystems c har acterised by high abiotic stress due to a higher number of mutualistic interactions permitting species to exist in harsher environments than would otherwise be possible (Travis et al. 2005, Hernandez et al. 2021. Ho w e v er, the dependence on mutualism makes such communities sensitive to perturbations and r educes network stability, particularl y in the case of low modularity networks (Hernandez et al. 2021 ). Apart from the dominance of positive correlations, there were clear differences in network topology and properties between the brown Sós-ér, the two drying (Kelemen-szék and Zab-szék) and the two non-drying turbid pans (Böddi-szék and Pan no.60). The networks of the latter pans reflected the most consistent and stable seasonal succession processes, with distinct hubs corresponding to different seasons (i.e. spring and summer-autumn) and a high number of mainly negative global (i.e. SSCC) time-shifted associations. For the drying pans, seasonal network modularity was less clear and dominance of sync hr onous local (i.e. LS) associations was the most c har acteristic, indicating lesser importance of community dynamics overarc hing the entir e study period. Meanwhile, the networks of Sósér had three hubs, with one corresponding to the late summer period c har acterised by Nodularia bloom. Ov er all, the inter action networks of the soda pans reflected low community stability in these high stress habitats that was further exacerbated by local str ess e v ents, suc h as drying-r e w etting c ycles or c y anobacterial blooms.

Conclusions
By integrating network analyses, assessment of k e ystone taxa and the consideration of microeukaryotic and bacterial communities separ atel y, we gained novel and comprehensive insights into the unique seasonal dynamics of shallow soda lakes. Our findings revealed that despite significant environmental changes and subsequent community shifts, the studied soda pans were primarily inhabited by a common core microbiome and share certain c har acteristics of their seasonal trends. Ho w ever, the extent of the shar ed micr obiome was curtailed among pans of different habitat subtype (i.e. brown or turbid pan) and local stress events like desiccation and refillment, modified common seasonal patterns. In gener al, stable envir onmental conditions during spring fostered stable micr obial comm unities gov erned by tr ophic inter actions. Conv ersel y, the dynamic c hanges of summer and earl y autumn, coupled with local str ess e v ents, exerted str ong selectiv e pr essure that instigated varied response mechanisms in microeukaryotic and bacterial communities. Species recruitment from the cor e comm unity played a k e y r ole in the ada ptation to seasonal changes for both microeukaryotes and bacteria, while for micr oeukaryotes r esponse to str ess e v ents involv ed a lar ge part of the non-core communities, suggesting the influence of dispersal limitation in their r ecov ery. By contr ast, bacterial comm unities were, to a great extent, resistant to stressors and adapted to extr eme conditions thr ough species sorting fr om the cor e comm unity and competitive exclusion. To the best of our knowledge, this provides the first evidence in extreme aquatic habitats in support of the hypothesis that microeukaryotic communities exhibit higher sensitivity to str ess e v ents compar ed with bacteria. Furthermor e, our study underscor es the importance of separ atel y considering microeukaryotic and bacterial communities when assessing the impacts of local stressors . T his is particularly relevant as e v ents like desiccation and r efillment become mor e fr equent in shallow aquatic ecosystems due to climate change.