Geometric Morphometrics of Bilateral Asymmetry in Eunotia bilunaris (Eunotiales, Bacillariophyceae) as a Tool for the Quantitative Assessment of Teratogenic Deviations in Frustule Shapes

: A number of pennate diatom genera typically have teratogenic deformations of their siliceous frustules due to the effects of environmental stress, such as high concentrations of heavy metals and low pH. However, the quantitative assessment of these deformations has rarely been applied. One species in which aberrations have frequently been reported is Eunotia bilunaris , which typically has bilaterally symmetric frustules with dorso-ventral differentiation. In this study, we aimed to illustrate the geometric morphometric analysis of symmetry as a tool for assessing the severity of teratogenic deformations. These were quantiﬁed by Procrustes superimposition of equidistant points placed along the valvar outlines in pairs of conﬁgurations based on their bilateral reﬂection symmetry. The shape deformations were mostly conﬁned to central parts of the ventral outlines and were captured both by the symmetric and asymmetric subspaces of the variation. The amount of bilateral asymmetry in individual cells was negatively related to frustule size via the allometric power law relationship, illustrating that asymmetry increased in the asexual diminution series. The presented analysis provides a framework for the quantitative assessment of frustule deformations in eunotioid diatoms that can be used for the comparative scoring of teratogenic deviations among cells, populations, or species.


