DOI: 10.5061/dryad.pg4f4qrtr
title: Globularia bisnagarica L. (Plantaginaceae) shows genetic uniformity throughout its disjunctive range

description: The paper presents the results of studying the morphological and genetic
 variation in populations of two species from genus Globularia, Globularia
 bisnagarica L. and Globularia trichosantha Fisch. and May. Both species
 are relicts and have disjunctive ranges. The analysis involved 25
 populations of G. bisnagarica and 4 populations of G. trichosantha growing
 far apart from the main parts of their ranges. Sequencing of some regions
 of cpDNA (trnL−trnF, rbcL, and rpS16) and nuclear rDNA (ITS barcoding)
 revealed an extremely low level of genetic variability in the populations
 of both studied species. At the same time, G. bisnagarica demonstrates
 genetic uniformity within the entire disjunctive range. A study of the
 leaf blade shape variation with the geometric morphometrics approach
 revealed no clear interpopulation differences for both species; in most
 cases, the shape varied widely at the intrapopulation level. However, in
 the case of G. bisnagarica, a subtle trend was found: in populations from
 habitats with southern locations, the leaves had a narrow, oblong shape,
 while the populations located to the north, were characterized by wider,
 spatulate shape of the leaf blade. The assessment of variation in size
 parameters of plant organs using classical linear morphometrics revealed a
 more noticeable separation of populations for both G. bisnagarica and G.
 trichosantha. Modeling with environmental factors, designed to explain the
 observed patterns of variation, showed that the linear sizes of plant
 organs were more affected by climatic and topographic factors than the
 leaf blade shape. The main drivers of morphological variation in
 morphometric parameters of G. bisnagarica populations were the elevation,
 mean diurnal range of temperature, and annual precipitation. Variation in
 morphometric parameters of G. trichosantha populations mostly were driven
 by the annual mean temperature. Morphological variation observed in
 populations of G. bisnagarica within the studied fragments of the
 disjunctive range has exclusively modified nature caused by natural
 climatic factors and local conditions of population growth.
