Morphometric Analyses of ( Hidden ) Directional Asymmetry in Leaf Blades

The deviation from perfect bilateral symmetry is a phenomenon actively exploring in evolutionary and environmental studies. The bilateral variation presents on different ecosystem level. The methods applied vary depending on the task and the final goal of study. The present study demonstrates the statistically significant presence of components (traces) of directional asymmetry (DA), fluctuating asymmetry (FA) and antisymmetry (AS) in leaf blade of Betula pendula on low level of ecosystem such as tree and leaf blade. The Generalized Procrustes Analysis was applied for testing asymmetry of shape in all data set and subgroups. The category ‘population’ revealed no significant value of DA (factor ‘side’). In the levels ‘tree’ and ‘leaf blade’ the factor ‘side’ was statistically significant as well as fluctuating asymmetry (factor interaction ‘leaf × side’).The principal component analysis showed visually the difference in PC scores between antisymmetry matrices of the left and right halves of the leaf blade. Covariate analysis matrices demonstrated the shape deviation from strict symmetry. The metric traits showed directional asymmetry in t-test in the leaf blade and the low kurtosis values. Permutation test of kurtosis values in geometric morphometric approach showed deviation from normality that could verified as a weak presence AS traces. The discriminant analysis results showed that traces of DA differed at the subpopulation level, as three of 10 populations revealed not significant factor ‘side’. Descriptive statistics of the metric train showed a correspondence to the trace of directional asymmetry in the shape of leaf blade. Fluctuating asymmetry in its pure form, at three levels of ecosystem was met only in single population of ten that should be taken into account testing developmental stability of birch and possibly other woody plants.

The morphogeometric method developed for object asymmetry was used here to take into account the fact that the left and right halves of leaves are connected mirror images of each other while the axis of symmetry passes through the entire structure [8], [10,11].Antisymmetry is a third kind of bilateral asymmetry, characterized by low kurtosis values of difference between left and right side values.The bimodal distribution of (R -L) value is typical for antisymmetry.This kind of asymmetry is widely spread in many taxa [12].There are many studies reported antisymmetry in plant.Nevertheless '… some studies have found that antisymmetry can be fairly difficult to demonstrate with geometric morphometrics' [13].Among plants some studies demonstrated the successful clusterization of leaves with antisymmetry [14].Subtle antisymmetry might be fairly widespread in plants, because leaf asymmetry relates to leaf phyllotaxis, the arrangement of leaves along the shoot.
The normalizing difference is a trivial approach for FA testing.An alternative method is considered to be the geometric morphometric methods (GMMs) which takes into account the labels placed on the bilaterally symmetric structures [9][10][11].The labels are rejected from the centroid points of consensus figure, which is drawn by the averaging of landmarks in Cartesian coordinates and the value of the FA shape of an organ (or the whole organism) is evaluated in a mixed-model two ANOVA.
According to Leamy et al, 2002Leamy et al, , 2005 and other many sources, the fluctuating asymmetry is encoded by several epistatic genes [13][14][15][16][17].In plants DA is tested approximately in 10% of some dimensional traits, for example in the leaf blades of woody plants, as silver birch [18].In English oak with heterogeneous structure of leaf blade the amount of such traits was higher and DA tested in GMMs revealed in all populations observed [19].
Phenogenetics can be characterised as a study of homologous phenes and their individual combinations in evolutionary-ecological aspects.Currently, population ecology (in the field of phenogenetic monitoring) is actively developing by employing phenotypic traits [20,21].Most of the work on phenogenetics employed the geometric morphometric methods in the study on populations of rodents and insects [26,27].Therefore, a detailed study of phenogenetics according to the program from the phenotype to the genotype in the kingdom Vegetabilia is expected.
Directional asymmetry is in focus of many researchers.Some work report the 2 way mixed-model ANOVA gives an undesirable bias in the value of the FA, for example depending on character of DA in the mixture FA and DA in the sample [6].Therefore, the phenotypic effect of fluctuating asymmetry is difficult to test [1,2,13], although some approaches are used to portion the value of FA and DA [6,7,22].
In previous studies the correlation obtained between the magnitude of the FA index found by linear measurements and the value of the FA index, produced by the method of geometric morphometrics.Such a correlation cannot be regarded as mandatory, and depends on the magnitude of linear traits, making a greater contribution to the shape of bilaterally symmetrical halves.The magnitude of the FA and the stability/instability of development depended on some factors including the value, the kind of pollution, and the height of the relief [19].
The study of the relation between genotypic and phenotypic effect in populations of the plants was carried out depending on the habitat and climatic peculiarities, with the use of traditional linear methods for FA and DA testing [22,23].The apparent simplicity of the methodology very often led to a distortion of the results because of difference in size and heterogeneity of variables caused by environmental factors [13,25].Recent studies of fluctuating asymmetry in birch leaves indicated the presence of a paradoxical non-monotonic effect and hormesis in which a relatively small toxic dose increases the asymmetry value [28].
The studies in geometric morphometrics in developmental stability of plants possess many aspects.They are morphological integration [29], phenotypic plasticity [30], shifts in shape symmetry and role of genes influencing symmetry and asymmetry [31], and differentiation in shape in context of evolution and taxonomy [8,14,29,32].Some studies indicate that the size and the shape of the edge of the leaf can play a significant role in the reproducibility of the results of the magnitude of the fluctuating asymmetry [29,34].The accuracy of measuring, visual inspection of laminas undoubtedly is also essential.The unbiased estimation of asymmetry, both on the population and individual level also is in focus [35].Nevertheless, the FA index as an indicator of stability of development remains a tool showing the level of environmental impact [36][37].
The purpose of this study was to examine the magnitude of the two kinds of asymmetry (FA and DA) in shape of a Betula pendula Roth leaf in relatively normal ambient conditions (no serious pollution sources) at different levels of ecosystem (population, tree, leaf blade).Conventional method of testing on the base of metric traits in comparison with geometric morphometric methods was applied.A model of the covariation matrix of the Procrustean distances in the shape of the leaf plate was used.The working hypothesis was the following: the increase in the accuracy of the measurement makes it possible to determine hidden structure of asymmetry by the method of geometric morphometrics.Under the structure the character of distribution asymmetric homologous attributes was understood.