Introduction
Diatoms are one of the most significant groups of protist organisms. They abound in phytoplankton and phytobenthos of aquatic habitats worldwide. It is estimated that they account for as much as 20% of primary production on Earth [1]. Their most important phenotypic synapomorphy is the rigid silica cell wall (frustule), which is composed of two complementary parts (thecae). The upper part of frustules, called the valve, is characterised by various types of symmetric patterns, such as the radial symmetry of the centric diatoms or the biradial and bilateral symmetry of the pennate diatoms [2]. Empty diatom frustules are well preserved in sediments. Therefore, they are a key model group in paleoecological research of aquatic ecosystems [3].
Diatoms are generally sensitive to the environmental conditions of their habitats, which is reflected in the changes of their species composition in natural communities and in the patterns of morphological variation of frustules [4,5]. Teratogenic changes in frustule morphology have been well documented in scientific literature, especially in relation to various stress factors, such as increased sublethal concentrations of heavy metals [6][7][8], acidification [9,10], increased radiocativity [11], and organic pollution [12]. These teratogenic changes affect different parts of frustules, but the most common are deformities of the frustule outlines compared to their species-specific symmetric shapes.
The incidence of such morphological deformations in the sedimentary record has frequently been used as proxy data for the identification of historical changes in the anthropogenic loading of various aquatic ecosystems [3,5,6,13]. In these studies, the stress rate is usually quantified as the proportion of deformed frustules relative to the "standard" shapes. Thus, the severity of deformation of individual frustules is usually not explicitly quantified in any way. However, in recent years, several studies have been published that have sought to quantify the outline deformations of diatom frustules using the techniques of geometric morphometrics [14][15][16]. In general, these studies have shown that the degree of teratogenic deviations can be sensitively reflected in the patterns of shape changes of the studied populations. However, they did not separate symmetric and asymmetric variation, but rather focused on the principal components of the overall shape dynamics within the studied datasets that typically correlate with the degree of teratogenic shape deviations. Thus, to the best of our knowledge, no study has yet explicitly quantified shape asymmetries between the corresponding parts of the configurations representing the frustule shapes with teratogenic morphology.
There is, however, a well-founded analytical apparatus based on the geometric morphometric analysis of the configurations, representing a finite symmetry group of the studied system [17,18]. Principal component analysis (PCA) of such Procrustes-aligned configurations separates symmetric variation among objects from morphometric asymmetry between individual symmetric parts of the studied system. In addition, the Procrustesaligned coordinates of the configurations, representing the symmetry group of a single object, can be used to quantify the shape asymmetry of that particular specimen [18,19]. Thus, an obvious advantage of such an analysis is the ability to quantify the deformation of each object separately, i.e., without considering the variation of other individuals. With this analytical algorithm, it might be possible to compare the severity of teratogenic deformation of individual frustules without having to re-analyse entire large datasets each time.
From an analytical point of view, however, it is necessary to apply a proper symmetry group corresponding to the geometric arrangement of analysed frustules for the morphometric analysis of symmetry and asymmetry in different diatom lineages. The pennate diatoms, which usually dominate the freshwater phytobenthos, are usually characterised by different types of two-dimensional biradial and bilateral symmetry in the valvar view of their frustules. For example, biradial symmetry, characterised by two perpendicular axes of symmetry dividing the valve into four ideally identical segments, is typical for the genera Frustulia, Navicula, Luticola, and Achnanthidium [20][21][22][23]. However, even in these biradial genera, their ideally symmetric arrangement is usually broken by slight systematic asymmetries of valve morphology, such as the structure of the raphe system or Voigt discontinuities [24][25][26]. In addition, the bilaterally symmetric valve patterns repeatedly evolved in several pennate diatom lineages [27]. One such evolutionarily significant example of bilaterally symmetric morphology is the eunotioid diatoms (Eunotiales), which are characterised by valves with a distinct asymmetric curvature across the apical axis. Thus, in the valvar view, their cells are distinctly dorsiventral [2]. In contrast, the symmetry across the transapical axis of eunotioid valves is typically preserved, and two cellular halves cannot be distinguished from each other among different individuals. Members of the genus Eunotia Ehrenberg, the most frequently distributed genus of this lineage, occur in a wide variety of habitats, including the freshwater phytobenthos and subaerial biofilms and soils [2,[9][10][11]28]. In this study, the eunotioid diatom morphotype was selected to explore the potential of geometric morphometric analysis of asymmetry in assessing the severity of teratogenic deformities. We aimed to identify an analysis that could easily be replicated with different datasets for the comparative analysis of asymmetry in the shapes of eunotioid diatom frustules.
The model system of our study was Eunotia bilunaris (Ehrenberg) Schaarschmidt, a diatom species complex occurring mostly in oligotrophic waters with a low pH, such as peatland aquatic habitats in boreal ecosystems. It is composed of several closely related and morphologically very similar lineages. The valves of E. bilunaris are typical for arcuate isopolar outlines with rounded polar ends that are slightly tapered ( Figure 1). Interestingly, in the specimens of this species, teratogenic deformations of valvar outlines have been repeatedly reported [3,6,7,9,10,13]. These deformities were largely linked to adverse environmental conditions, such as heavy metal pollution [6,7] or low pH [9,10]. Our study was based on the analysis of cells from a natural population of E. bilunaris with a relatively high proportion of visibly deformed cells.
Symmetry 2021, 13, x FOR PEER REVIEW 3 of 11 The model system of our study was Eunotia bilunaris (Ehrenberg) Schaarschmidt, a diatom species complex occurring mostly in oligotrophic waters with a low pH, such as peatland aquatic habitats in boreal ecosystems. It is composed of several closely related and morphologically very similar lineages. The valves of E. bilunaris are typical for arcuate isopolar outlines with rounded polar ends that are slightly tapered ( Figure 1). Interestingly, in the specimens of this species, teratogenic deformations of valvar outlines have been repeatedly reported [3,6,7,9,10,13]. These deformities were largely linked to adverse environmental conditions, such as heavy metal pollution [6,7] or low pH [9,10]. Our study was based on the analysis of cells from a natural population of E. bilunaris with a relatively high proportion of visibly deformed cells. Decomposition of the symmetric and asymmetric components was performed to illustrate the main patterns of shape variation in both these subspaces. Our main goal was to determine whether the different levels of frustule shape deformations could be detected and quantified by geometric morphometric analysis of their bilateral asymmetry. We also determined whether the shape asymmetry of frustules was related to their size, i.e., if the population was constrained by an allometric relationship that could be linked to the changes in cell size during a series of asexual cell divisions.

Sampling and Data Acquisition
Algae were sampled on 21 January 2021, in a Sphagnum-dominated peat bog ditch located in the mountainous plateau of Krušné hory/Erzgebirge Mts. in Northwest Bohemia, Czech Republic (GPS: 50.536616N, 13.251838E) at an altitude of 827 m above sea level. The pH level of the sampling locality was 4.1, and the conductivity of the water was 112 μS·cm −1 . Immediately after sampling, the diatoms were fixed with ethanol so that the final solution had a 10% concentration of fixative. The valves were cleaned by oxidation with 35% hydrogen peroxide and potassium dichromate for 15 min. Then, they were washed with distilled water and mounted onto permanent slides using Naphrax (Brunel Microscopes Ltd., Wiltshire, UK).
Two-dimensional images of valves were taken on a Leica DM2500 light microscope (Leica Microsystems, Wetzlar, Germany) at 1000× magnification using Bresser MikroCam Decomposition of the symmetric and asymmetric components was performed to illustrate the main patterns of shape variation in both these subspaces. Our main goal was to determine whether the different levels of frustule shape deformations could be detected and quantified by geometric morphometric analysis of their bilateral asymmetry. We also determined whether the shape asymmetry of frustules was related to their size, i.e., if the population was constrained by an allometric relationship that could be linked to the changes in cell size during a series of asexual cell divisions.

