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.

Morphometric analysis

The material was collected in 2021–2022 from twenty-five G. bisnagarica and four G. trichosantha 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 G. bisnagarica: Volga Upland (VU), Sokskiye Yary (SYA), Bugulma-Belebey Upland (BBU), Obshchy Syrt (OS), Zilair Plateau (ZP), and Stavropol Upland (SU). G. trichosantha 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 G. bisnagarica and 120 images for G. trichosantha. 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 Globularia 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 bilat.simmetry and g.pagen functions from the ‘geomorph’ package (Adams et al., 2021).

The intrapopulation morphological variability of the leaf blade shape was calculated using the morphol.disparity 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 Globularia 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 procD.lm 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 powerTransform and bcPower 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 lmer function from the ‘lme4’ package (Bates et al., 2015). The model parameters were calculated by REML estimation. We calculated conditional and marginal R² (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.

Molecular genetics analysis

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 trnL−trnF intron, the rbcL gene, and the rpS16 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 (G. bisnagarica and G. trichosantha) 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 Globularia species available from GenBank were included in the analysis.