2-1-Sites for Collection
The silver birch (warty birch, or East Asian white birch) has a very wide area in Eurasia; in Russia its areal extends from the tundra to the steppes.Vladimirskay oblast, European part of Russia (56°5′0″ N; 40°37′0″ E) was chosen for study.The sites for leaf collection included as very close (2-5 km in limit of cities) as well remote ones up to 70-90 km.So "populations" considered as conditional ones and pronounced better as cenopopulations, as they included other species and forms of plants.The collection of leaf blades (laminas) was carried out in 2016, September.50 leaf blades from each population of 10 trees were selected using trees of the same age and the same generative stage of development.In a whole the leaves were sampled from trees of age 15-25 years according to the method developed by Kryazheva et al. [36] and adopted by Gelashvili et al. for populations under different environmental stress level [38].The trees were located on distance 2-5 m of each other.From each tree 5 laminas were taken from the shortened shoots randomly on the height 1.5-2 m under conditions of relatively the same sun lightening.To reduce the possible allometric measurement error, leaves with a maximum leaf half-width 3-3.5 cm were selected.The sites for collecting differed in height of the relief, and therefore differed in the physical and chemical properties of the soil.The leaves were gathered evenly around the lower part of the tree crown and storage dry under the paper press for two weeks.

2-2-Imaging
The images of leaves were taken using a Panasonic DMC-FZ100 camera (resolution 14 mega pixels), and JPG file format was used.Each lamina was photographed separately, the grid lines (option 'flower') were used as guides in order to obtain the maximal vertical image.The midvein (the costa) was aligned in parallel with the grid line.The back, not glossy leaf side was used.Each plate was photographed twice to calculate measurement error.The paired bilateral vessel is a traditionally used trait in study of FA B. pendula [18,20,34,36].Nevertheless the most defined vessels being chosen and 12 landmarks were labelled in the same order on each picture, after setting a scale factor.Files for data manipulation were created using the TPS software package [39,40].The landmarks were digitised twice and were classified as homological landmarks type I as represented by pair labels on the endpoints of the lateral vessels (Figure 1).The method of Procrustes ANOVA as analogue of two-way mixed-model ANOVA (individual × side) which is used for FA value test in metric and in meristic traits was employed.Procrustes fit presents a space in limit of centroid size.Accordingly this method, the right and the left point were aligned along with the mirror-reflected landmarks.The Procrustes alignment included the original and the mirrored configurations of a sample combined, and superimposed all of them simultaneously.For averaging consensus formation the method of least squares was used.In detail this method is observed in many basic studies, for example [9,39].
In present study the MorphoJ1.06dpackage was used, available on web site www.morphometrics.org[41].The total TPS file for population sample presented a single file with all data (50 leaf blades × 4 times repeats).The sampling procedure resulted in a dataset, with leaves nested within trees and trees nested within populations.Thus the whole plan consisted of table 2000 × 6, first column included coordinates and served as identificator, others columns included coded values for classifiers: population, tree, leaf, image and digitizing.After creating a Procrustes fit the generalized Procrustes analysis was conducted for all data set and for each ecosystem level.