Sampling and Data Acquisition
Algae were sampled on 21 January 2021, in a Sphagnum-dominated peat bog ditch located in the mountainous plateau of Krušné hory/Erzgebirge Mts. in Northwest Bohemia, Czech Republic (GPS: 50.536616N, 13.251838E) at an altitude of 827 m above sea level. The pH level of the sampling locality was 4.1, and the conductivity of the water was 112 µS·cm −1 . Immediately after sampling, the diatoms were fixed with ethanol so that the final solution had a 10% concentration of fixative. The valves were cleaned by oxidation with 35% hydrogen peroxide and potassium dichromate for 15 min. Then, they were washed with distilled water and mounted onto permanent slides using Naphrax (Brunel Microscopes Ltd., Wiltshire, UK).
Two-dimensional images of valves were taken on a Leica DM2500 light microscope (Leica Microsystems, Wetzlar, Germany) at 1000× magnification using Bresser Mikro-Cam 5.0 MP (Bresser, Rhede, Germany) digital photographic equipment. Micrographs representing a total of 64 valves were randomly selected, cropped, and rotated in GIMP ver. 2.8.14 [29], with their apical axis horizontal and the dorsal edge of the valves pointing downwards. The outline forms were digitised in TpsUtil, ver. 2.22 [30]. Two fixed landmarks were placed at the apices of each valve. The ventral and dorsal outlines were registered by 61 semi-landmarks placed along each valve side using the semi-automated background_curves tool of TpsDig, ver. 2.22. The equidistant position of these points along the ventral and dorsal valve outlines was achieved using the function digit.curves of the geomorph package, ver. 4.0.0, in R, ver. 4.0.5 [31]. This procedure decreased the number of points in each outline curve by one [32]. Thus, in the final dataset, each valve was registered by 122 two-dimensional outline points. These primary data are available online at https://doi.org/10.5281/zenodo.5732459 (accessed on 16 December 2021).