description: <strong>Morphometric
 analysis</strong> The material was collected in 2021–2022 from
 twenty-five G<em>.</em> <em>bisnagarica</em> and
 four <em>G</em>. <em>trichosantha</em> populations
 (Table 1, Fig. 1). Previously, the spatial structure of some of these
 populations was studied (Kondratieva et al., 2021).
 The studied populations were
 separated into several groups depending on the type and geographical
 location of the topographic features, to which their habitats are
 confined. Six groups were identified for <em>G</em>.
 <em>bisnagarica</em>: Volga Upland (VU), Sokskiye Yary (SYA),
 Bugulma-Belebey Upland (BBU), Obshchy Syrt (OS), Zilair Plateau (ZP), and
 Stavropol Upland (SU). <em>G</em>.
 <em>trichosantha</em> populations were divided into 2 groups:
 Greater Caucasus (GC) and Crimean Mountains (CM).
 For the linear morphometrics
 approach, the sixteen quantitative morphometric parameters of vegetative
 and generative organs were measured for 17–30 individuals in each
 population. These parameters included: number of rosettes, diameter of
 rosettes, number of leaves per rosette, length of rosette leaf, width of
 rosette leaf, thickness of rosette leaf, number of generative shoots per
 individual, number of generative shoots per rosette, length of generative
 shoot, diameter of generative shoot, number of leaves on generative shoot,
 length of the leaf on generative shoot, width of the leaf on generative
 shoot, thickness of the leaf on generative shoot, diameter of
 inflorescence, and height of inflorescence. 
 From the same plants, one
 rosette leaf per individual was cut off and herbarized for cameral treatment
 and further analysis by geometric morphometrics methods. Each herbarium
 specimen was photographed with a Nikon D3100 camera with the same focal
 distance against a white background with a ruler. Three fixed landmarks
 and two outlines were placed on the photographs, using the tpsUtil and
 tpsDig2 programs (Rohlf, 2015). Landmark No. 2 represented the leaf blade
 apex; No. 1 and 3 were the points where the leaf blade
 starts to expand. Each outline (a curve, describing the form of left and
 right sides of the leaf blade) consisted of 50 equidistant semilandmarks,
 limited by fixed landmarks at the ends (Supporting information). The total
 sample size amounted to 724 digitized images for <em>G</em>.
 <em>bisnagarica</em> and 120 images for
 <em>G</em>. <em>trichosantha</em>. The Euclidean
 coordinates of landmarks and semilandmarks were used as input data for the
 subsequent shape analysis. In order to study morphological variability in
 more detail, we conducted two complete independent morphometric analyses
 separately for each species. We applied the Generalized Procrustes Analysis
 (GPA) to extract shape data by removing the extraneous information
 about size, location and
 orientation (Zelditch et al., 2012). Because of the presence of
 semilandmarks GPA involved an additional step of iterative sliding of
 these points along the outline tangents to achieve their final position
 yielding the smoothest possible deformation of the actual configuration
 from the mean shape of the analyzed data set. The positions of these
 points were defined by the lowest possible bending energy between each
 configuration and the mean shape (Neustupa and Woodard,
 2021). Since the leaves of
 <em>Globularia</em> are bilaterally symmetrical objects, the
 configurations of landmarks that describe their shape contain not only the
 information about the differences between individuals but also about
 individual asymmetry (differences between the right and left side of the
 leaf). The nature of such data structure has to be taken into account in
 further morphometric analysis (Klingenberg et al., 2002). GPA considering
 the bilateral symmetry was performed using the
 <em>bilat</em>.<em>simmetry</em> and
 <em>g</em>.<em>pagen</em> functions from the
 ‘geomorph’ package (Adams et al., 2021).
 The intrapopulation
 morphological variability of the leaf blade shape was calculated using the
 <em>morphol</em>.<em>disparity</em> function from
 the ‘geomorph’ package as Procrustes variance, which was the sum of the
 diagonal elements of the group covariance matrix divided by the number of
 observations in the group (Zelditch et al., 2012). Additionally, the
 contribution of each population to the overall morphological variation
 (morphological partial disparities) was assessed.
 Principal component analysis
 (PCA) was used to identify the main patterns of morphological variation in
 the leaf blade shape in <em>Globularia</em>
 populations. Next, we
 studied the relationships between leaf blade shape and its size. The
 centroid size (CS) – the square root of the sum of the squared distances
 from the center of the leaf blade to each of the landmarks – was used to
 characterize the size of the leaf blade (Vasil'ev et al., 2018). The
 possible dependence of the leaf blade shape changes on its size
 (allometry) was estimated by the Procrustes ANOVA, using the
 <em>procD</em>.<em>lm</em> function performed in
 ‘geomorph’ package, with the matrix of Procrustes distances between the
 samples as a dependent variable, and the logarithm of the centroid size
 (log(CS)) as predictor. The same method was used to determine the
 significance of interpopulation differences in the leaf blade
 shape.  To investigate the relationships between the
 leaf blade shape and environmental variables, multiple linear mixed-effect
 regression models were used. The scores of the first two principal
 components (PCs), reflecting the main patterns of morphological variation,
 were used as dependent variables in the modeling procedure (Innangi et
 al., 2020). The fixed effects included 19 bioclimatic variables derived
 from temperature and precipitation, amount of solar radiation, water vapor
 pressure, elevation and slope angle (Supporting information). Climate and
 elevation data were obtained from the WorldClim database with 30 s (1 km)
 spatial resolution (Fick and Hijmans, 2017). Each variable was power
 transformed with optimal lambda using the
 <em>powerTransform</em> and <em>bcPower</em>
 functions from the ‘car’ package (Fox and Weisberg, 2019) for normality of
 distribution. All variables were standardized to zero mean and unit
 variance. The grouping factor “population” was considered as a random
 intercept. To avoid misspecification of predictors of real
 interest due to their high collinearity, we used climate and topographic
 variables to construct a correlation matrix (based on Spearman's
 correlation coefficients), further excluding one variable from the pairs
 for which |r| ≥ 0.7. The
 models were built using the <em>lmer</em> function from the
 ‘lme4’ package (Bates et al., 2015). The model parameters were calculated
 by REML estimation. We calculated conditional and marginal
 R<sup>2</sup> (Nakagawa, 2017) to evaluate the model fit using
 the package ‘performance’ (Lüdecke et al., 2021). Visualization of models output and checking the
 main assumptions were presented with the help of ‘sjPlot’ package
 (Lüdecke, 2022). The preprocessing of linear morphometrics data
 included the calculation of Spearman correlation coefficients to identify
 strongly correlated morphometric parameters (|r| ≥ 0.9). Since the analysis of the correlation matrices
 did not reveal strongly intercorrelated variables, all of them were
 retained for subsequent analysis. Interquartile range (IQR) was used as a measure
 of intrapopulation variability of morphometric parameters. In addition, a
 non-parametric Fligner-Killeen test for homogeneity of group variances
 based on ranks was performed (Conover et al., 1981). The procedures of
 identifying the main patterns of morphological variation and its
 relationships with environmental factors were carried out similarly to
 those for geometric morphometrics data.
 We set a significance level of
 0.05 for all computations. We used R version 4.1.2 (R Core Team, 2021) for
 all statistical analyses and data manipulations.
 <strong>Molecular
 genetics analysis</strong> DNA for molecular genetic studies was extracted
 from leaves dried in silica gel using the NucleoSpin® Plant II kit
 (MACHERY-NAGEL, Germany) according to the manufacturer's
 protocol. To analyze the entire sample (27 populations,
 three samples from each population), an intergenic transcribed ribosomal
 DNA spacer (ITS1-5.8S-ITS2) was used, which was amplified using primers
 ITS-1 and ITS-4 (White et al., 1990). The
 <em>trn</em>L−<em>trn</em>F intron, the
 <em>rbc</em>L gene, and the <em>rpS</em>16 gene
 were used as markers for the analysis of cpDNA. The choice of markers was
 due to the availability of sequences of these regions in GenBank for
 objects of our interest (<em>G</em>.
 <em>bisnagarica</em> and <em>G</em>.
 <em>trichosantha</em>) from other parts of the species'
 areas (Germany, Turkey, France, and Croatia).
 Polymerase chain reaction
 (PCR) was carried out in a volume of 50 µl. The reaction mixture contained
 10 µl of the prepared MaGMix reaction mixture (200 µM of each dNTP, 1.5 mM
 MgCl2, 1.5 units of SmarTaqDNA polymerase and buffer; Dialat Ltd., Moscow,
 Russia), 35 µl of deionized water, 3.4 pmol of each primer, and 5 µl
 original DNA. PCR was performed in a Mastercycler gradient amplifier
 (Eppendorf, Germany) according to the following algorithm: initial
 denaturation for 5 minutes at 95°C, then 35 cycles of 30 seconds at 95°C,
 30 seconds at 55°C, and 2 minutes at 72°C, with a final elongation of 10
 minutes at 72°C. PCR products were purified on agarose gel and eluted
 using NucleoSpin® Gel and PCR Clean-up kit (MACHERY-NAGEL,
 Germany). Sequencing was performed on ABI PRISM 3130 XL
 sequencer using the BIG DYE TERMINATOR kit ver. 3.1 according to the
 manufacturer's protocol at the company “Synthol” (Moscow, Russia).
 The forward and reverse sequences were edited and aligned manually using
 BioEdit 7.0.5.3. (Hall, 1999). The obtained DNA sequences were deposited
 in the GenBank database (ON406420, ON416865-67).
 Data matrices were analyzed by
 the maximum parsimony method using the algorithm described in Templeton et
 al. (1992) and implemented in the TCS v. 1.21 (Clement et al., 2000). The
 method evaluates the unrooted network of haplotypes and all links between
 haplotypes with a probability of 95%. Only nucleotide substitutions were
 taken into account in the analysis, indels were considered as missing
 data. Sequences of other <em>Globularia</em> species available
 from GenBank were included in the analysis.