2-3-Evaluation of Fluctuating and Directional Asymmetry and Antisymmetry
Using metric traits the FA was obtained on formula normalizing difference FA= |R -L|/ (R + L), where R and L is respectively the right and the left trait value.Conventional 2 way mixed-model ANOVA was used as well.
A Generalized Procrustes Analysis was applied to all landmark configurations to extract shape by removing effects of scale, orientation and rotation (references to be included for the Procrustes fit) [5,8,10,11,42].The value of FA is determined via a two-way mixed-model analysis of variance from the magnitude of the mean square variance components of two factors interaction: 'sample" (random) and 'side" (fixed).The first factor is denoted by the code values corresponding to the level of population variability, to the individual (tree) or to the organ (leaf blade).Factor 'side" is denoted only by two code values ('right' and 'left').The value of the mean square of the factor 'side' indicated the presence of directional asymmetry.
The magnitude of the fluctuating asymmetry was determined by the mean square MS and the value of the F-Goodall criterion evaluating the interaction of one of the random factor 'population', 'tree', or 'leaf' by fixed factor 'side'.
The directional asymmetry of the linear characteristics was verified by the t-test with the verification of the null hypothesis H0, on the equality of the right and left trait value.Also the two-way ANOVA was performed with next treatment data in Excel with mixed-model significance level obtaining [1].
The statistical significance of DA by method of geometric morphometrics was tested on the value of MS of factor 'side' and on the value of F-Goodall criterion.The difference in eigenvalue variance between asymmetry matrix for left and for right side was also observed.
To test the antisymmetry the value of kurtosis was observed from descriptive statistics and t-test was used [1].The normality of PC scores for each score was checked by Shapiro-Wilks test.The permutation multiplication of samples normalizing their distribution was tested with the value, statistics and p of normality for kurtosis (Mardia tests).

2-4-Visualization
For visual acceptance the covariance matrix of asymmetry was generated.The asymmetry covariance matrix reflected the variability of the bilateral structures with the testing of the null hypothesis that there is no difference in the spatial position of the bilateral paired landmarks.
The asymmetry covariance matrix was generated from left and right landmarks separately and showed the picture of distribution principal components for left and for right sides in the comparative aspect.Partial least squares analysis (PCA) was used for analysis of a scatter plot between left and right side asymmetry matrices.Multivariate regression of asymmetry components left and right landmarks was provided for testing dependency: centroid size -shape.Canonical variate analysis (CVA) is well optimized to find difference among samples therefore.CVA was used for visual acceptance of difference between overall shapes of subgroup asymmetry component.Discriminant function was applied for testing distinguishing between subgroup levels.