Geometric Morphometrics
The valvar outline of an eunotioid diatom, registered by a series of two-dimensional points, represents a configuration with bilateral object symmetry [33]. Such a configuration consists of two halves that can be reflected across the axis of symmetry, ideally producing two identical configurations. Any difference between the original and the reflected configuration is due to the presence of asymmetry in the position of the individual points. However, this asymmetry cannot be decomposed into directional and fluctuating components, as in the standard procedure of a two-factor, mixed-effect analysis of variance (ANOVA) model that is frequently used in geometric morphometric studies of bilateral symmetry in multicellular organisms [18]. Consequently, fluctuating asymmetry could not be used in these diatoms as an explicit measure of individual stress or developmental instability among specimens [34] because frustules of eunotioid diatoms typically lack any unambiguous directional orientation of their two symmetrical halves delimited by the transapical axis. In a set of such objects, the individual halves cannot be unambiguously designated as "right" or "left." This situation is not at all exceptional and is even typical of most unicellular protist organisms [35][36][37][38][39][40]. The analysis of a dataset must then be limited to the decomposition of the overall asymmetry and symmetric differences among superimposed objects [17].
In geometric morphometrics, the optimal superimposition of individual configurations is typically achieved using generalised Procrustes analysis (GPA), which optimises the relative position of corresponding points by removing any differences caused by their different position, rotation, and size from the set of configurations [41]. This is accomplished by rotating each configuration to minimise the summed squared distances of landmarks from an iterated estimate of the mean shape. Principal component analysis (PCA) of such Procrustes-aligned coordinates obtained by the GPA of the original and reflected configurations of all objects in the studied dataset yields a set of principal components (PCs) that each describe either purely symmetric shape variation among objects or patterns of bilateral shape asymmetry between their two corresponding halves [17,18]. The ratio of shape variation spanned by the asymmetric and purely symmetric PCs indicates the overall shape asymmetry of the analysed dataset.
Alternatively, pairs of configurations representing a symmetry group of a single eunotioid object, i.e., the original and reflected configurations of a single frustule, can be superimposed separately [18,19]. The residual distances between corresponding landmarks then represent the overall shape asymmetry between the two halves of one object (frustule). An important feature of such a procedure is that it is independent of the other analysed objects, and each individual frustule can be characterised by its own indicator of shape asymmetry [19]. However, it is not possible to visualise the patterns of symmetric and asymmetric variation within the explored morphospace, which is an inherent property of PCA based on the Procrustes aligned coordinates representing the complete symmetry group of the entire dataset.
In this study, we first performed a joint GPA of the entire dataset of original and mirrored configurations. Semi-landmarks were kept in their equidistant positions along the outlines [42]. Although differences between the GPA based on fixed equidistant positions of semi-landmarks and their sliding along the tangents to the outline curve usually only marginally influence the actual shape representation of the objects, introducing an additional step of semi-landmark sliding still slightly alters the actual morphology of the analysed shapes, which may introduce an artificial signal into the data [43]. This is because semi-landmark sliding along the tangents to the outline curve typically leads to different densities of the aligned semi-landmarks in different parts of the outline [44]. This then leads to increased shape differences between its original and mirrored configurations, thereby increasing the shape asymmetry of individual objects and the entire dataset. The Procrustes aligned configurations were subjected to PCA, and the marginal occupied position on PCs that illustrated the largest proportion of symmetric and asymmetric shape variation in frustules were plotted using the deformation grids of individual configurations from the mean shape of the entire dataset.
To obtain explicit values of the bilateral shape asymmetry of individual frustules, the original and mirrored configurations of each object were separately aligned by GPA ( Figure 2). The average Euclidean distances between corresponding landmarks from the original and mirrored configurations of each object were used to illustrate the relative differences in the shape asymmetry of different parts of the analysed frustules. The Procrustes distance (i.e., the square root of the sum of squared differences in the landmark positions of two shapes) between the original and reflected configurations was used to describe the shape asymmetry of individual objects. The same procedure involving a series of GPAs based on the original and mirrored configurations was also performed with the configurations solely spanning the dorsal and ventral valve sides. These analyses were used to identify different levels of asymmetry between these parts of the frustules.
Symmetry 2021, 13, x FOR PEER REVIEW 5 of In this study, we first performed a joint GPA of the entire dataset of original a mirrored configurations. Semi-landmarks were kept in their equidistant positions alo the outlines [42]. Although differences between the GPA based on fixed equidistant po tions of semi-landmarks and their sliding along the tangents to the outline curve usua only marginally influence the actual shape representation of the objects, introducing additional step of semi-landmark sliding still slightly alters the actual morphology of t analysed shapes, which may introduce an artificial signal into the data [43]. This is becau semi-landmark sliding along the tangents to the outline curve typically leads to differe densities of the aligned semi-landmarks in different parts of the outline [44]. This th leads to increased shape differences between its original and mirrored configuration thereby increasing the shape asymmetry of individual objects and the entire dataset. T Procrustes aligned configurations were subjected to PCA, and the marginal occupied p sition on PCs that illustrated the largest proportion of symmetric and asymmetric sha variation in frustules were plotted using the deformation grids of individual configur tions from the mean shape of the entire dataset.
To obtain explicit values of the bilateral shape asymmetry of individual frustules, t original and mirrored configurations of each object were separately aligned by GPA (F ure 2). The average Euclidean distances between corresponding landmarks from the or inal and mirrored configurations of each object were used to illustrate the relative diff ences in the shape asymmetry of different parts of the analysed frustules. The Procrust distance (i.e., the square root of the sum of squared differences in the landmark positio of two shapes) between the original and reflected configurations was used to describe t shape asymmetry of individual objects. The same procedure involving a series of GP based on the original and mirrored configurations was also performed with the config rations solely spanning the dorsal and ventral valve sides. These analyses were used identify different levels of asymmetry between these parts of the frustules. The frustule size was evaluated by the centroid size of the original landmark conf urations. Centroid size, defined as the square root of the sum of squared distances of ea landmark from the centroid of each configuration, provides an unbiased morphomet size measure of analysed objects that can be used for subsequent comparison with th shape characteristics [45]. In this study, the centroid size of each frustule was compar with their shape asymmetry evaluated by the Procrustes distance between the origin and mirrored configurations. To investigate the relationship between asymmetry a frustule size, four alternative functions were compared (Table S1). The linear function w based on the ordinary least squares regression. Then, the quadratic, exponential, a power law functions were also used to evaluate the relationship between the size a asymmetry measures of individual frustules. The Akaike information criterion (AIC) w used to compare these four models and determine which one might be the best fit for t observed data [46].
GPA was implemented using the function procGPA in the shapes package, ver. 1.2 [47], in R, ver. 4.0.5 [48]. PCA of the dataset consisting of the Procrustes aligned origin and mirrored configurations was conducted in TpsRelw, ver. 1.65 [29]. The deformati The frustule size was evaluated by the centroid size of the original landmark configurations. Centroid size, defined as the square root of the sum of squared distances of each landmark from the centroid of each configuration, provides an unbiased morphometric size measure of analysed objects that can be used for subsequent comparison with their shape characteristics [45]. In this study, the centroid size of each frustule was compared with their shape asymmetry evaluated by the Procrustes distance between the original and mirrored configurations. To investigate the relationship between asymmetry and frustule size, four alternative functions were compared (Table S1). The linear function was based on the ordinary least squares regression. Then, the quadratic, exponential, and power law functions were also used to evaluate the relationship between the size and asymmetry measures of individual frustules. The Akaike information criterion (AIC) was used to compare these four models and determine which one might be the best fit for the observed data [46].
GPA was implemented using the function procGPA in the shapes package, ver. 1.2.6 [47], in R, ver. 4.0.5 [48]. PCA of the dataset consisting of the Procrustes aligned original and mirrored configurations was conducted in TpsRelw, ver. 1.65 [29]. The deformation grids illustrating shapes typical for individual PCs were also obtained using this software. Analyses examining the relationship between size and shape asymmetry of frustules were performed in PAST, ver. 4.04 [49]. The R scripts used for the analyses are available online at https://doi.org/10.5281/zenodo.5732459 (accessed on 16 December 2021).

