Miriam Beck
email: miriam.beck@imev-mer.fr

Caroline Cailleton

Lionel Guidi

Corinne Desnos

Laetitia Jalabert

Amanda Elineau

Lars Stemmann

Sakina-Dorothée Ayata

Jean-Olivier Irisson

Morphological diversity increases with decreasing resources along a zooplankton time series

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

seldom explicitly considered, except for body size. We describe morphological diversity of marine zooplankton seasonally and over 12 years using quantitative imaging of weekly plankton samples, in the North-Western Mediterranean Sea. We extract 45 morphological features on greater than 800 000 individuals, which we summarize into four main morphological traits (size, transparency, circularity and shape complexity). In this morphological space, we define objective morphological groups and, from those, compute morphological diversity indices (richness, evenness and divergence) using metrics originally defined for functional diversity. On both time scales, morphological diversity increased when nutritive resources and plankton concentrations were low, thus matching the theoretical reciprocal relationship. Over the long term at least, this diversity increase was not fully attributable to taxonomic diversity changes. The decline in the most common plankton forms and the increase in morphological variance and in extreme morphologies suggest a mechanism akin to specialization under low production, with likely consequences for trophic structure and carbon flux.

KEYWORDS

diversityzooplankton community morphologyoligotrophy oligotrophymorphology quantitative imaging traits zooplankton communitydiversity

Numerous studies indicate that in aquatic and terrestrial ecosystems biodiversity is positively linked to productivity, for example, reporting increasing primary production due to increased resource use efficiency [START_REF] Setubal | Effects of functional diversity and salinization on zooplankton productivity: an experimental approach[END_REF][START_REF] Tilman | Diversity and productivity in a long-term grassland experiment[END_REF]. Biodiversity is in turn influenced by the overall production of a system, as shown by the linear decline in taxonomic richness of phyto-and zooplankton when nutrient concentrations increase [START_REF] Declerck | Plankton biodiversity along a gradient of productivity and its mediation by macrophytes[END_REF]. Still, in arctic soft bottom communities, the relationship seem to be more complex: while invertebrate diversity decreased at high production levels, it was also low at very low production levels, resulting in a unimodal shape on a local scale [START_REF] Witman | The relation between productivity and species diversity in temperate-arctic marine ecosystems[END_REF]. Also, a plankton community study [START_REF] Korhonen | Productivity-diversity relationships in lake plankton communities[END_REF] and a meta-analysis [START_REF] Mittelbach | What is the observed relationship between species richness and productivity?[END_REF] revealed varying slopes in this relationship, depending on the spatial scale and organism type, underlining the contextdependence and complex mechanisms acting on it.

Depending on the question addressed, biodiversity is quantified by various indices describing, for example, taxonomic, phylogenetic or functional diversity. The latter, relying on functional traits instead of taxonomic identity, captures the range of ecological strategies present at a site and is thus directly linked to ecosystem functioning [START_REF] Schmera | Characterising functional strategies and trait space of freshwater macroinvertebrates[END_REF][START_REF] Vallina | Phytoplankton functional diversity increases ecosystem productivity and stability[END_REF]. This allows insights often not captured when relying solely on taxonomic indices. For instance, in winter, a phytoplankton community can remain functionally diverse despite declining taxonomic richness [START_REF] Weithoff | Comparing seasonal dynamics of functional and taxonomic diversity reveals the driving forces underlying phytoplankton community structure[END_REF] and plant assemblages can experience a substantial loss in functional diversity over the long term although species richness increases [START_REF] Nielsen | More is less: net gain in species richness, but biotic homogenization over 140 years[END_REF]. The taxonomic and functional approaches thus complement each other and provide a better understanding of community structuration by identifying ecological processes such as competitive exclusion [START_REF] Leruste | Complementarity of the multidimensional functional and the taxonomic approaches to study phytoplankton communities in three Mediterranean coastal lagoons of different trophic status[END_REF] or functional homogenization [START_REF] Magurran | Rapid biotic homogenization of marine fish assemblages[END_REF].

Among the commonly used functional traits, morphological traits, describing part of the phenotype of an organism, are scarcely considered, except for body size. Indeed, the latter is considered a 'master trait', since most physiological rates, e.g. metabolism [START_REF] Brown | Toward a metabolic theory of ecology[END_REF][START_REF] Kleiber | Body size and metabolism[END_REF], growth or feeding [START_REF] Hansen | Zooplankton grazing and growth: scaling within the 2-2,000μm body size range[END_REF] scale with size. Size itself is affected by environmental conditions and, for example, warming and nutrient enrichment locally increased the proportion of small phyto-and zooplankton, macroinvertebrates and fish, at both species and community levels [START_REF] Latli | Long-term trends in trait structure of riverine communities facing predation risk increase and trophic resource decline[END_REF][START_REF] Rasconi | Increasing water temperature triggers dominance of small freshwater plankton[END_REF][START_REF] Strecker | Effects of experimental greenhouse warming on phytoplankton and zooplankton communities in fishless alpine ponds[END_REF]. Such shifts in the size-structure of a community, with consequences on its functioning, can occur with no change in taxonomic diversity. This was documented for benthic macroinvertebrates for which a lower proportion of large organisms decreased decomposition rates without affecting Shannon diversity and evenness [START_REF] Dossena | Warming alters community size structure and ecosystem functioning[END_REF]. The positive relationship between predator size and the intensity of trophic cascades is another example of how body size influences food web structure and energy flux [START_REF] Delong | The body size dependence of trophic cascades[END_REF]. Using morphological features beyond size, Marinello & Bernhard [START_REF] Marinello | Wing morphology of Neotropical bats: a quantitative and qualitative analysis with implications for habitat use[END_REF] showed that wing shape of bats can be a proxy for their flight pattern and feeding strategy, with carnivores having rather short, round wings of a large surface allowing them to lift and carry prey. Sonnet et al. [START_REF] Sonnet | Length, width, shape regularity, and chain structure: time series analysis of phytoplankton morphology from imagery[END_REF] found a negative correlation of phytoplankton shape complexity and size with nutrient concentration, which suggests that considering phytoplankton shape can improve our understanding of the seasonal succession.

Introduction

providing insight into the functioning of planktonic ecosystems during the ice break up in the Arctic. In mosquitofish, measurements of morphological traits including eye position, head size and elongation confirmed expected differences between individuals in high and low-predator environments, which enhance swimming performance under predator pressure [START_REF] Langerhans | Predator-driven phenotypic diversification in gambusia affinis[END_REF].

Ecosystem functioning worldwide is, and will continue to be, affected by climate change (IPBES and IPCC; [START_REF] Pörtner | IPBES-IPCC co-sponsored workshop report on biodiversity and climate change[END_REF]). Beyond temperature rise, indirect consequences of warming on marine environments include reduced nutrient and oxygen availabilities in surface layers [START_REF] Kwiatkowski | Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections[END_REF] and a global decline in primary production and animal biomass [START_REF] Bryndum-Buchholz | Twenty-first-century climate change impacts on marine animal biomass and ecosystem structure across ocean basins[END_REF][START_REF] Lotze | Ensemble projections of global ocean animal biomass with climate change[END_REF]. Many marine ecosystems will experience drastic species turnover [START_REF] Cheung | Projecting global marine biodiversity impacts under climate change scenarios[END_REF] with smaller organisms and/or taxa likely favoured in most cases [START_REF] Daufresne | Global warming benefits the small in aquatic ecosystems[END_REF][START_REF] Lefort | Spatial and body-size dependent response of marine pelagic communities to projected global climate change[END_REF][START_REF] Moullec | An end-to-end model reveals losers and winners in a warming Mediterranean Sea[END_REF]. Increasing nutrient concentrations often cause functional and/or taxonomic homogenization through the dominance of generalist taxa [START_REF] Chihoub | Coastal eutrophication as a potential driver of functional homogenization of copepod species assemblages in the Mediterranean Sea[END_REF][START_REF] Li | Eutrophication and predation mediate zooplankton diversity and network structure[END_REF][START_REF] Zhang | Nutrient enrichment homogenizes taxonomic and functional diversity of benthic macroinvertebrate assemblages in shallow lakes[END_REF], suggesting community shifts towards specialists and increasing diversity with oligotrophication, as already observed for a stream phytoplankton community over 34 years [START_REF] Abonyi | Functional diversity of phytoplankton highlights long-term gradual regime shift in the middle section of the Danube River due to global warming, human impacts and oligotrophication[END_REF].