3-1-Measurement Error
The leaf blade represented a true replication, because each plate was measured twice.The author takes into account the point of view Klingenberg [13] about two types of errors: digitization (positioning) and measurement (labelling).Therefore, a repeated survey of each plate was carried out, with a double marking on each image.The measurement error was the sum of the errors in photographing and in marking.The imaging and digitizing were performed as separate effects in the model with one nested into the other [42] (Table1): The error of imaging and digitizing were less in leaf blade level in comparison to the level of population and level of tree.

3-2-Visualizes Shape Patterns
Left and right covariance matrix generated, consisted the data were pooled from population and from leaf (for comparison).Matrices in this data included right and left centroid size and Procrustes coordinates.Regression between left and right side centroids size (only option CS in regression analysis) showed several outliers and they were removed (option 'find outliers in data set').The principal component analysis for both covariance matrices was fulfilled separately.Six principal components were generated for each sample.After outliers swap the Procrustes fit was performed again and Procrustes ANOVA was repeated for each population.The principal component analysis showed different picture of distribution PC scores for left (A) and right (B) covariate antisymmetry matrices, the Procrustes coordinates were pooled by leaf (Figure 2).Partial least squares analysis of asymmetry component of covariation revealed the coefficient overall strength of association between blocks RV = 0.57.The outliers were notable on the graph.The procedure of outliers deleting was provided.From 3 up to 5 landmarks were swept in each population data set.This procedure controlled manually in Box plot graphs (in PAST).The regression indirectly confirmed the absence of serious presence of antisymmetry.

3-3-Effect of Size on Shape
Effect of size on shape was tested visually.To check a possible deviation in shape associated with size difference a multivariate regression was provided.The shape assumed as an asymmetry matrix pulled by leaf and played role of dependent variable.Centroid size (independent variable) was transformed to its natural logarithm to increase the fit of the model (Figure 4).On the picture a cloud line was perpendicular to the size OX axis to zero point.This meant the shape was variable but did not depend on size linearly.12 landmarks were removed and Procrustes ANOVA was performed.

3-4-Generalized Procrustes ANOVA
Procrustes ANOVA was fulfilled in series for individual as a population, as a tree, and as a leaf (3 set of analysis).
Results are placed in one Table 2.The mean squares for FA and individual variation exceeded the error component by more than 10 fold, indicating that measurement error was negligible relative to the biological shape variation.The magnitude of the F-Goodall criterion of the source 'side' increased noticeably from the population to the leaf.
Association between centroid size and asymmetry component matrix (leaf level) was not detected p = 0.113).For population level a weak difference was detected (p = 0.03).This meant the weak genetic difference in populations on asymmetry component of shape dependence from size of leaf.