Results
Variation among individual frustules of E. bilunaris accounted for 92.3% of the total shape variability that was described by symmetric PCs yielded by PCA of the original and mirrored configurations of frustules. The remaining 7.7% of the variation was accounted for by asymmetric axes that described the bilateral asymmetry of configurations.
The single dominant feature of shape variation was a bilaterally symmetric pattern distinguishing cells with narrowly elongated shapes from those that were typically relatively shorter and wider with a central projection on the ventral side (Figure 3). This trend accounted for 80.5% of the total variation, as illustrated by PCA of the Procrustes aligned data of the original and mirrored configurations. Changes in the shape of the central part of the ventral outlines were also a prominent feature of the symmetric variation patterns described by PC2 and PC4, which accounted for 5.8% and 2.1% of the total shape variability, respectively. In these cases, however, the changes in the shape of this central part of the ventral outline were not associated with any notable morphological changes in other parts of the valves (Figure 3). Two main patterns of bilateral shape asymmetry were illustrated by PC3 and PC5, accounting for 2.9% and 1.5% of the variation, respectively. Interestingly, both patterns largely involved asymmetric deviations of the ventral sides ( Figure 3). grids illustrating shapes typical for individual PCs were also obtained using this software. Analyses examining the relationship between size and shape asymmetry of frustules were performed in PAST, ver. 4.04 [49]. The R scripts used for the analyses are available online at https://doi.org/10.5281/zenodo.5732459 (accessed on 16 December 2021).