In the Mediterranean Sea, a hotspot of climate change, many regions endure multiple stressors such as increasing temperature and salinity, extreme events and continuing oligotrophication [START_REF] Garrabou | Marine heatwaves drive recurrent mass mortalities in the Mediterranean Sea[END_REF][START_REF] Ramírez | Spatial congruence between multiple stressors in the Mediterranean Sea may reduce its resilience to climate impacts[END_REF]. Investigating ecosystem responses to climate change in this area hence both important and informative. In the North-Western basin, temperature increased while nutrients declined over the last two decades [START_REF] Feuilloley | Concomitant changes in the environment and small pelagic fish community of the Gulf of Lions[END_REF]. We thus expect zooplankton in this region to respond with increased morphological diversity due to species turnover and/or intraspecific variations as the community shifts towards individuals with more specialized morphologies under increasingly oligotrophic conditions. To test this hypothesis, we quantified zooplankton morphological diversity on seasonal and inter-annual time scales, along a 12-year times series collected at a site in the NW Mediterranean. We synthesized the morphological features from greater than 800 000 images of individual zooplanktonic organisms and defined groups of morphologically similar individuals. Based on these 'morphs' we characterized morphological diversity and investigated the impact of environmental drivers. . Samples were fixed with formalin buffered with borax before further treatment. At the same site, the following environmental variables: water temperature (°C), salinity (psu), dissolved oxygen (m1 1 ), nitrate ( ; µmol l ), nitrite ( ; µmol l ), phosphate ( ; µmol l ), chlorophyll a (μg l ), pH, particulate organic carbon and particulate organic nitrogen (POC, PON; μg l ), were measured once per week, in the morning, at six depths, as part of the SOMLIT network following standardized protocols (http://somlit.fr). To match biological and environmental samples in frequency and sampling depth, zooplankton samples were pooled per week while environmental variables were integrated between 75 m and the surface.

The workflow summarizing the steps from the imaged individuals to the quantification of diversity are displayed in figure 1.

Zooplankton samples were fractionated on a 1 mm filter to separate large and small objects. Each fraction was subsampled (at 1/2, 1/4, 1/8 etc.) using a Motoda splitting box [START_REF] Motoda | Devices of simple plankton apparatus[END_REF] until approximately 1000 objects could be digitized using the ZooScan, a scanner designed for zooplankton. Scanned images (1 pixel = 10.6 µm) were automatically processed to extract and measure all objects with the ZooProcess software [START_REF] Gorsky | Digital zooplankton image analysis using the ZooScan integrated system[END_REF]. In total, 45 morphological features were extracted (electronic

Material and methods

(a) Sampling

3 -1 NO - 3 -1 NO - 2 -1 PO 3- 4 -1 -1 -1
Workflow from individual images to the quantification of diversity. Green font indicates the object handled in the next step; red a variable; blue the method or software. (b) Quantitative imaging supplementary material, S1; figure 1). They were related to size (e.g. perimeter, area, size of smallest rectangle or ellipse surrounding the object), shape (e.g. circularity, elongation, perimeter complexity), and the distribution of grey levels (e.g. contrast, minimum grey, modal grey). Images showcasing processing artefacts such as undersegmentation (two touching objects considered as one) were removed. Individual objects and their associated measurements were uploaded to the EcoTaxa platform for taxonomic identification [START_REF] Picheral | EcoTaxa, a tool for the taxonomic classification of images[END_REF]. Identifications were reviewed by several experts for the entire time series, to ensure consistency. In total, 44 taxa were identified with taxonomic levels ranging from species (only few) to broader taxonomic groups (e.g. Amphipoda, Ostracoda; full list in electronic supplementary material, S2). Juvenile and larval stages were treated separately when appearing in high numbers and otherwise were grouped with adults of their taxa. Zooids of colonial taxa were treated as the identified taxon. Largely incomplete objects (i.e. body parts) and rare taxa (appearing with n < 100 out of a total of n = 845 812 individuals over the whole time series) were excluded from the analysis. Samples with high numbers of scanned objects (n > 2000 for the sum of the two scans) were randomly subsampled (n = 2000) to avoid differences in imaging effort possibly influencing the assessment of diversity. The concentration of each taxon in each sample (ind m ) was calculated based on the image count, the Motoda fraction and the filtered water volume.

The 45 morphological features were summarized through an observation-weighted Principal Component Analysis (PCA; figure 1) computed on the 845 812 imaged individuals. Each image was weighted by its 'individual' concentration to be representative of its importance in the water rather than on the scan, hence erasing differences in Motoda fractionating rate or volume sampled at sea among weeks. The PCA synthesized their correlation and allowed placing each individual in a reduced 'morphological space'. To ensure the applicability of PCA histograms of each morphological feature were checked prior to the analysis. When necessary, outliers were removed by censoring the lower and/or upper 0.1% of the feature values. For not-normal distributions, a log transformation was applied to increase symmetry. Electronic supplementary material, table S1 details what was applied to each feature. After censoring, individuals with more than five removed values (n = 1448) were considered too extreme and excluded from the analysis. Other missing values (affecting n = 20,986 individuals) were replaced by the mean of the corresponding feature, which is neutral in PCA. The number of PCA axes worthy of further examination was determined by judging the scree plot. Each of the selected PCA axes was then interpreted as a morphological trait, with a meaning defined according to the morphological features most associated with it. To follow the variation in morphology within the community over time, for each date, we calculated the mean, standard deviation, 0.25 and 0.75 quantile of individuals' coordinates along each axis.

Finally, groups of morphologically similar individuals (i.e. 'morphs') were defined, to be then used to quantify morphological diversity: Observation-weighted k-means clustering [START_REF] Hartigan | Algorithm AS 136: a K-Means Clustering Algorithm[END_REF] was applied on the individuals' coordinates along the selected PCA axes (figure 1); weights were again the 'individual' concentrations. The number of clusters was set to 200, which was identified as a good compromise between the resolution in the definition of morphs and their generality; defining 400 morphs only marginally changed the results. Morphs differed from taxa, since one morph comprised on average 15 taxa (ranging from 5 to 30) and each taxon appeared on average in 70 morphs (ranging from 5 for Pyrosomatida to 197 for Copepoda). For each date, the concentration of each morph was calculated as the sum of the concentrations of all individuals belonging to that morph.

Villéger et al. [START_REF] Villéger | New multidimensional functional diversity indices for a multifaceted framework in functional ecology[END_REF] developed three indices describing different and complementary aspects of functional diversity: richness, regularity and divergence. Using 'morphs' and their average morphological traits, instead of species and their average functional traits, we used those indices to quantify morphological diversity. Each index was calculated for each date, based on the position of the barycentre of each morph in the reduced morphological space (= the centroids of the k-means clusters) and its concentration (= sum of the 'individual' concentrations of cluster members) at that date. Morphological richness (MRic) represents the proportion of the total morphological diversity present at a given date. It is calculated as the volume of the convex hull enclosing all the morph's centroids present at that certain date. Thus, it does not take the morphs' concentrations into account and rather depends on the most extreme morphs. Morphological evenness (MEve), or regularity, is a measure of the homogeneity of the differences in appearances between morphs present at a given date. It is high if the appearance changes gradually between morphs; it is low if some morphs are very different from other morphs that are similar to each other. Numerically, the minimum spanning tree between the present morphs calculated and the equidistance along the branches of this tree defines the morphological regularity. Finally, morphological divergence (MDiv) quantifies morphological dissimilarity within the fraction of diversity present at a given date. MDiv is high when several morphs of very different appearance are abundant; it is low when the abundant morphs are morphologically similar. First, the gravity centre of morphs present at the target date was computed. Then, the distances from each morph to this gravity centre as well as their average were computed. The index is based on the sum of the anomalies of theses distances compared to their mean. Since everything is computed based only on the morphs present at a given date, the index quantifies morphological divergence within a certain level of diversity (i.e. richness).

Taxonomic diversity was quantified by calculating taxonomic richness (i.e. number of taxa present, TRic), Pielou Evenness (TEve ; [START_REF] Pielou | An introduction to mathematical ecology[END_REF]) and Shannon Diversity (TDiv; [START_REF] Shannon | A mathematical theory of communication[END_REF].

All time series were processed equally: morphological and taxonomic diversity indices, environmental variables, total zooplankton concentration and morphological change in the community (mean, standard deviation, 0.25 and 0.75 quantiles). Values missing over short periods (1-2 weeks; n = 2-3 per year) were filled by linear interpolation of the original variable. For gaps spanning a period ≥3 weeks (n = 1 in begin of 2017) missing values were interpolated by the analysis technique described thereafter. The seasonal component was extracted by Seasonal-Trend Decomposition by Loess (STL; [START_REF] Cleveland | STL: a seasonal-trend decomposition procedure based on loess[END_REF]). Since different de-seasonalizing methods within STL provided very similar results, we kept the simplest one and assumed seasonal component to be identical over all years ('periodic' signal).

After removing this seasonal signal, the significance of the long-term trend was tested using a generalized least squares (GLS) linear regression with a first-order autoregressive correlation structure in the residuals [START_REF] Fox | An R companion to applied regression[END_REF]. This level of autocorrelation was determined by examining the partial auto-correlogram of the residuals of a simple linear model. The assumptions of the linear model (normality, homoscedasticity) were then tested on the de-correlated residuals. When the assumptions were not met, the process was repeated on transformed data (e.g. log(n + 1) for zooplankton concentrations). While changes are not likely to be completely progressive and linear over the whole period (they may rather occur more intensely at some periods), this approach was used as a parsimonious way of statistically testing for overall changes.

The definition of a morphological space and morphs, quantification of diversity and time-series analyses were repeated considering only the copepods within the community (n = 649 825 images), to investigate whether the results persisted in a more taxonomically homogeneous community.

All analyses were performed using R 3.4.4 [START_REF]R: A language and environment for statistical computing[END_REF] and the following packages: ade4 1.7-10 [START_REF] Dray | The ade4 Package: implementing the duality diagram for ecologists[END_REF] for weighted PCA, vegan 2.4-6 [53] for calculating the Pielou and Shannon indices, pastecs 1.3-18 [START_REF] Grosjean | Pastecs: Package for analysis of space-time ecological series[END_REF] for interpolation and calculation of moving averages on time series, stlplus 0.5.1 [START_REF] Hafen | Package 'stlplus' Enhanced Seasonal Decomposition of Time Series by Loess. Version 05 1 Repos CRAN Nocopyright R Found Core Team[END_REF] for STL decomposition, nlme 3.1-131.1 [START_REF] Pinheiro | Package 'nlme. Linear nonlinear mixed effect models[END_REF] for GLS regression, dplyr 0.7.4 [START_REF] Wickham | dplyr: a grammar of data manipulation[END_REF] for data manipulation, and ggplot2 2.2.1 [START_REF] Wickham | ggplot2. R package version 2[END_REF] and factoextra 1.0.7 [START_REF] Kassambara | Factoextra: extract and visualize the results of multivariate data analyses[END_REF] for graphics. The representation of the morphological space (figure 2) was coded in https://github.com/jiho/morphr/ and the computation of observation-weighted k-means was coded in https://github.com/jiho/wkmeans. The morphological diversity indices were calculated after optimizing the functions of Sébastien Villéger (http://villeger.sebastien.free.fr/Rscripts.html). The entire code used for the analyses is available for reproducibility, at https://github.com/jiho/ptb_morphodiv.

(e) Times-series analysis (f) Software

Zooplankton morphological space defined by the first four PCA axes, computed from 45 variables measured on 845 812 individual images. For each plane, only the 12 most contributing features are drawn. Features are coloured according to the morphological facet they capture: red for shades of grey, blue for size, orange for symmetry and green for other aspects of the shape. The variance captured by each axis and its interpretation in terms of synthetic morphological trait are indicated in the axis label. Images representative of a region of the morphological space were synthesized by combining three images randomly selected in the neighbourhood and the result was drawn on the plot.

Figure 2.

The morphological space was built by the first four axes of the PCA on morphological features, which explain together 75.39% of the total variance of zooplankton morphology (37.86%, 21.29%, 9.4% and 6.84%, respectively) and represent the main morphological traits of the community. The first axis (PC1) mainly captured changes in body size ( figure 2). The second axis (PC2) can be interpreted as transparency separating light, i.e. transparent, individuals with homogeneous grey levels from darker and more contrasted individuals. The third axis (PC3) separated individuals by their circularity and symmetry, opposing elongated individuals with a marked bilateral symmetry and more circular ones. Finally, the fourth axis (PC4) quantified shape complexity and separated individuals with a simple perimeter from those with more complex contours.

Taxonomically, Limacinidae (Gastropoda) and Bivalvia dominated the lower end of the size continuum and Salpida and

Results (a) Description of the morphological space

Chaetognatha dominated the upper end (electronic supplementary material, S2). The transparency axis separated gelatinous taxa, such as Appendicularia (Fritillariidae, Oikopleuridae) and Chaetognatha from opaque, carapace-bearing (Eumalacostraca, Copepoda, Amphipoda) or shell-bearing (Gastropoda like Cavolinia inflexa) zooplankton. Chaetognatha and Creseis acicula (Gastropoda) expressed the highest elongation, while Salpida, Phaeodarea and pieces of Abylopsis tetragona (Siphonophorae) belonged to the most circular individuals. Taxa with simpler outlines were Salpida, Chaetognatha and Diphyidae (Siphonophorae), while the most complex shapes were observed for Echinodermata larvae, Phaeodarea and Copepoda with antennae in open position.

The seasonal component showed the typical succession of a plankton bloom: Chlorophyll a concentration peaked in March, followed by zooplankton concentration peaking in the beginning of April (figure 3). Temperature increased towards summer, while nutrient and particle concentrations declined, starting just before the chlorophyll peak for some (e.g. , ) or at the zooplankton peak for others (PON; figure 3; electronic supplementary material S3). The pattern for taxonomic and morphological diversity was similar for all indices (richness, evenness, divergence): diversity was low during spring (March, April) before sharply increasing until July and remaining stable at a high level over summer.

In fall, taxonomic diversity indices started to decrease while their morphological counterparts remained at a high level until the next spring (figure 3). Diversity, which was minimal in spring, was negatively correlated to zooplankton concentration, which was maximal at that time (electronic supplementary material, S4).

Seven of the eleven environmental variables showed significant trends over the 12-year period. Average temperature between 0 and 75 m increased significantly (+0.59°C in 12 years, assuming a linear increase, with p < 0.05), so did salinity (+0.11 PSU) and oxygen (+0.43 ml l ; table 1). Concentrations of phosphate (-0.02 mg l ), particulate organic carbon (-43.2 mg l ) and nitrogen (-2.28 mg l ) significantly declined (table 1). The concentrations of nitrite, nitrate and silicate, as well as density, did not show a significant monotonous trend over the 12 years. The concentration of chlorophyll a significantly increased (+0.16 mg l ) while the total zooplankton concentration significantly declined during this period.

zooplankton environment variable slope R2 sig. variable slope R2 sig.

(c) Long-term trends in zooplankton and environmental conditions

-1 -1 -1 -1 -1
Trends in the zooplankton community (concentration, morphological and taxonomic diversity indices, values of morphological traits) and environmental variables over the 12-year period (2009-2020) obtained from generalized least squares (GLS) regression. For each index and variable, the slope describes the direction and intensity of the trend per year, the R shows how much variance this linear model captures, and asterisks indicate significance with p < 0.001 (***), p < 0.01 (**) and p < 0.05 (*).

Table 1. In this context, taxonomic richness (TRic) significantly increased and the community became more morphologically diverse, both in terms of absolute (MRic) and relative (MDiv) diversity, i.e. the morphs present were more variable ( figure 4 and table 1). There was no significant trend in morphological evenness (MEve) ( p > 0.05, figure 4) and in taxonomic evenness (TEve) or diversity (TDiv).

The morphological traits of the community (i.e. the positions on the axes of the morphological space defined by the PCA) changed notably (table 1): The variance in size increased and the smaller individuals became smaller; the average and highly transparent organisms became more transparent and the overall variance in transparency increased; organisms became more elongated, both in average and at the extremes. Only shape complexity did not change significantly (table 1).

Analysing only Copepoda, the traits defining the morphological space were similar to those found for the entire zooplankton community: PC1 (36.46% of explained variance) and PC2 (19.77%) were size and transparency, respectively. The morphological interpretation of PC3 (8.94%) was less clear, but individuals with low scores on PC3 were more symmetric. PC4 (7.34%) was elongation (electronic supplementary material, S5).

As for the whole plankton community, copepod concentration peaked in April and then strongly declined (electronic supplementary material, S6). The three morphological diversity indices followed the opposite pattern: Their values were PO 3- lowest in April and increased over the year until July/August, fluctuating-more strongly than in case of the whole community-around intermediate values in late summer and fall.

Over the 12 years, morphological richness and divergence of Copepoda increased significantly while their concentration declined (electronic supplementary material, S7), as observed for the entire zooplankton community. There was no significant trend in morphological evenness (p > 0.05). Significant trends in morphological traits were more variable sizes with notably more abundant smaller individuals, an increase in variance of transparency, and increases in average and extreme elongation values. The slope of these trends was similar or lower than those found for the entire community.

Using morphological features derived from images, we successfully defined an interpretable morphological space comprising four synthetic traits: size, transparency, circularity and shape complexity. After defining objective groups of morphologically similar zooplankton, we detected both seasonal and long-term changes in morphological diversity indices computed from these groups. Over the spring-summer transition and over the 12-year study period, the total zooplankton concentration declined while the community became morphologically more diverse.

ZooScan imaging of preserved net samples allows the temporal decoupling of sampling and imaging, but comes with some limitations. First, the usual limitations of net-sampling still apply and the destruction of particularly fragile organisms can lead to an incomplete community description. In the target size range (200 µm to few centimetres), some Rhizaria and small jellies were likely missed. Second, the conservation of organisms with formalin can cause a temporary shrinkage of tissues in the days following fixation. Here, we scanned samples within weeks after collection (except for 2009 samples, which were processed in 2010) to minimize this bias and avoid further biases due to different storage periods. Finally, identifying organisms on images allows to process more of them than through manual microscopy but at the cost of a lower taxonomic resolution. In our study, taxa had to be aggregated to broader levels than what is possible to achieve from images to ensure consistency over the series. This limited the detection of taxonomic shifts that may have been visible at a finer taxonomic level. Nevertheless, it was sufficient to detect the seasonal shift in taxonomic diversity and an increase in taxa richness over the long term. A finer taxonomic resolution was available over the four most recent years. This was too short to perform a trend analysis, but the seasonal signal was very similar to that detected at the lower taxonomic resolution; only its amplitude was slightly affected (electronic supplementary material, S8).

On a seasonal scale, diversity and production were negatively correlated: all morphological and taxonomic diversity indices increased towards summer, while zooplankton concentration and nutrients declined. Besides a more diverse morphology, both absolutely (MRic) and relatively (MDiv), there was an increase in transparency and elongated shapes over the springsummer transition. Classically, the spring peak of zooplankton following the phytoplankton bloom denotes a period of high production. It is generally dominated by Copepoda, and taxonomic-and as we showed also morphological-diversity is low. With the decline of Copepoda and the increase of other taxa notably gelatinous Appendicularia and Chaetognatha [START_REF] Calbet | Annual Zooplankton Succession in Coastal NW Mediterranean Waters: The Importance of the Smaller Size Fractions[END_REF][START_REF] Feuilloley | Temporal fluctuations in zooplankton size, abundance, and taxonomic composition since 1995 in the North Western Mediterranean Sea[END_REF][START_REF] Fullgrabe | Zooplankton dynamics in a changing environment: a 13-year survey in the northwestern Mediterranean Sea[END_REF] diversity increased to stay on a high level later in the year. This increase in gelatinous plankton matches the morphological changes that we detected in the plankton community. Besides taxonomic changes, the ontogenetic development of organisms also likely increased morphological diversity towards summer: while copepodite stages of Copepoda are morphologically similar across species, their appearances as adults differ. This could explain the increase in average size and in the variability in morphological traits that we noticed (electronic supplementary material, S3). Thus, it seems that both increasing taxonomic diversity and the life cycle of dominating taxa caused the increase in morphological diversity towards the low-production, summer period.

Under an impoverishment in nutrients (especially nitrate) and organic matter, the zooplankton concentration declined over the 12-year period while the community became morphologically more diverse. This resembles the trend observed from spring to summer, albeit with a lower amplitude and likely caused by different ecological mechanisms.

Here also, the increase in morphological diversity could have resulted from increasing taxonomic diversity. This has been observed for functional diversity indices, in Copepoda and aquatic macroinvertebrates [START_REF] Becker | Latitudinal gradient of copepod functional diversity in the South Atlantic Ocean[END_REF][START_REF] Schmera | Functional diversity: a review of methodology and current knowledge in freshwater macroinvertebrate research[END_REF]. Also, the morphological changes (increases in transparency and elongation) correspond to the taxonomic shifts in the zooplankton community at point B in recent years [START_REF] Feuilloley | Temporal fluctuations in zooplankton size, abundance, and taxonomic composition since 1995 in the North Western Mediterranean Sea[END_REF]: after a drop in zooplankton concentration in 2015, following a sudden temperature increase in summer that year, gelatinous zooplankton (notably Appendicularia and Chaetognatha, both with rather elongated body However, while both taxonomic and morphological richness increased significantly, only morphological but not taxonomic diversity increased. This discrepancy may be partially attributable to the rather low taxonomic resolution in our study. Still, as discussed above, a higher resolution did not change the seasonal signal much. Also, when considering only Copepoda, a group morphologically more homogeneous but taxonomically diverse, we obtained morphological traits and trends in morphological diversity very similar to those of the whole community. Taxonomy therefore was unlikely the only driver. Finally, the fact that each morph comprised several taxa, and conversely each taxon comprising several morphs also indicates that morphology is not a mere mirroring of taxonomy. Overall, this underlines the relevance of considering diversity descriptors beyond taxonomy to detect changes in ecosystems from different perspectives [START_REF] Villéger | Functional homogenization exceeds taxonomic homogenization among European fish assemblages[END_REF].

In our case, a mechanism similar to that of taxonomic or functional specialization could also explain the observed pattern: diversity increases when a community shifts from generalists to more specialized taxa. Generalist or opportunistic taxa tolerate a large range of environmental conditions and/or express high growth rates, allowing them to establish their dominance in resource-rich conditions. One example are Copepoda (mostly Calanoida) during the spring bloom, but it has also been described for grassland moths or benthic macroinvertebrates under (agricultural) eutrophication [START_REF] Zhang | Nutrient enrichment homogenizes taxonomic and functional diversity of benthic macroinvertebrate assemblages in shallow lakes[END_REF][START_REF] Mangels | Diversity and trait composition of moths respond to land-use intensification in grasslands: generalists replace specialists[END_REF]. Such communities, dominated by generalists, often display low taxonomic and functional diversity [START_REF] Mangels | Diversity and trait composition of moths respond to land-use intensification in grasslands: generalists replace specialists[END_REF][START_REF] Villéger | Contrasting changes in taxonomic vs. functional diversity of tropical fish communities after habitat degradation[END_REF]. Under oligotrophic conditions, specialized taxa more efficiently exploit the scarce resources, each within a narrow range of conditions, and outcompete generalists [START_REF] Clavel | Worldwide decline of specialist species: toward a global functional homogenization?[END_REF]. In that situation, either a larger number of taxa coexist and taxonomic diversity increases, or a similar number of functionally more diverse taxa coexist and functional diversity increases, or both. Such mechanism has often been reported to cause taxonomic or functional homogenization with increasing eutrophication, along a spatial gradient of sites, for a range of aquatic [START_REF] Chihoub | Coastal eutrophication as a potential driver of functional homogenization of copepod species assemblages in the Mediterranean Sea[END_REF][START_REF] Li | Eutrophication and predation mediate zooplankton diversity and network structure[END_REF][START_REF] Jeppesen | Trophic structure, species richness and biodiversity in Danish lakes: changes along a phosphorus gradient[END_REF] and terrestrial [START_REF] Chisté | Intensive land use drives small-scale homogenization of plantand leafhopper communities and promotes generalists[END_REF] habitats. Benthic communities for example were more taxonomically heterogeneous under low nutrient and chlorophyll a levels, both within-and among lakes [START_REF] Donohue | Nutrient enrichment homogenizes lake benthic assemblages at local and regional scales[END_REF]. Our results suggest that this also happens in time and for morphological heterogeneity.

Over the study period, increases in temperature and salinity were accompanied by an impoverishment of nutrients, thereby intensifying the already oligotrophic conditions of the Mediterranean Sea. These observations match climate models predictions [START_REF] Adloff | Mediterranean Sea response to climate change in an ensemble of twenty first century scenarios[END_REF][START_REF] Pagès | Projected effects of climate-induced changes in hydrodynamics on the biogeochemistry of the Mediterranean sea under the RCP 8.5 Regional Climate Scenario[END_REF][START_REF] Soto-Navarro | Evolution of Mediterranean sea water properties under climate change scenarios in the Med-CORDEX ensemble[END_REF]. Increased stratification of the water column reduced the depth of vertical mixing in winter in recent years, thus preventing nutrient supply to the surface layer, as observed at two stations in the same region [START_REF] Margirier | Abrupt warming and salinification of intermediate waters interplays with decline of deep convection in the Northwestern Mediterranean Sea[END_REF]. Contrary to those expectations, chlorophyll a concentration, a proxy for autotrophic phytoplankton biomass, significantly increased. However, add a nearby site phytoplankton pigment measurements by high performance liquid chromatography revealed a concomitant shift of the phytoplankton community towards smaller taxa (Migon [Q1]pers. comm.), of lesser nutritious quality for zooplankton [START_REF] Jagadeesan | Copepod grazing and their impact on phytoplankton standing stock and production in a tropical coastal water during the different seasons[END_REF]. Thus, even though chlorophyll a was increasing, the amount of food actually available to the zooplankton community, notably to Copepoda, may have decreased.

In both scenarios (taxonomic changes or morphological specialization), the increase in morphological diversity can be explained by a decline of the dominant morphologies, here associated with Copepoda, both seasonally and over the long term. Twelve years of are not sufficient to disentangle inter-annual variability from climate-related effects. Still, the observed environmental trends were mostly consistent in a longer 20-year period (electronic supplementary material, S9), a time scale becoming relevant to study climate changes. Thus, the relationships between environmental forcing and zooplankton morphological diversity may extend further in time.

The observed increase in morphological diversity was caused by increasing variance in size and transparency as well as shifts towards a higher proportions of very small, transparent and elongated organisms, on both temporal scales. A full analysis of potential changes in ecosystem functioning requires detailed (e.g. species level) taxonomic identifications, to link each species to its functional traits and describe their change through time. This is not possible with our data, but the literature suggests some relationships between morphology and ecological functions that transcend species boundaries, as discussed in the following.

Body size was the most variable morphological trait in the zooplankton community (i.e. appearing as the first principal component). Although no significant change in average body size occurred, the proportion of very small organisms and the overall variation in size increased, both for the entire community and within Copepoda. This partially supports the hypothesis that warming would 'benefit the small' and is in line with previous works reporting or projecting increasing abundances of small taxa, from phyto-and zooplankton to fish [START_REF] Daufresne | Global warming benefits the small in aquatic ecosystems[END_REF][START_REF] Moullec | An end-to-end model reveals losers and winners in a warming Mediterranean Sea[END_REF][START_REF] Avaria-Llautureo | Historical warming consistently decreased size, dispersal and speciation rate of fish[END_REF][START_REF] Uszko | Smaller species but larger stages: warming effects on inter-and intraspecific community size structure[END_REF]. Although interspecific interactions are complex, large proportions of small body sizes, i.e. steep slopes of community size spectra, are usually associated with low transfer efficiency through the food web [START_REF] Atkinson | Increasing nutrient stress reduces the efficiency of energy transfer through planktonic size spectra[END_REF][START_REF] Toruan | Response of Zooplankton size structure to multiple stressors in urban lakes[END_REF]. Conversely, a wider size range of (meso-) zooplankton (through increased (c) Morphological shifts and possible functional implications size variability) might increase diet niche partitioning and increase trophic transfer [START_REF] García-Comas | Prey size diversity hinders biomass trophic transfer and predator size diversity promotes it in planktonic communities[END_REF][START_REF] Ye | Increasing zooplankton size diversity enhances the strength of top-down control on phytoplankton through diet niche partitioning[END_REF].

Transparency is an efficient camouflage in open water (Johnsen [START_REF] Johnsen | Hidden in Plain Sight: The Ecology and Physiology of Organismal Transparency[END_REF]) and a common strategy in plankton. With increased water clarity through oligotrophication transparency likely becomes more relevant to hide from visual predators [START_REF] Williamson | Toward a more comprehensive theory of zooplankton diel vertical migration: Integrating ultraviolet radiation and water transparency into the biotic paradigm[END_REF]. This is consistent with the increase in transparent organisms during a decrease in organic particles observed in our study, on both temporal scales. Here, one explanation of this increase, and of the trend towards more elongated body shapes, is the shift towards Appendicularia and Chaetognatha. Given the well-defined feeding strategies of both taxa (Chaetognatha are carnivores, Appendicularia are filter-feeders), those changes could hint at functional changes in the zooplankton community. Furthermore, increasing transparency of mesozooplankton might be negatively correlated to carbon export [START_REF] Perhirin | Morphological and taxonomic diversity of mesozooplankton is an important driver of carbon export fluxes in the ocean[END_REF], possibly because of similar taxonomic and functional changes in the community with consequences for the ocean carbon flux. This flux is affected by both zooplankton taxonomic composition and the morphology of marine snow particles [START_REF] Taucher | In situ camera observations reveal major role of zooplankton in modulating marine snow formation during an upwelling-induced plankton bloom[END_REF][START_REF] Trudnowska | Marine snow morphology illuminates the evolution of phytoplankton blooms and determines their subsequent vertical export[END_REF], making it interesting to investigate the relationship between zooplankton morphology and particle morphology and, in consequence, carbon export.

More generally, it would be relevant to investigate how these phenotypic variations as such relate to functional diversity and ecosystem functioning, e.g. whether increasing morphological diversity has the same beneficial effects as taxonomic or functional diversity on ecosystem functions such as nutrient transfer [START_REF] Ye | Functional diversity promotes phytoplankton resource use efficiency[END_REF], resource use efficiency [START_REF] Vallina | Phytoplankton functional diversity increases ecosystem productivity and stability[END_REF], primary or secondary production [START_REF] Setubal | Effects of functional diversity and salinization on zooplankton productivity: an experimental approach[END_REF][START_REF] Korhonen | Productivity-diversity relationships in lake plankton communities[END_REF]. A morphologically diverse community might, for example, harbour more specialized feeding strategies because of the large morphological range of prey, and changes in transparency and elongation might alter or indicate a modification of the trophic structure.

The four synthetic morphological traits we derived for a zooplankton community in the Mediterranean Sea were close to those previously identified for Copepoda in the Arctic [START_REF] Vilgrain | Trait-based approach using in situ copepod images reveals contrasting ecological patterns across an Arctic ice melt zone[END_REF] and phytoplankton in the NE Atlantic [START_REF] Sonnet | Length, width, shape regularity, and chain structure: time series analysis of phytoplankton morphology from imagery[END_REF]. They were not purely driven by ZooScan imaging or limited by the measurements performed on those images since (i) the two studies above used other instruments (the Underwater Vision Profile and the Imaging Flow Cyto Bot) and (ii) the four traits were similar for a morphological space built with features computed by a deep learning model. Therefore, they seem to be important phenotypic characteristics of planktonic communities, which strengthens the hypothesis that they have a functional role.

A quantitative imaging dataset of zooplankton allowed us to synthesize traits that provide insight into a widely unused facet of biological diversity. Our approach is generic and applicable to other taxonomic groups or ecosystems, provided that standardized images of individual organisms are available. To encourage its reuse, we share general-purpose code libraries. We demonstrated that common patterns observed and/or assumed for taxonomic and functional diversity, i.e. increased diversity through specialization under less productive conditions, were also evident in morphology and on different temporal scales: seasonally and over 12 years. To our knowledge, this is only the second time this supposedly wellknown link between production and diversity is demonstrated along a time series [START_REF] Abonyi | Functional diversity of phytoplankton highlights long-term gradual regime shift in the middle section of the Danube River due to global warming, human impacts and oligotrophication[END_REF]. Over the long term at least, the changes in morphology did not seem attributable to taxonomy only, suggesting a direct functional role of morphological traits. If these phenotypic changes are indeed caused by the impoverishment of Mediterranean waters, this trend likely continues in the future according to climate model predictions [START_REF] Cos | The Mediterranean climate change hotspot in the CMIP5 and CMIP6 projections[END_REF].

  Zooplankton and environmental variables were sampled at 'point B', near the entrance of the bay of Villefranche-sur-mer (France, North-Western Mediterranean Sea, 43°41′06.0″ N 7°18′56.4″ E). The bay is largely open to waters and plankton from the open sea since it is the tip of a 2000 m deep submarine canyon. Zooplankton was sampled in mornings and afternoons of each working day, from 1 July 2009 to 12 December 2020. Samples were taken with a WP2 net (200 µm mesh size), hauled vertically from 75 m depth to the surface (approx. 18.75 m of seawater filtered per haul)

Figure 1 .

 1 Figure 1.

3

 3 Seasonal component of environmental variables temperature, , , POC, PON), zooplankton and chlorophyll a concentration (a proxy for the biomass of autotrophic phytoplankton), taxonomic and morphological diversity indices (richness, evenness, diversity/divergence). This component was derived by Seasonal-Trend Decomposition by Loess over 2009-2020. Curves were smoothed with a moving average (window size = 5) and all variables were scaled and centred (mean = 0, variance = 1) for display.

Figure 3

 3 Figure 3.

  of environmental variables (temperature, POC, PON, , chlorophyll a) and zooplankton concentration, taxonomic and morphological diversity indices (richness, evenness, diversity/divergence) over the years 2009-2020. Black lines are the de-seasonalized series, red lines are the significant trends (with p < 0.05) obtained from GLS regression.

Figure 4 .

 4 Figure 4.

  diversity and specialization over the long term shapes) recovered, but crustacean densities remained rather low.

(c) Morphological space and definition of 'morphs' (d) Quantifying morphological and taxonomic diversity

(d) Intra-taxon analysis on Copepoda

Acknowledgements

The authors would like to thank the PIQv and CCPv platforms of EMBRC-France, a national Research Infrastructure supported by ANR, under the reference ANR-10-INSB-02, which handled the processing and storing of samples. They acknowledge the financial support to the zooplankton monitoring program by the host institution in Villefranche (previously Observatoire Océanologique de Villefranche-OOV, now Institut de la Mer de Villefranche-IMEV) and G Gorsky for initiating it. Sampling was carried out on the local vessel operated by the French Oceanographic Fleet. Environmental data and HPLC measurements were provided by the SNO SOMLIT; we thank C Migon for his advice on the HPLC data. Finally, we also thank two anonymous reviewers whose comments helped to improve the manuscript.

Funding

This work was funded mainly by our salaries as French state employees and therefore by French taxpayers. It was initially funded by the ModelOmics project of the émergence program of Sorbonne Université. Additional support was provided by the Institut des Sciences du Calcul et des Données (ISCD) of Sorbonne Université (SU) through the support of the sponsored junior team FORMAL (From ObseRving to Modelling oceAn Life), especially through the post-doctoral contract of MB. SDA also acknowledges the CNRS for her sabbatical year as visiting researcher at ISYEB in 2018-2020 and the support of the French Agence Nationale de la Recherche (ANR) under grant ANR-22-CE02-0023-1 (project TRAITZOO). Data sorting and accessibility was eased by the EcoTaxa application, the improvement of which was funded by the Belmont Forum through project WWW.PIC (ANR-18-BELM-0003) led by JOI. This project benefited from the support of HIRCOM funded by the SAM 'Société des explorations de Monaco'.

Data accessibility

Plankton occurrences and concentrations are available at https://doi.org/10.14284/473. The images and object level features are available from https://ecotaxa.obs-vlfr.fr; access will be provided after registration and upon request to the corresponding author. The environmental data are accessible under https://www.somlit.fr/demande-de-donnees. The code is available on GitHub https://github.com/jiho/ptb_morphodiv and Zenodo (10.5281/zenodo.10053905).

The data are provided in electronic supplementary material [START_REF] Beck | Morphological diversity increases with decreasing resources along a zooplankton time series[END_REF].

Ethics

This work did not require ethical approval from a human subject or animal welfare committee.

Declaration of AI use

We have not used AI-assisted technologies in creating this article.

Authors' contributions

M.B.: formal analysis, investigation, methodology, writing-original draft, writing-review and editing; C.C.: conceptualization, formal analysis, investigation, writing-review and editing; L.G.: investigation, writing-review and editing; C.D.: data curation, writing-review and editing; L.J.: data curation, writing-review and editing; A.E.: data curation, writing-review and editing; L.S.: investigation, writing-review and editing; S.A.: conceptualization, funding acquisition, investigation, methodology, writing-review and editing; J.I.: conceptualization, formal analysis, funding acquisition, investigation, methodology, writing-review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.


Nicolas Faroß
email: faross@math.uni-sb.de

Simon Schwarz
email: sschwarz@mpi-inf.mpg.de

Gröbner Bases for Boolean Function Minimization

Keywords: Multi-level Logic Optimization, Boolean Function Synthesis, Gröbner Bases

Boolean function minimization techniques try to find, for a given formula, a smaller equivalent formula. In this work, we present a novel technique for heuristic boolean function minimization. By using an algebraic encoding, we embed the minimization problem into an algebraic domain, where algorithms for computing Gröbner bases are applicable. A Gröbner basis usually forms a compact representation of our encoded function. From the Gröbner basis, we then reconstruct an equivalent, more compact boolean formula. Our approach is the first to use Gröbner bases for function minimization. Combined with advances of algebraic Gröbner bases in satisfiability checking, this motivates further research on applications of Gröbner bases in the context of boolean logic.

Introduction

Boolean function minimization algorithms find, given a boolean formula, a small equivalent boolean formula. This, for example, is a central problem in circuit synthesis [START_REF] Brayton | Logic minimization algorithms for VLSI synthesis[END_REF][START_REF] Rudell | Multiple-valued logic minimization for PLA synthesis[END_REF].

In general, the exact boolean minimization problem is NP-hard, it even is Σ 𝑃 2 -complete [START_REF] Buchfuhrer | The complexity of boolean formula minimization[END_REF]. Hence, exact approaches quickly become infeasible. Thus, approaches in this area often have trade-offs between result quality and runtime. Previous work often focused on minimizing disjunctive normal forms (DNF) or used syntactical heuristics to simplify functions (Section 1.1).

Our contribution is a novel heuristic approach for multi-level boolean function minimization. Given an input formula in DNF, our approach encodes it as a system of equations over F 2 , using a modification of previous encodings [START_REF] Condrat | A Gröbner basis approach to CNF-formulae preprocessing[END_REF]. This representation allows computing a Gröbner basis with common algorithms [START_REF] Buchberger | A theoretical basis for the reduction of polynomials to canonical forms[END_REF]. The resulting Gröbner basis, then, can again be interpreted as an empirically smaller, equivalent boolean formula. This process is repeated recursively to further minimize the formula. Contrary to other minimization algorithms, our encoding over F 2 can efficiently represent the exclusive-or (⊕) operator, allowing exponentially more compact representations compared to DNF. Thus, our approach works particularly well for formulas which contain many exclusive-or operations. In practice, such formulas occur, for example, in microarchitectural hash functions [START_REF] Mccalpin | Mapping addresses to L3/CHA slices in Intel processors[END_REF]. Some proprietary hash functions, e.g. used for cache slice indexing in modern microprocessors have been reversed (i.e. a compact representation has been found) by using our approach (work under submission). In practice, we observe more compact representations than a minimized DNF even on random formulas, Section 3. Furthermore, for some classes of formulas, nearly optimal results are provable (Theorem 2).

In addition to previous applications of Gröbner bases in logic [START_REF] Condrat | A Gröbner basis approach to CNF-formulae preprocessing[END_REF][START_REF] Nguyen | Combinations of boolean Gröbner bases and SAT solvers[END_REF] that focus on preprocessing for satisfiability checking, our approach is a novel, complementary application of Gröbner bases for boolean logic. This motivates further research in this direction.

Related Work

Boolean Minimization Usually, a boolean logic minimizer is given a formula in DNF, which can be easily obtained from a truth table. Then, it produces a small, equivalent formula. Existing tools can be classified by the shape of the output formula:

Two-level logic optimization tools produce again a formula in DNF. Classically, Quine and McCluskey [START_REF] Quine | The problem of simplifying truth functions[END_REF][START_REF] Mccluskey | Minimization of boolean functions[END_REF] focus on optimal two-level minimization. However, their optimal approach quickly becomes infeasible [START_REF] Brayton | Logic minimization algorithms for VLSI synthesis[END_REF]. Thus, later approaches such as ESPRESSO [START_REF] Brayton | Logic minimization algorithms for VLSI synthesis[END_REF] make use of a heuristic search for producing a compact DNF. Still, representing a formula in DNF can have exponential overhead. For example, any DNF of 𝑓 (𝑥 1 , . . . , 𝑥 𝑛 ) = 𝑥 1 ⊕ • • • ⊕ 𝑥 𝑛 is exponentially larger in 𝑛 than 𝑓 . This motivates more general minimization approaches.

Multi-level logic optimization allows minimized formulas of arbitrary shape. Still, some tools produce only formulas of a fixed structure. For example, EXORCISM [START_REF] Mishchenko | Fast heuristic minimization of exclusive-sums-of-products[END_REF] always produces formulas of the shape "exclusive-sum-of-products". Other tools produce formulas of arbitrary depth and structure, for example MIS [START_REF] Brayton | MIS: A multiple-level logic optimization system[END_REF] or LSS [START_REF] Darringer | LSS: A system for production logic synthesis[END_REF]. Most multi-level optimization tools rely on two-level optimization and apply syntactic transformations such as subexpression recognition and replacement to the resulting formula [START_REF] Perkowski | Multi-level logic minimization[END_REF]. More recently, approaches for multi-level synthesis based on satisfiability [START_REF] Soeken | Practical exact synthesis[END_REF][START_REF] Haaswijk | SAT-based exact synthesis: Encodings, topology families, and parallelism[END_REF] have been studied. They provide optimal solutions but are only feasible for small instances. Our presented approach provides a heuristic method for multi-level logic optimization which, contrary to previous work, does not rely on syntactic recognitions.

Gröbner Bases are special generating sets for ideals and are widely used in computer algebra [START_REF] Becker | Gröbner bases[END_REF], for example for ideal membership testing of a polynomial. However, we focus on applications of Gröbner bases in logic. For example, Gröbner bases are used in the context of satisfiability checking and model counting, as well as in verification and SMT solving.

In satisfiability checking, Gröbner bases can be used for pre-processing clause sets in conjunctive normal form (CNF) [START_REF] Condrat | A Gröbner basis approach to CNF-formulae preprocessing[END_REF][START_REF] Nguyen | Combinations of boolean Gröbner bases and SAT solvers[END_REF]. Concretely, it is possible to encode a (sub-)set of clauses of a CNF formula as a system of polynomials. Then, a Gröbner basis for this system is computed. The resulting system of polynomials is then, again, interpreted as a set of clauses. The resulting, usually more compact set of clauses is equivalent to the original set. Hence, it is possible to replace sets of clauses with more compact, equivalent sets. Satisfiability checking on the preprocessed clauses is usually faster than on the original CNF [START_REF] Condrat | A Gröbner basis approach to CNF-formulae preprocessing[END_REF]. However, the cost of computing the Gröbner bases outweighs the benefits. A similar encoding can be used for model counting for boolean formulas [START_REF] Gao | Counting zeros over finite fields with Gröbner bases[END_REF]. Given a set of clauses reduced by the above technique, it becomes (computationally) easy to count satisfiable assignments. In contrast to the above approaches, we apply Gröbner bases to formulas in disjunctive normal form (DNF), and in the context of logic minimization.

Other applications of Gröbner basis are verification of arithmetic gates [START_REF] Ritirc | Improving and extending the algebraic approach for verifying gate-level multipliers[END_REF][START_REF] Kaufmann | Verifying large multipliers by combining SAT and computer algebra[END_REF][START_REF] Kaufmann | Incremental column-wise verification of arithmetic circuits using computer algebra[END_REF] and SMT solving over finite fields [START_REF] Hader | SMT solving over finite field arithmetic[END_REF][START_REF] Ozdemir | Satisfiability modulo finite fields[END_REF] or real numbers [START_REF] Junges | On Gröbner bases in the context of satisfiabilitymodulo-theories solving over the real numbers[END_REF]. However, note that both applications do not directly encode boolean formulas as polynomials, but use Gröbner bases for theory solving. Still, combined with our work, this suggests that Gröbner bases have multiple applications in logic and could be worthwhile to investigate further.

Gröbner Bases for Logic Minimization

Encoding Boolean Formulas as Polynomials

In the following, we will identify the truth-values {false, true} with the set F 2 = {0, 1}, where false corresponds to 0 and true corresponds to 1. Note that F 2 has the structure of a field with multiplication given by ∧ and addition is given by ⊕. Hence, we can view a polynomial

𝑓 = ∑ 𝑖 ∏ 𝑗 𝑦 𝑖𝑗 ∈ F 2 [𝑥 1 , . . . ,

⋀

𝑗 𝑦 𝑖𝑗 . Conversely, we can encode any boolean formula 𝜙 as a polynomial over F 2 by additionally using the identities ¬𝑥 = 𝑥 ⊕ 1 and 𝑥 ∨ 𝑦 = ¬(¬𝑥 ∧ ¬𝑦), which will be encoded as 𝑥 + 1 and 1 + (1 + 𝑥) • (1 + 𝑦) = 𝑥𝑦 + 𝑥 + 𝑦 respectively. We call this the algebraic encoding of 𝜙.

The previous relationship between boolean formulas and polynomials can be further extended to a correspondence between equivalent formulas and ideals in F , where we identify a literal 𝑦 𝑖𝑗 with a polynomial via our algebraic encoding. This can be verified by using Theorem 1 and observing that 𝜙(𝑥) = 0 if and only if 𝑓 𝑖 (𝑥) = 0 for all 𝑓 𝑖 , which is equivalent to 𝑓 (𝑥) = 0 for all 𝑓 ∈ 𝐼. Similarly, an ideal 𝐼 = (𝑓 1 , . . . , 𝑓 𝑚 , 𝑥 2 1 + 𝑥 1 , . . . , 𝑥 2 𝑛 + 𝑥 𝑛 ) corresponds to the boolean formula ⋁ 𝑖 𝑓 𝑖 , where the polynomials 𝑓 𝑖 are again identified with a boolean formula. Note that the polynomials 𝑥 2 𝑖 + 𝑥 𝑖 represent the idempotency law 𝑥 𝑖 ∧ 𝑥 𝑖 = 𝑥 𝑖 and allow the elimination of all higher powers of 𝑥 𝑖 during the Gröbner basis computation. Thus, there exist only 2 𝑛 different leading monomials besides 𝑥 2 𝑖 , which implies that the size of a reduced Gröbner basis is bounded by 2 𝑛 + 𝑛. Further, it is possible to compute this Gröbner basis in time 2 𝑂(𝑛) using Buchberger's algorithm, see Proposition 4.1.1 in [START_REF] Gao | Counting zeros over finite fields with Gröbner bases[END_REF]. This provides a significant improvement compared to the usual double exponential bound for Gröbner bases.

Gröbner Basis Minimization

Empirically, Gröbner bases are relatively small if the ideal is not too complex. Hence, we simplify a boolean formula by converting it to an ideal 𝐼. Then, we compute a Gröbner basis of 𝐼 to obtain a set of generators 𝐺. 𝐺 is then converted back to a formula of the form ⋁ 𝑔∈𝐺 𝑔 ∼ = ⋁ 𝑔∈𝐺 ¬¬𝑔. We can apply this approach recursively to the terms ¬𝑔 to further reduce the size and obtain a formula of the form

⋁ 𝑖 ¬ ⋁ 𝑗 ¬ ⋁ 𝑘 . . . ∼ = ⋁ 𝑖 ⋀ 𝑗 ⋁ 𝑘 .
. ., as depicted in the following algorithm:

1: function Minimize(𝜙) 2:
𝐺 ← GröbnerBasis(Ideal(𝜙))

3:

for 𝑔 ∈ 𝐺 do 4:
if 𝑔 is linear or recursion limit is reached then 5:

𝜙 𝑔 ← Formula(𝑔) 6:

else 7:
𝜙 𝑔 ← ¬Minimize(DNF(¬Formula(𝑔)))

8:

return ⋁ 𝑔∈𝐺 𝜙 𝑔 Note that it is not guaranteed that our algorithm produces a smaller formula and by the previous bound the Gröbner bases can have a size exponential in the number of variables. However, empirically we observed very compact formulas. Furthermore, for special classes of formulas, better bounds can be proven.

Theorem 2.

Let 𝜙 be a formula in DNF which is equivalent to the be matrix given by the non-zero rows in the reduced row echelon form of 𝐴, which can be computed using Gaussian elimination. Then one can check that a reduced Gröbner basis of 𝐼 is given by the terms ∑ 𝑛 𝑗=1 𝐴 ′ 𝑖𝑗 𝑥 𝑖 corresponding to the rows of 𝐴 ′ and the terms 𝑥 2 𝑖 + 𝑥 𝑖 for which 𝑥 𝑖 is not a leading monomial of one of the first terms. Since a reduced Gröbner basis is unique, it will also be computed by our minimization algorithm. Further, the terms of the form 𝑥 2 𝑖 + 𝑥 𝑖 are redundant and will be removed, such that a formula of the form ⋁ 𝑚 𝑖=1 ⨁ 𝑛 𝑗=1 𝐴 ′ 𝑖𝑗 𝑥 𝑗 is returned. By counting the number of operators and variables, we obtain that the size of such a formula is always less than 2𝑛𝑚.

⋁ 𝑖 ⨁ 𝑗 𝐴 𝑖𝑗 𝑥 𝑗 for a matrix 𝐴 ∈ F 𝑛×𝑛 2 . Then |Minimize(𝜙)| < 2 • 𝑛 • rank(𝐴),
Further Improvements In addition to our main algorithm, we implemented the following two improvements: First, minimizing the DNF with a two-level minimization algorithm during preprocessing reduces the size of the input. This does not change result quality, but can significantly speed up the algebraic computations. Second, some elements of the Gröbner basis can be redundant in the sense that a subset of the Gröbner basis already generates the corresponding ideal. Finding a minimal generating subset of the Gröbner basis can be formulated as a weighted set-cover problem, which we solve in a post-processing step to improve result quality. Eliminating these redundant terms directly during the computation of the Gröbner basis is subject to further work.

Evaluation

A first prototype of this algorithm has been implemented in SageMath [START_REF]The Sage Developers[END_REF], using SINGULAR [START_REF] Decker | Singular -A computer algebra system for polynomial computations[END_REF] for computing Gröbner bases, as well as ESPRESSO [START_REF] Brayton | Logic minimization algorithms for VLSI synthesis[END_REF] for DNF preprocessing and GLPK [START_REF][END_REF] for set-cover computations. Note that Gröbner bases are computed with respect to the degree reverse lexicographic order. A basic benchmark in Figure 1 shows the average size of a reconstructed formula, and the runtime of the minimization process. This is compared to a minimized formula in DNF from ESPRESSO [START_REF] Brayton | Logic minimization algorithms for VLSI synthesis[END_REF] and a minimal formula produced with ABC's 'exact' command [START_REF] Brayton | ABC: An academic industrial-strength verification tool[END_REF][START_REF] Soeken | Practical exact synthesis[END_REF]. The benchmark ran on 200 random formulas for each size. The formulas were sampled uniformly at random on a syntactic level, i.e. one syntax tree is sampled uniformly at random out of all possible syntax trees for a specific size. Note that exact, optimal synthesis with ABC is only possible for up to 7-8 variables, while our approach can deal with up to 20 variables. Our resulting formulas are always more compact than a DNF produced by ESPRESSO. However, our approach has a significantly greater runtime. We are currently working on reverse-engineering microarchitectural hash functions used e.g. for cache indexing in processors with our algorithm (to appear in [START_REF] Gerlach | Efficient and generic microarchitectural hash-function recovery[END_REF]).

Conclusion

We have presented a novel technique for heuristic multi-level formula minimization based on an algebraic encoding of formulas and Gröbner basis computations. Currently, our algorithm can handle functions with up to 20 input bits, whereas exact synthesis can only handle up to 8 input bits. Empirically, our approach produces smaller formulas than two-level minimization algorithms such as ESPRESSO. In particular, our idea works well for formulas with many exclusive-or operations. Our algorithm shows a new application of Gröbner basis computations in logic minimization. In addition to previous work on Gröbner bases in satisfiability [START_REF] Condrat | A Gröbner basis approach to CNF-formulae preprocessing[END_REF][START_REF] Nguyen | Combinations of boolean Gröbner bases and SAT solvers[END_REF], this result suggests that the use of our algebraic encoding, as well as Gröbner bases, brings benefits to boolean logics and could be worthwhile to investigate further.

Figure 1 :

 1 Figure 1: Comparison of output size and runtime.

  Given a boolean DNF 𝜙 = ⋁ 𝑖 ⋀ 𝑗 𝑦 𝑖𝑗 with literals 𝑦 𝑖𝑗 ∈ {𝑥 1 , . . . , 𝑥 𝑛 , ¬𝑥 1 , . . . , ¬𝑥 𝑛 }, we want to find an equivalent smaller formula. Our idea is to switch from boolean formulas 𝜙 to ideals 𝐼 in the polynomial ring F 2 [𝑥 1 , . . . , 𝑥 𝑛 ]. Then, Gröbner bases of 𝐼 correspond empirically to smaller equivalent formulas of 𝜙.

  𝑥 𝑛 ] with 𝑦 𝑖𝑗 ∈ {1, 𝑥 1 , . . . , 𝑥 𝑛 } as a boolean formula of the form ⨁

	𝑖

  2 [𝑥 1 , . . . , 𝑥 𝑛 ]. Recall that an ideal 𝐼 in a polynomial ring 𝑅 is a subset 𝐼 ⊆ 𝑅 which can be written as 𝐼 = { ∑ 𝑚 𝑖=1 𝑔 𝑖 𝑓 𝑖 | 𝑔 1 , . . . , 𝑔 𝑚 ∈ 𝑅} for some polynomials 𝑓 1 , . . . , 𝑓 𝑚 ∈ 𝑅. Such polynomials are called generators of 𝐼, and we write 𝐼 = (𝑓 1 , . . . , 𝑓 𝑚 ). There is a bijection between equivalence classes of boolean formulas [𝜙] and ideals 𝐼 ⊆ F 2 [𝑥 1 , . . . , 𝑥 𝑛 ] containing 𝑥 2 1 +𝑥 1 , . . ., 𝑥 2 Hence, there is a bijection between equivalent formulas and subsets of F 𝑛 2 . Similarly, any ideal containing 𝑥 2 1 + 𝑥 1 , . . . , 𝑥 2 𝑛 + 𝑥 𝑛 is uniquely determined by its zero set {𝑥 ∈ F 𝑛 2 | ∀𝑓 ∈ 𝐼, 𝑓 (𝑥) = 0} by Hilbert's Nullstellensatz for finite fields [24]. Additionally, every subset 𝑌 ⊆ F 𝑛 2 is finite and hence a zero set of an ideal by elementary results from algebraic geometry. Further, this ideal can always be assumed to contain 𝑥 2 1 + 𝑥 1 , . . . , 𝑥 2 𝑛 + 𝑥 𝑛 , such that we have a bijection between ideals containing 𝑥 2 1 + 𝑥 1 , . . . , 𝑥 2 𝑛 + 𝑥 𝑛 and subsets 𝑌 ⊆ F 𝑛 2 . By combining this bijection which the first one, we obtain the statement of the theorem. Importantly, we can explicitly convert a boolean formula to generators of the corresponding ideal and vice versa. Let 𝜙 = ⋁ 𝑖 ⋀ 𝑗 𝑦 𝑖𝑗 be a formula in DNF, then its corresponding ideal 𝐼 is generated by 𝑥 2 1 + 𝑥 1 , . . ., 𝑥 2 𝑛 + 𝑥 𝑛 and 𝑓 𝑖 = ∏ 𝑗 𝑦 𝑖𝑗

	Theorem 1.

𝑛 +𝑥 𝑛 , such that 𝜙(𝑥) = 0 ⇐⇒ ∀𝑓 ∈ 𝐼, 𝑓 (𝑥) = 0. Proof. Let [𝜙] be an equivalence class of boolean formulas in 𝑛 variables. Then it is uniquely determined by its zero set {𝑥 ∈ F 𝑛 2 | 𝜙(𝑥) = 0}. Conversely, any set 𝑌 ⊆ F 𝑛 2 is the zero set of the formula ⋀ 𝑦∈𝑌 ⋁ 𝑛 𝑖=1 (𝑥 𝑖 ⊕ 𝑦 𝑖 ).

  where | • | denotes the number of boolean operators and variables in a formula. . Then the corresponding ideal 𝐼 from Theorem 1 is generated by 𝑥 2 𝑖 + 𝑥 𝑖 and the linear terms ∑ 𝑛 𝑗=1 𝐴 𝑖𝑗 𝑥 𝑖 defined by the rows of 𝐴. Denote with 𝑚 the rank of 𝐴 and let 𝐴 ′ ∈ F 𝑚×𝑛

	Proof. Consider a boolean formula 𝜙 which is equivalent to	⋁ 𝑛 𝑖=1	⨁ 𝑛 𝑗=1 𝐴 𝑖𝑗 𝑥 𝑗 for a matrix
	𝐴 ∈ F 𝑛×𝑛 2
	2

Acknowledgements: We thank our anonymous reviewers for their constructive comments.