3-5-Testing of Kurtosis and Skewness
In a whole the observed samples of matrices grouping data by leaf included information on traces of three types of symmetry.The distribution of row data R and L (where R and L are value of the right and respectively the left metric trait) were not normally.The linear traits were checked on extreme value of kurtosis (R -L) in t-test.The highest and lowest values of kurtosis were tested by sequential Bonferroni correction.The values were not obtained as significant.Even lowest value -1.65 (population №4, 3 d trait) did not show significant antisymmetry presence.The tabulated data of critical values of the kurtosis [1] test statistic for deviations of frequency distributions showed a serious presence of platykurtosis (with broad-peaked or bimodal distribution).Table №3 displays the results of testing metric traits, including t test for kurtosis value and a results FA and DA testing in Procrustes ANOVA for individual as a population including mixed-model significance level of DA through F ratio, as MSS/MSIS [1].Normalised index FA2 showed the highest value in population №10.This sample characterized by rather high value of factor 'side' in tree level, highest ratio MSS/MSI×S.The highest DA in the population corresponded to the high DA in the trees, for example: F = 7.52 for 'tree'; F = 15.09 for 'leaf' (population №9).Thus, with a change in the level of the biosystem in the direction from the population to a lower level, there was an increase in the magnitude of the directional asymmetry (df = 10 in Procrustes ANOVA).FA2 was crucially biased by presence DA in population level and trace of DA in level leaf.
The values of kurtosis in Mardia test were on the same level 104.9-106.2 in all samples.The values of skewness were also similar 42.1-45.3.Shapiro-Wilks and other similar tests showed normal or close to normal distribution of PC score for each Principal component (p >0.05).Skewness and kurtosis were close in samples (skewness is 15.03-15.05,kurtosis is 31.51).So low value p suggested the high kurtosis displayed AS presence in morphometric data of Procrustes coordinates.
The platykurtic distribution (negative value of kurtosis) pointed out on weak presence of antisymmetry on the background of FA and DA that meant the mixture of three types of asymmetry.Even in population with clear FA the histogram pointed out the platykurtosis (broad-peaked, bimodal, or many modal).Accordingly Palmer and Strobeck [1], the leptokurtic distribution is a result of mixture platykurtosis and normal distribution.If so the common leptokurtic distribution included the element or trace of antisymmetry.
Only the first population, showed a 'clear' fluctuating asymmetry of the leaf blades.In this population the mixedmodel ANOVA showed a value FA10 = 0.004 as a non-directional asymmetry variance after removing measurement error, which possessed 23.9% of total between side variance.It corresponded to FA2 = 0,103 for metric trait №3, the largest FA value trait.
Four populations from ten revealed significant DA (№№ 7-10).Population №7, №9 had not DA in linear trait, DA (traces of DA) presented in leaf shape only.The classical method of normalizing difference showed a statistically significant presence of directional asymmetry (t-test, H0: R = L, where R and L are linear trait value) in population №8 and №10, last exposed DA in two traits.Thus, these metric traits, including the angle between the middle and second veins (in 10 th population), demonstrated a relationship between the linear characteristics and the directional asymmetry of shape.
MANOVA was used for determine significance of DA ('side') of asymmetry component shape variation (parametric) in covariate matrix.The Pillai' trace criterion of asymmetry component shape variation was significant on both level 'tree' and 'leaf', correspondingly, p = 0.2 and 0.09 (both p<0.0001).The level "population" showed the insignificance as well as in hierarchal Procrustes ANOVA (Pillai' trace 0.46; p = 0.66).
Beside of Generalized Procrustes ANOVA, a modified Procrustes ANOVA was fulfilled.The leaf was as individual random effect, tree and population played a role of the extra main effect.The results was close to one performed in GPA, degree of freedom was higher as expected, p-level was the same (<0.0001)for individual variation of leaf category, for factor side and for factors interaction.For tree as an individual and population as an extra main effect results was also significant with the same p-level (p < 0.0001).

3-6-Clustering Population
In order to get clear the difference between matrices of asymmetry the population with lowest and highest value 'ind× side' were compared.