Results
Variation among individual frustules of E. bilunaris accounted for 92.3% of the total shape variability that was described by symmetric PCs yielded by PCA of the original and mirrored configurations of frustules. The remaining 7.7% of the variation was accounted for by asymmetric axes that described the bilateral asymmetry of configurations.
The single dominant feature of shape variation was a bilaterally symmetric pattern distinguishing cells with narrowly elongated shapes from those that were typically relatively shorter and wider with a central projection on the ventral side (Figure 3). This trend accounted for 80.5% of the total variation, as illustrated by PCA of the Procrustes aligned data of the original and mirrored configurations. Changes in the shape of the central part of the ventral outlines were also a prominent feature of the symmetric variation patterns described by PC2 and PC4, which accounted for 5.8% and 2.1% of the total shape variability, respectively. In these cases, however, the changes in the shape of this central part of the ventral outline were not associated with any notable morphological changes in other parts of the valves (Figure 3). Two main patterns of bilateral shape asymmetry were illustrated by PC3 and PC5, accounting for 2.9% and 1.5% of the variation, respectively. Interestingly, both patterns largely involved asymmetric deviations of the ventral sides ( Figure 3).  five principal components (PC) that described most variation in the shape space of the original and mirrored configurations of 64 analysed valves. The shape dynamics associated with PC1, PC2, and PC4 involved both symmetric halves of the frustules equally. In contrast, PC3 and PC5 described the bilaterally asymmetric shape changes within the dataset.
When looking at the average asymmetric variability of the individual frustules, the ventral side of the valve outlines was clearly more asymmetric than the dorsal side (Figures 4 and 5). Moreover, most of the asymmetry of the frustules was due to shape deviations on their ventral sides. The highest levels of asymmetry were concentrated in the central parts of the ventral curve ( Figure 5), except for the central portion itself, where the shape variability was primarily manifested in the symmetric components.
When looking at the average asymmetric variability of the individual frustules, the ventral side of the valve outlines was clearly more asymmetric than the dorsal side (Figures 4 and 5). Moreover, most of the asymmetry of the frustules was due to shape deviations on their ventral sides. The highest levels of asymmetry were concentrated in the central parts of the ventral curve ( Figure 5), except for the central portion itself, where the shape variability was primarily manifested in the symmetric components.  Comparison of the different models for the analysis of the relationship between frustule size and asymmetry (Table S1) favoured the power law relationship with a negative scaling coefficient of −0.38 (95% C.I. (−0.32, −0.43)). This was corroborated by a strongly significant linear log-log relationship between these variables (r = −0.80, R 2 = 0.64, p = 0.0001), indicating that the relatively large frustules within the dataset were usually considerably less asymmetric than the smaller cells. Additionally, 20% of the smallest frustules in the sample accounted for approximately half the range of the observed bilateral asymmetry ( Figure 6). When looking at the average asymmetric variability of the individual frustules, the ventral side of the valve outlines was clearly more asymmetric than the dorsal side (Figures 4 and 5). Moreover, most of the asymmetry of the frustules was due to shape deviations on their ventral sides. The highest levels of asymmetry were concentrated in the central parts of the ventral curve ( Figure 5), except for the central portion itself, where the shape variability was primarily manifested in the symmetric components.  Comparison of the different models for the analysis of the relationship between frustule size and asymmetry (Table S1) favoured the power law relationship with a negative scaling coefficient of −0.38 (95% C.I. (−0.32, −0.43)). This was corroborated by a strongly significant linear log-log relationship between these variables (r = −0.80, R 2 = 0.64, p = 0.0001), indicating that the relatively large frustules within the dataset were usually considerably less asymmetric than the smaller cells. Additionally, 20% of the smallest frustules in the sample accounted for approximately half the range of the observed bilateral asymmetry ( Figure 6). Comparison of the different models for the analysis of the relationship between frustule size and asymmetry (Table S1) favoured the power law relationship with a negative scaling coefficient of −0.38 (95% C.I. (−0.32, −0.43)). This was corroborated by a strongly significant linear log-log relationship between these variables (r = −0.80, R 2 = 0.64, p = 0.0001), indicating that the relatively large frustules within the dataset were usually considerably less asymmetric than the smaller cells. Additionally, 20% of the smallest frustules in the sample accounted for approximately half the range of the observed bilateral asymmetry ( Figure 6).

Discussion
The analysed deformations in the shape of the valves of E. bilunaris were visually similar to those illustrated in the genus Eunotia in several previous studies [3,7,9,10,50,51]. These deformations primarily affected the ventral part of the frustules. In addition, our