4-Discussion
This article represents the first contribution to the study of asymmetry in subpopulation level in context of relationship between asymmetry in shape and in metric traits.As a species Betula pendula showed a high phenotypic plasticity from absence of DA (clear FA popul.№1) to highly statistically significant level of DA traced in leaf biosystem level.The results showed that the presence of traces directional asymmetry in the leaf level shape is a common type of asymmetry in the birch leaf blades.
A similar study was conducted in 2015 without repeats of leaf blade images.Procrustes ANOVA showed the presence of DA in populations with a level of statistical significance p = 0.007 (F-Goodall = 3.19).Accordingly, the 'tree' showed a higher statistical significance (p < 0.001; F-Goodall = 4.97.The leaf blade showed highest significance (F-Goodall = 12.17; p <0.001), that confirmed the presence of directional asymmetry trace in shape of leaf blades a year earlier, under relatively similar environmental conditions.Only one population from seven showed 'clear' fluctuating asymmetry [19].In the current study, after increasing df, as the number of measurements was increased and the measurement error reduced, DA was not obtained at the population level and was not statistically significant.
The study showed that the variability of linear characteristics effect the asymmetry of the shape leaf blade in plants.Variation in FA of linear traits significantly influences the shape of the leaf blades that was typically for Ouercus robur.The method of geometric morphometric could be seemed as preferable for evaluating fluctuating asymmetry and stability of development because method takes into account the shape of the photosynthesis organ, i.e. lamina.The method of the normalizing difference takes into account only the sum of the values of the FA individual bilateral traits anyone of which cannot be free of directional asymmetry or antisymmetry.If populations are genetically differing we can say they carry individual genetic component highly heterogeneous for Betula pendula.
There is a controversial means on technogenic stress factors influencing developmental stability of birch species.Some factors (and their interaction) may decrease or increase the magnitude of the FA [37,43,44].At least in not stressful ambient the FA of shape highly differed among populations (up to 20 times).The presented study confirmed the heterogeneity of DA in populations of the Betula pendula species.As follows from the reports [22,45], high heterogeneity value of DA contribute to the shift of the FA value on individual level in a geometric morphometric context.It is of interest to study with an increase in the number of repetitions of measurements and a greater decrease in the measurement error.It would be allow testing individual (tree) and intra individual (leaf) variability in shape asymmetry including variability in genotype component of DA.
Betula pendula is wind-blown pan mixed species.It could not be in interest for discrimination of geographical differences.Sooner B.pendula is interest of environmental stress assess.The traces of DA and traces of AS both presented in tree and leaf level.The last was confirmed by platykurtic type of kurtosis evenly presented in asymmetry of shape.The outliers generated a heterogenic structure of asymmetry.Thus, a working hypothesis is confirmed about the joint presence of both components of variability: genotypic and phenotypic in the asymmetry of the shape of the leaf blade.Fluctuating asymmetry in its pure form was met only in 10% population studied, that should be taken into account in assessing the stability of development of the birch and possibly other woody plants.

5-Acknowledgments
The author thanks VlSU employees and students for their assistance in collecting herbarium material and statistical processing of primary data.

Figure 1 .
A: 5 pair of landmarks and metric traits: 1-half width of the leaf, 2-the length of the veins of the second order, 3-the distance between the bases of the veins, 4-the distance between the ends of the veins, 5-the angle between the midrib and lateral veins; B: Procrustes fit of left landmarks №№1 to 5, and right landmarks №№8-12; C: Procrustes fit of all amount landmarks №№1-12, asymmetry component.Black dots show the landmarks after alignment by superimposition.The matrix of asymmetry (C) reflects the variation of landmarks on the left and the right sides in overall shape.

Figure 2 .
Figure 2. PC scores for left (A) and right (B) covariate matrices antisymmetry.OXprincipal component; OY -% variance First principal component of left covariate matrix included less (27.9%) overall variance, right first component included (35.1%).This affirmed the difference in variances left and right data.The asymmetry covariance matrix was checked on variation in asymmetry component levels: population, tree, and leaf.The analysis of a scatter plot of differences between left and right sides demonstrated that the analyzed samples (level of leaf) display the pair wise correlation of PLS score between left (№№2, 3, 4, 5) and right (№№7, 8, 9, 10) homologues landmarks was r = 0.38 (p<0.0001).PLS1 possessed 60.3% of total covariation score (Figure3).

Figure 3 .
Figure 3. Regression between two blocks (block 1 is left and block 2 is right landmarks)

Figure 4 .
Figure 4. Regression asymmetry component (asymmetry matrix pooled by leaf) from centroid size

Table 1 .
Measurement error for leaf blade Note: MSmean squares; SSsum squares; dfdegree of freedom; F -F-Goodall criterion, sums of squares and mean squares are in units of squared Procrustes distance.Measurement error 1 (positioning error) explained 10.7% of sum square individual variation of Procrustes distances.Measurement error 2 (error of digitizing) explained 8.9% of sum square individual variation of Procrustes distances.

Table 3 .
Results linear measurement and morphometrics

Table 4 .
Procrustes distances and T square for three populations (№1, 8 and 10) Discrimination in geographical context did not show the strong clustering according to spatial isolation.