Discussion
The analysed deformations in the shape of the valves of E. bilunaris were visually similar to those illustrated in the genus Eunotia in several previous studies [3,7,9,10,50,51]. These deformations primarily affected the ventral part of the frustules. In addition, our morphometric analysis also showed that the overall rate of these deformations increased significantly with decreasing cell size. In diatoms, this basically means that the cell outlines were more distorted as they underwent more vegetative divisions because the cell size decreased with each binary division throughout the asexual vegetative cycle [2,52]. However, not all shape deformations were asymmetric with respect to the bilaterally symmetric Bauplan of eunotioid diatoms. The most significant shape trend in our data, which was captured by PC1 as illustrating more than 80% of the total symmetric and asymmetric variation, was clearly related to the allometric shape-to-size relationship in which cells progressively shrink during a series of vegetative divisions while also changing shape symmetrically in both halves of the frustule. This phenomenon of shape allometry associated with a diminution series during the vegetative cycle has been illustrated in pennate diatoms by several examples [19,21,23,53,54]. In E. bilunaris, it mainly involves a gradual shortening of the cells, while largely maintaining their relatively stable width [52,55]. In terms of shape variation, this process also included the production of a central bump on the ventral side of the valves, as shown by the shape changes associated with PC1. However, these allometric shape changes were detected in E. bilunaris both in natural populations and in clonal cultures under optimal growth conditions and are therefore unlikely to be a manifestation of teratogenic morphogenesis [52,55]. Thus, the formation of a central bump in small cells is probably an intrinsic part of the allometric shape pattern in this species.
The asymmetric deviations from bilaterally symmetric shapes captured by our morphometric analysis also visually resembled those previously recorded in various species of the genus Eunotia in recent populations or in fossil material associated with various types of environmental stress, such as acidification or heavy metal contamination [3,7,9,10,13,56]. Our analysis showed that the regions of the ventral outline near the central part of the valves were the most prone to asymmetric deformations, while the average degree of asymmetry decreased significantly towards the apical cell poles. The allometric pattern of markedly increasing frustule asymmetry with decreasing size, observed in E. bilunaris, represents a distinctly different phenomenon from the more or less stable rate of subtle asymmetric variation and individual fluctuating asymmetry of frustules during the asexual diminution series that were previously illustrated in various species of the genera Navicula, Luticola, and Sellaphora [19,22,57]. These results confirm that when a teratogenic deformity arises in a population, it typically persists and becomes increasingly amplified through the asexual reproduction series [9,58]. The strong allometric relationship of the asymmetric deformations in E. bilunaris, contrasting with the absence of similar patterns in the subtle shape fluctuation asymmetry in the aforementioned genera grown under optimal conditions, provides indirect evidence that the observed variation represented a teratogenic response at the population level that was accentuated in the course of the vegetation cycle.
Several previous studies and reviews have stressed that assessing the severity of frustule teratology is an important step for further advancing research on the relationship between these morphogenetic deformities and their causative environmental factors [4,5,58,59]. We suggest that the presented geometric morphometric framework for the analysis of bilateral asymmetry across the transapical axis of eunotioid diatoms represents one such approach for quantifying shape teratogeny in this diatom group. Previous geometric morphometric studies that addressed teratogenic deformations of diatom frustules were mostly based on the analysis of overall shape variability in the valves, i.e., without explicitly distinguishing symmetric and asymmetric subspaces of the variation [11,[14][15][16]. This obviously resulted in a situation when most of the detected patterns (described by the first few PCs) highlighted symmetric differences among individual frustules. Importantly, the specific position of frustules in relation to the average shape also necessarily depended on the structure of the particular dataset under study [60].
In contrast, the geometric morphometric analysis of the asymmetry of individual frustules allowed us to quantify their individual degree of deviation from the ideal bilaterally symmetric pattern regardless of other analysed cells. With the same number of equidistant points registering the valve outlines, the results could be compared among populations, studies, or even different taxa with eunotioid frustule outlines. When investigating diatoms with different types of valvar symmetry, the design of the asymmetry analysis would have to be adjusted to their respective symmetry groups. For example, three perpendicular geometric components of shape asymmetry can be decomposed in configurations with biradial symmetry that are typical for most araphid diatoms and several raphid genera, such as Frustulia, Luticola, Navicula, and Achnanthidium [19,20,23,57]. In addition, teratogenic shape asymmetry in genera with heteropolar valves, such as Gomphonema, Licmophora, and Didymosphenia [61], needs to be quantified based on the reflection of configurations across the apical axis. However, in all cases, the Procrustes-based decomposition of the shape asymmetry could represent a well-applicable framework for scoring the severity of diatom frustule teratology in relation to different types of environmental stress, such as the metal contamination, acidification, eutrophication, or increased radioactivity [5].
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/sym14010042/s1, Table S1: Comparison of the different models for the analysis of the relationships between frustule size and shape asymmetry of E. bilunaris.
Author Contributions: K.W. and J.N. conceived and designed the research; K.W. took the microphotographs and digitised the data; J.N. conducted the geometric morphometric analysis; K.W. and J.N. wrote the paper. All authors have read and agreed to the published version of the manuscript.