Uncovering the genomic basis of phenological traits in Chouardia litardierei (Asparagaceae) through a genome-wide association study (GWAS)

Chouardia litardierei (Asparagaceae) is a non-model, perennial species characterized by exceptional ecological plasticity. In this research, we studied the genetic architecture underlying several phenological traits in selected ecologically diverged populations of this species. We conducted a genome-wide association study (GWAS) to identify genomic regions linked to the following populations-specific phenological traits: Beginning of Sprouting (BOS), Beginning of Flowering (BOF), Flowering Period Duration (FPD), and Vegetation Period Duration (VPD). Combining phenological data from a common garden experiment with an SNP dataset obtained through the ddRAD-seq approach, we identified numerous loci associated with these traits using single- and multi-locus GWAS models. Narrow-sense heritability estimates were high for all traits, with the VPD trait showing the highest estimate (86.95%), emphasizing its importance for local adaptation. Functional annotation of associated genomic regions revealed key protein families involved in flowering time regulation, vegetative growth timing, and stress adaptation. These findings provide insights into the molecular mechanisms of local adaptation in C. litardierei’s populations from different habitats, emphasizing the role of genetic factors in phenological trait variation and ecological divergence across populations.

Introduction

Understanding the genetic basis of phenotypic variation is essential for evolutionary biology, as it elucidates mechanisms underlying speciation, biogeographical distributions, and fitness in natural populations (Savolainen et al., 2013; Mckown et al., 2014). Natural selection acts on allele frequencies within populations, shaping their variation and promoting adaptive traits that enhance survivability and reproductive success (Hu et al., 2020; Walter et al., 2022; Lee et al., 2023). As populations undergo local adaptation, ecological speciation may lead to the emergence of new ecotypes (Turesson, 1922; Todesco et al., 2020) — genetically distinct populations of the same species well-adapted to specific ecological niches (Rundle and Nosil, 2005; Cortés et al., 2018). Although the role of ecotypes in the speciation process remains debated (Lowry, 2012; Fernández-Meirama et al., 2022), several studies highlight their importance in driving genetic divergence along ecological gradients (Lowry et al., 2008; Brandrud et al., 2017; Cortés et al., 2018; Bakhtiari et al., 2019). Rapidly evolving lineages in heterogeneous environments offer valuable insights into the genetic mechanisms driving adaptation and speciation (Feder et al., 2011; Cortés et al., 2018).

Phenology is one of the key features of plants as sessile organisms. It determines the timing of life cycle phases and the duration of growth and reproduction (Schwartz, 2003). Although other factors like photoperiod (Adole et al., 2019; Wang et al., 2020), water availability (Zhou et al., 2024), or selection by pollinators (Sandring and Ågren, 2009) may play an important role as well, temperature is considered to be the environmental element with the most substantial impact on various phenological traits (Schwartz, 2003; Cook et al., 2012). Matching the growth and especially reproduction periods with the optimal environmental conditions is of exceptional evolutionary importance and is strongly influenced by natural selection (Duputié et al., 2015). Among phenological traits, flowering time is particularly sensitive to environmental factors, marking a critical transition from vegetative growth to reproduction (Hill and Li, 2016; Gaudinier and Blackman, 2020). In seasonally variable habitats, where the timing and duration of the vegetational season differ across landscapes, plants must initiate sprouting and flowering within a constrained annual timeframe (Anderson et al., 2012). Therefore, the regulation of flowering time emerges as a frequent target of evolutionary processes (Gaudinier and Blackman, 2020). Ecologically divergent taxa in numerous lineages often have different flowering times (e.g., Heslop-Harrison, 1964; Grant, 1981; Levin, 2000) suggesting that some niche shifts were predicated upon temporal change (Levin, 2006). Consequently, alterations in flowering schedules may allow populations to better exploit different groups of pollinators (e.g., Waser, 1983; Goldblatt and Manning, 1996; Johnson et al., 1998), while movements into new pollinator niches are accompanied by changes in floral attributes (Levin, 2006). Natural selection generally favours bigger individuals at maturity; however, the timing of flowering presents a trade-off between maximizing fecundity and ensuring reproductive completion before adverse conditions, such as drought or winter, occur (Anderson et al., 2012). Species facing water limitations often adjust their flowering phenology to align with peak moisture availability, taking advantage of optimal conditions (Settele et al., 2016). For example, Schmalenbach et al. (2014) found that late-flowering Arabidopsis plants coped better with drought by compensating for early growth losses with later recovery, while early-flowering plants, which may flower sooner to exploit available moisture before drought, exhibited lower fitness under the same conditions. High salinity also impairs plant growth in Arabidopsis, acting as a suppressive factor that delays flowering time (Li et al., 2007; Lee et al., 2023). Coupled with variation in mating opportunity, temporal variation in sexual phases of individual flowers may have a significant impact on reproductive success in dichogamous plants (Sargent and Roitberg, 2000). Since phenological traits display extensive variations in plants and are often related to local adaptation (Rathcke and Lacey, 1985), the analysis of their genetic background presents a great opportunity to study the mechanisms of the adaptive divergence process.

Investigating the genomic underpinnings of specific traits within the framework of environmental dynamics is essential for uncovering the mechanisms driving local adaptation and elucidating the complex relationship between phenological traits and adaptive responses (Bernatchez et al., 2023). Although much of our understanding of flowering regulation and vegetation duration derives from studies on model organisms such as A. thaliana (Engelmann and Purugganan, 2006; Kinmonth-Schultz et al., 2021), significant advancements have also been made in agriculturally important species (e.g., Molla, 2022; Vicentini et al., 2023; Flohr et al., 2017; Song et al., 2023). However, broadening this research beyond model organisms could increase our understanding of the diverse genetic mechanisms governing phenological variation in populations of wild, non-model species facing different ecological pressures in their habitats.

Here, we investigated the genetic basis of selected phenological traits in the amethyst meadow squill, Chouardia litardierei (Breist.) Speta; a small, bulbous, perennial species belonging to the Asparagaceae family [following the APG III system (Bremer et al., 2009)]. Being a typical geophyte, C. litardierei plants undergo a dormancy period, which usually stretches from mid-summer to late autumn or early spring, depending on the individual season’s properties. During the spring, soon after the development of young leaves, inflorescence emerges. From late April to early June, depending on the population’s location, the flowering phenophase will unfold, shortly followed by fruiting, which marks the beginning of dying back to an underground perennating organ, i.e., a bulb. C. litardierei produces a large racemose inflorescence, typically consisting of several dozen radially symmetrical flowers, without any apparent morphological adaptations for specific pollination mechanisms. While this has not been formally studied, it is expected to be an open-pollinated species (pers. obs.). In addition to sexual reproduction, it propagates clonally through the formation of bulbs surrounding the central bulb. C. litardierei populations are found across the Dinaric Alps in the western parts of the Balkan Peninsula (Ritter-Studnička, 1954; Gaži-Baskova, 1962). Throughout this region, populations inhabit highly contrasting habitats, thus indicating a very pronounced ecological plasticity of the species (Figure 1).

Figure 1

Figure 1. Habitat types of the studied Chouardia litardierei populations, shown from top to bottom: (A) Karst poljes meadows (locality of Budoške Bare population), (B) Dry mountainous grasslands with exposed bedrock (Lovćen), and (C) Saline coastal marshes (Vrana Lake).

In terms of habitat types, the most substantial contrast can be observed between southernmost populations, which are found on patches of exposed dolomite bedrocks or dry mountainous grasslands with very thin and sparse soil cover, on one side, and populations occupying lush meadows of karst poljes, enclosed depressions with deep and fertile soils, abundant in water, on another. These groups of populations cope with very different types of challenges. For the first group of the populations, the most substantial adaptation pressure is expected to come from limited resource availability accompanied by pronounced seasonality in water availability and temperature, which are usual for such a habitat (Mota et al., 2021). At the same time, the second group faces seasonal flooding that can last up to seven months each year (Mihevc et al., 2010; Bonacci, 2014). In addition to these two prevailing groups of populations, based on a habitat type, a third and the smallest group can be recognized, the one inhabiting deep-soiled marshes along the coastline in western parts of the species distribution range. These populations situated in proximity to the seashore are experiencing different climates [i.e., Cfa and Csb climate types according to Köppen classification (Köppen, 1918)] than other inland meadow-habitat populations, which, for the most part, are found in habitats characterized by Cfb type of climate. In addition, these seashore populations are exposed to periodical sea flooding, which causes an increase in soil salinity, one of the major factors in plant ecology (Bui, 2013). Nonetheless, it is essential to note that although this issue was already addressed by Šilić (1990), no clear differentiation, either phenotypic or genetic, among these groups of populations has yet been recognized.

To learn as much as possible about the genetic background of phenological traits of selected C. litardierei populations from across its distribution range and from different habitats, results from a common garden experiment and genotyping were processed through a set of comprehensive single- and multi-locus genome-wide association (GWA) models. Functional annotation of recognized candidate loci was further performed, thus enabling us to deepen our understanding of the complex genetics behind the phenological aspect of adaptive divergence and to analyze the extent to which differentiation of the studied populations has advanced.

Methods

Plant material, common garden experiment, and phenotyping

To establish the common garden experiment, 214 individuals were relocated from nine chosen populations of C. litardierei. Three populations were selected to represent each of the three presumed groups of populations from different habitat types, as illustrated in Figure 1.

During the sampling expeditions, 22 to 25 individuals were selected from each population, ensuring a minimum distance of 10 meters between them, following the 1:20 rule (Wagner, 1995). The geographic coordinates of the sampling locations are listed in Supplementary File 1. Leaf material from each individual was collected for DNA extraction and desiccated using silica gel. Each sampled individual (represented by a single bulb) was transplanted into a separate two-litre plastic container filled with a mixture of soil, sand, and perlite. The containers were placed in raised beds outdoors, creating a common garden setup that exposed the plants to a temperate continental climate (Cfb climate type) (Köppen, 1918; Zaninović et al., 2008). No additional interventions, such as supplemental watering or pesticide use, were applied, allowing the plants to grow under natural, undisturbed environmental conditions.

After two vegetative seasons of acclimatization, four phenological traits were selected for further research (Table 1): (i) Flowering Period Duration (FPD), calculated as the number of days from the appearance of the first flower to the last; (ii) Vegetation Period Duration (VPD), measured as the time from sprouting to the opening of the first capsule with ripened seeds, also in days; (iii) Beginning of Flowering (BOF), recorded using the earliest plant flowering dates a reference; and (iv) Beginning of Sprouting (BOS), noted by referencing the sprouting date of the first individual. All traits examined were measured on an individual genotype level and were considered polygenic.

Table 1

Table 1. Descriptive statistics of the Chouardia litardierei phenological traits examined in the study.

To assess differences in phenological traits across individual populations and three groups of populations originating from different habitat types, Kruskal-Wallis tests were implemented in the PAST software (Hammer et al., 2001), were performed. We further performed pairwise comparisons using Mann-Whitney post-hoc tests with Bonferroni correction to identify significant trait variations. Since none of the variables followed a normal distribution, Spearman’s correlation analysis was conducted to examine the relationships between FPD, VPD, BOF, and BOS variables using the “stats” package in R (R Core Team, 2016).

Sequencing, genomic data processing, and population genetic structure

DNA isolation was carried out using the GenElute™ Plant Genomic DNA Miniprep Kit (Sigma–Aldrich^®). DNA concentrations were measured with the Qubit™ Fluorometer (Thermo Fisher Scientific, Wilmington, DE, USA), and samples were subsequently diluted to a concentration of 20 ng/μL.

For genotyping the studied C. litardierei populations, a ddRAD-seq approach was utilized (Peterson et al., 2012). DNA was initially digested with two restriction enzymes, AseI and NsiI (NEB #R0526L and #R0127L, respectively). The resulting fragments were then ligated with barcoded i5 and i7 adapters, allowing all samples to be multiplexed. Final amplification was carried out after nick repair using DNA polymerase I (NEB #M0209L). The resulting DNA libraries were double-sequenced (150 bp paired-end) on the Illumina HiSeq X platform.

The initial sequencing data underwent preprocessing for quality trimming and adapter removal using Trim Galore (Martin, 2011). Post-trimming, BAM files were generated by aligning the reads to the C. litardierei reference genome (Radosavljević et al., 2023) using the Burrows-Wheeler Aligner (Li and Durbin, 2009). SNP identification was performed with the Stacks software package v1.48 (Catchen et al., 2013). The ref_map.pl wrapper module was utilized, following Paris et al. (2017) recommendations, the pstacks module was executed to extract loci previously aligned to the reference genome, with a minimum coverage depth of three reads to ensure a reliable representation of loci across samples and reduce low-confidence genotype calls. The cstacks module then constructed a comprehensive catalogue of loci across populations, allowing a maximum of four mismatches among sample loci to minimize alignment errors. Subsequently, the populations module calculated population-level summary statistics. To ensure high data quality, loci were retained only if present in all nine populations and at least 70% of individuals within each population, with a maximum observed heterozygosity of 0.70. Additional filtering criteria included retaining only one SNP per locus and excluding loci with minor allele frequencies (MAF) below 1%. This stringent filtering approach focused on common and well-represented genetic variants, reducing the risk of inaccuracies due to sequencing or sampling errors. The resulting dataset, comprising high-quality genetic markers, was exported in .vcf format for downstream analysis.

To assess the neutral population genetic structure of the studied populations, we used a model-based clustering method implemented in ParallelStructure (Pritchard et al., 2000; Besnier and Glover, 2013). To overcome the issue of this analysis’s high computational demands and lengthy processing time for such a large number of SNPs, we constructed a subset of 5,000 randomly selected SNPs. The analysis comprised ten runs for each of the ten clusters (K). Each run consisted of a burn-in period of 50,000 steps, followed by 500,000 Monte Carlo Markov Chain (MCMC) replicates. We used the StructureSelector online software (Li and Liu, 2018) to obtain the most likely number of clusters (K) following Evanno’s method (Evanno et al., 2005) as well as to retrieve the final data through the clustering and averaging of the runs following the Clumpak algorithm (Kopelman et al., 2015). The obtained results were processed using CorelDRAW X7 v.17.1.0.572 software (Corel Corp., Ottawa, Canada) for improved visualization.

Genome-wide association analyses

Figure 2 illustrates a schematic representation of the methodological approach used in this study. All traits were treated as polygenic and GWAS analyses were carried out assuming an additive genetic model. Variants with a minor allele frequency (MAF) below 1% were excluded using the BCFtools software (Danecek et al., 2021). Two distinct statistical approaches were employed for each association analysis: the frequentist single-locus approach and the Bayesian multi-locus approach. In the frequentist single-locus approach, two distinct models were applied. A standard linear mixed model (LMM) was fitted using GEMMA 0.98.5 (Zhou and Stephens, 2012) for all four traits, keeping in mind that this approach assumes a normal trait distribution. Additionally, all traits were analyzed using GMMAT 1.4.2 (Chen et al., 2019), applying a Poisson generalized linear mixed model (GLMM), to account for their count-based distributions. The Poisson GLMM in GMMAT was selected because it effectively accounts for the non-normal distribution of count data, providing a complementary approach to the LMM analysis performed in GEMMA.

Figure 2

Figure 2. A schematic outline of the methodological approach employed to study the genetic basis of phenological traits in Chouardia litardierei.

In the Bayesian multi-locus approach, a Bayesian sparse linear mixed model (BSLMM) (Zhou et al., 2013) was simultaneously fitted for all traits under analysis. Significant SNPs for each trait were identified by first intersecting the sets of significant SNPs obtained from GLMM and LMM, and then further intersecting the resulting set with those identified by BSLMM, ensuring consistency across both the frequentist and Bayesian approaches (Figure 2). Additionally, a multivariate linear mixed model (mvLMM) was performed in GEMMA to simultaneously analyze significantly correlated traits (FPD and VPD, as well as BOF and VPD) to identify shared association signals between them.

The results were visualized using Manhattan plots generated with the R package “qqman” (Turner, 2018) and “CM plot” (Yin et al., 2021). An ad hoc threshold of 1×10⁻³ was used for the frequentist GWAS analyses (GLMM, LMM, and mvLMM).

Generalized linear mixed model using a poisson distribution

The generalized linear mixed model (GLMM) with a Poisson distribution was applied using GMMAT, and the model is expressed as follows (Equations 1–3):

$\begin{array}{ll}log({\mu }_i)=W_i\alpha +x_i\beta +u_i& (1)\end{array}$

$\begin{array}{ll}u\sim MV{N}_n(0,\lambda K)& (2)\end{array}$

$\begin{array}{ll}{y}_i\sim Poisson({\mu }_i)& (3)\end{array}$

In this model, y_i​ represents the observed count for the i-th individual, while μ_i denotes the mean count, modeled as the exponential of the linear predictor. ​W_i is the i-th row of an n × c matrix of covariates (fixed effects), α is the corresponding vector of coefficients for these covariates, x_i represents the genotype of the i-th individual, and β denotes the effect size of the genetic marker. The random effects u are assumed to follow a multivariate normal distribution MVN_n (0,λK), where K is the relatedness matrix of size n × n, and λ represents the ratio of variance components. The observed data y_i​ is assumed to follow a Poisson distribution with μ_i​. This model incorporates individual-level random effects and a genetic relationship matrix K to account for population structure and relatedness. When assuming a normal distribution and an identity link function for continuous traits, GMMAT conducts association tests using linear mixed models (LMMs).

Linear mixed model

The standard LMM was applied using GEMMA 0.98.5. in the following form:

$\begin{array}{ll}y=W\alpha +x\beta +u+\text{ ϵ}& (4)\end{array}$

$\begin{array}{ll}u\sim MV{N}_n(0,\lambda {\tau }^{-1}K)& (5)\end{array}$

$\begin{array}{ll}ϵ\sim MV{N}_n(0,{\tau }^{-1}{I}_n)& (6)\end{array}$

Here, y represents a vector of trait values for 214 individuals, and W is an n × c matrix of covariates (fixed effects), which, in this case, consists of a column of 1s. Let α represent a c-vector of the intercept, x be an n-vector of marker genotypes, and β denote the effect size of the marker. Additionally, u is an n-vector of random effects, ϵ is an n-vector of errors, τ^-1 denotes the variance of the residual errors, and λ is the ratio between the two variance components. K is the known n × n relatedness matrix, and I_n is an n × n identity matrix. MVN_n refers to the n-dimensional multivariate normal distribution. The effect sizes indicate the change in trait values associated with each additional effect allele in the genotypes of individuals.

Bayesian framework

The LMM (Equations 4–6) implemented in GEMMA evaluates the alternative hypothesis H₁: β ≠ 0 against the null hypothesis H₀: β = 0 for each SNP individually. Extensions of the LMM that account for the effects of variants across multiple loci simultaneously could improve the power to identify causal variants. Bayesian LMMs can model all markers simultaneously by assigning different prior distributions to the marker effects and sampling from their posterior distribution. These Bayesian models, designed for estimating SNP effect sizes, start with a basic linear model that links genotypes X to phenotypes y:

$\begin{array}{ll}y=1_n\mu +X\beta +ϵ& (7)\end{array}$

$\begin{array}{ll}ϵ\sim MV{N}_n(0,\tau -1I_n)& (8)\end{array}$

we let y be a vector of phenotypes observed on n individuals, and X be an n × p matrix of genotypes for these same n individuals at p genetic markers. The vector β represents the effects of genetic markers, 1_n is an n-vector of 1s, µ is a scalar representing the mean phenotype, and ϵ is an n-vector of error terms with variance τ^-1. Our aim was to estimate the parameter β, which corresponds to the effects of the genetic markers. However, because the number of genetic markers in our study (p = 23,315) far exceeds the number of individuals (n = 214), certain modeling assumptions regarding SNP effect sizes β had to be made. These assumptions range from the infinitesimal (or polygenic) model, which posits that all SNPs have non-zero effects, to the sparse model, which assumes that only a small subset of SNPs affect the phenotype. The success of the model relies on the true genetic architecture of the trait being studied, although this is typically unknown. The most widely used polygenic model assumes that all SNPs impact the phenotype (i.e., have non-zero effects) with normally distributed effect sizes:

$\begin{array}{ll}\beta \sim N(0,{\sigma }_{\beta }^2)& (9)\end{array}$

When Equations 7–8 are combined with the normality assumption (Equation 9) for effect sizes b, they result in the previously described LMM, as it incorporates a random effect term that represents the combined genetic effects.

Bayesian sparse linear mixed model

A more general assumption, which includes both polygenic and sparse modeling scenarios, suggests that effect sizes come from a mixture of two normal distributions.

$\begin{array}{ll}{\beta }_i \sim \pi \text{N}(0,\frac{{\sigma }_a^2+{\sigma }_b^2}{\text{pτ}})+(1-\pi )\text{N}(0,\frac{{\sigma }_b^2}{\text{pτ}})& (10)\end{array}$

In this model, π represents the proportion of SNPs with large effects, while σ²_β and σ²_α correspond to the variances of small and large effects, respectively. The resulting BSLMM model combines polygenic and sparse effects in the prior distribution of effect sizes, allowing it to adapt to various genetic architectures of the traits being studied. BSLMM addresses population structure and relatedness by incorporating a genomic kinship matrix as a random effect term, and it accounts for linkage disequilibrium (LD) by estimating SNP effect sizes β while controlling for other SNPs in the model. The model uses a Markov chain Monte Carlo algorithm to sample from the posterior distribution and estimate SNP effect sizes. Unlike LMM, which provides p-values, BSLMM outputs a posterior inclusion probability (PIP) for each SNP, reflecting the likelihood that a marker is associated with the trait based on the data. This PIP is calculated as the proportion of chain iterations in which the SNP exhibits a large effect. SNPs with high PIPs are considered the most likely functional variants influencing the analyzed traits. We applied BSLMM to the same dataset (214 individuals and 23,315 variants) used in our primary frequentist association analysis to compare single-SNP and multi-SNP approaches and reduce false positives. The BSLMM chain was run with 1,000,000 sampling steps and 100,000 burn-in iterations. We used the estimated PIPs from BSLMM for additional fine-mapping of genomic regions identified in the frequentist analysis.

SNP heritability estimation

The proportion of variance in phenotypes accounted for by all available genotypes (PVE), also referred to as narrow-sense heritability (h²), along with the proportion of genetic variance explained by variants with large effects (PGE), was estimated for the traits shown in Table 1. This estimation was based on the assumption that SNP effect sizes follow a mixture of two normal distributions (Equation 10), as implemented in GEMMA BSLMM.

Multivariate genome-wide association analyses

To identify common variants associated with the trait pairs showing the strongest statistically significant correlations, multivariate genome-wide association analyses were performed using a multivariate linear mixed model (mvLMM) in GEMMA. Specifically, multivariate GWAS was conducted for the VPD and BOS traits, as well as for the VPD and BOF traits, which exhibited the strongest statisticaly significant correlations. This approach enabled the simultaneous analysis of genetic effects on both trait pairs of traits by treating them as dependent variables. The mvLMM method accounts for population structure and relatedness among individuals, ensuring accurate identification of genetic variants contributing to the observed phenotypic variation in these traits.

Candidate genes prediction

After identifying phenotypic evidence for local adaptation in distinct C. litardierei populations and conducting GWAS analysis, efforts focused on pinpointing associated candidate genes. Using the reference genome, sequences were extracted spanning a total of 50 kilobases – including 25 kilobases upstream and downstream of each significant SNP identified through both statistical models, using SAMtools (Danecek et al., 2021). Functional annotations for these sequences were then obtained through the eggNOG-mapper v2 database, applying an e-value threshold of < 1 × 10⁻² (Huerta-Cepas et al., 2019).

Results

Phenotyping

Figure 3 illustrates the phenological variations observed among C. litardierei populations in the common garden experiment.

Figure 3

Figure 3. The horizontal bar plot illustrates the durations of Vegetation period Duration (VPD) and Flowering Period (FPD) across populations of the Chouardia litardierei during one vegetational season. The x-axis represents the days of the year, while the y-axis lists the populations being compared. The dark magenta bars indicate the FPD, which represents the duration from the date of the first to the last flower for each genotype. In contrast, the grey bars represent the VPD, denoting the duration from the genotype sprouting to the opening of the first capsule. Additionally, the figure provides a visual reference for the Beginning of Flowering (BOF) and the Beginning of Sprouting (BOS), where BOF and BOS are calculated relative to the individual that flowered or sprouted first, respectively.

Out of the 214 individuals sampled across nine populations, 204 flowered successfully. Consequently, all traits related to flowering [FPD, VPD (since its ending is related to the start of the fruiting phenophase), and BOF] were measured and subsequent analyses were performed on the set of 204 individuals, while the remaining 10 were discarded. At the same time, the BOS trait was analyzed across all 214 individuals. The FPD and the VPD ranged from 9 to 25 days, with a median of 17 days (Q1 – Q3: 15 – 18), and 55 to 162 days, with a median of 97 days (Q1 – Q3: 88 – 107), respectively. The BOF and BOS traits ranged from 1 to 33 days, with a median of 11 days (Q1 – Q3: 10 – 13), and 1 to 88 days, with a median of 56 days (Q1 – Q3: 52 – 63), respectively. All the obtained data are summarised in Table 1. Supplementary File 2 contains the results of Kruskal-Wallis and Mann-Whitney post-hoc tests for the studied phenological traits, showing significant differences at the population level and between the assumed population groups. The distribution of these phenological traits is visually represented using box plots in Figure 4.

Figure 4

Figure 4. Box plots illustrate the obtained phenological results from a common garden experiment, depicting four phenological traits: (A) Flowering Period Duration (FPD) (top left), (B) Vegetation Period Duration (VPD) (top right), (C) Beginning of Flowering (BOF) (bottom left), and (D) Beginning of Sprouting (BOS) (bottom right) per genotype. Each box represents the interquartile range (IQR), with the horizontal line inside the box indicating the median. Whiskers extend to data points within 1.5 times the IQR, while dots represent outliers.

A correlation analysis revealed several significant associations among the studied traits (Table 2). A weak positive correlation was observed between FPD and VPD, while a strong positive correlation was found between VPD and BOF.

Table 2

Table 2. Spearman’s correlation coefficients and p-values for the four C. litardierei phenological traits: FPD, VPD, BOF, and BOS.

Sequencing, genomic data processing, and population genetic structure

The sequencing process generated a total of 1,284,680,304 reads. After filtering the raw sequences and mapping them to the reference genome, 1,278,409,966 reads were retained. SNP identification and filtration were performed using the Stacks software, resulting in the detection of 24,660 SNPs. Following the application of the BCFtools MAF filter with a 1% threshold, 23,315 SNPs were kept for subsequent analysis.

The cluster analysis based on the Bayesian model implemented in the ParallelStructure software revealed that the most likely number of genetic clusters was two (Supplementary File 3). One cluster corresponded to the group of populations from the dolomite bedrock habitat, while the remaining populations formed the other cluster (Supplementary File 4). Such structuring reflects the environmental preferences of the studied populations only to some extent, as populations from seashore and meadow habitats remained grouped without any differentiation among them.

Genome-wide association analyses

The analysis of the FPD trait using LMM identified 48 significant SNPs, while GLMM detected 8. An overlap of these results revealed 8 SNPs that were significant across both methods. Further validation using BSLMM confirmed 3 of these SNPs as significant, with one located on each of chromosomes 10, 7, and 11. For the VPD trait, LMM identified 26 significant SNPs, while GLMM detected 54. Fourteen SNPs were found to overlap between the two methods. Subsequent analysis with BSLMM confirmed 2 of these SNPs as significant, located on chromosomes 4 and 12. In the case of the BOF trait, LMM and GLMM identified 17 and 29 SNPs, respectively, with 8 overlapping SNPs. BSLMM analysis confirmed 1 significant SNP located on chromosome 2. For the BOS trait, LMM identified 34 significant SNPs, while GLMM detected 162. Seven SNPs overlapped between the two methods, and BSLMM analysis confirmed 1 significant SNP on chromosome 12. All SNPs passing the genome-wide significance threshold (1 × 10⁻³) in both LMM and GLMM single-SNP LMM analysis are listed in Table 3. The results from the single-SNP association analysis conducted in GMMAT and GEMMA are presented together in Manhattan plots in Figure 5.

Table 3

Table 3. SNPs passing the genome-wide significance threshold (1 × 10⁻³) in both GMMAT and GEMMA single-SNP LMM analyses for Chouardia litardierei traits: FPD, VPD, BOF, and BOS.

Figure 5

Figure 5. Manhattan plots of single-SNP association mapping of FPD, VPD, BOF, and BOS traits. Single-SNP analysis was conducted using (A) GMMAT (top row) and (B) GEMMA (bottom row) for each trait, where the x-axis represents the chromosomal positions of SNPs and the y-axis shows the −log10 (p-values) from the LMM analysis. The red horizontal line denotes the genome-wide significance threshold (p = 1 × 10⁻³). Each point on the Manhattan plot corresponds to a SNP, with stronger associations appearing higher due to lower p-values. Green dots indicate SNPs identified in both analyses.

In the Bayesian association analysis, two SNPs were identified as having a major sparse effect on the FPD trait. These SNPs were estimated to have a sparse effect in at least 10% of the BSLMM chain iterations (posterior inclusion probability, PIP ≥ 0.099). Additionally, both SNPs showed a sparse effect in over 16% of the iterations (PIP ≥ 0.165), further highlighting their significance. In contrast, for the VPD trait, 75 SNPs displayed a major sparse effect in ≥10% of BSLMM chain iterations (PIP ≥ 0.095). In addition, the top four SNPs displayed a major sparse effect in more than 44% of iterations (PIP ≥ 0.447). Concerning the BOF trait, three SNPs were identified with a major sparse effect in ≥10% of iterations (PIP ≥ 0.098) and the top SNP had a major sparse effect in over 17% of iterations (PIP ≥ 0.172). Similarly, for the BOS trait, 26 SNPs exhibited a sparse effect in ≥10% of BSLMM chain iterations (PIP ≥ 0.095), with the top two SNPs showing a strong effect in over 82% of iterations (PIP > 0.829). The data outlined above is reported in Supplementary File 5.

A total of 7 SNPs passed the genome-wide significance threshold (1 × 10⁻³) in the single-SNP LMM analyses and the posterior inclusion probability threshold (PIP ≥ 10%) in the Bayesian multi-SNP BSLMM analysis and are listed in Table 4. Manhattan plots from the BSLMM analysis are provided in Supplementary File 6.

Table 4

Table 4. SNPs passing the genome-wide significance threshold (1 × 10⁻³) in the single-SNP LMM analyses and the posterior inclusion probability treshold (PIP ≥ 10%) in the Bayesian multi-SNP BSLMM analysis.

SNP heritability estimation

The BSLMM analysis, performed using 23,315 SNPs, provided estimates of narrow-sense heritability (PVE) for the phenological traits studied, along with the proportion of genetic effect (PGE) and the count of variants with a major effect (n.gamma), as detailed in Table 5. The PVE estimate for the FPD revealed that 20.26% of the phenotypic variation was explained by all available genotypes, with 47.22% attributed to 60 SNPs exhibiting significant phenotypic effects. Similarly, the PVE estimate for the VPD indicated that 86.95% of the phenotypic variation was explained by all genotypes, with 65.72% attributed to 111 SNPs exhibiting notable phenotypic effects. Moreover, the BSLMM analysis revealed that 66.03% of the phenotypic variation in BOF was explained by all genotypes, with 25.86% of this variation accounted for by 47 SNPs with significant effects. The PVE estimate for the BOS revealed that 76.05% of the phenotypic variation was explained by all available genotypes, with 63.19% attributed to 52 SNPs exhibiting significant phenotypic effects. Supplementary File 7 contains the means, medians, and 95% equal tail posterior probability intervals (95% ETPPIs) of the hyperparameters derived from the BSLMM.

Table 5

Table 5. Genetic architectures of Chouardia litardierei phenological traits identified using a BSLMM.

Multivariate GWAS analysis

In the multivariate GWAS analysis, 113 SNPs surpassed the genome-wide significance threshold (p = 1 × 10^-3) for the model with BOS and VPD traits as dependent variables (Supplementary File 8). This indicates shared genetic factors influencing these phenological traits across multivariate and univariate analyses. Five SNPs were significant in both LMM and GLMM univariate analyses for the BOS trait, and these same five were also significant for the VPD trait, along with an additional eight SNPs that were significant only for VPD, bringing the total to 13 (Table 6). In the multivariate GWAS analysis for the model with VPD and BOF traits as dependent variables, 36 SNPs exceeded the same threshold (Supplementary File 9). Among these, 10 SNPs were significant in LMM and GLMM univariate analyses for the VPD trait, while 4 showed significance for the BOF trait (Table 6). The multivariate GWAS findings for BOS and VPD, and BOF and VPD are plotted in Manhattan plots in Figure 6. The frequencies of effect alleles across populations for the significant SNPs (shown in Tables 4, 6) are depicted in a plot provided in Supplementary File 10.

Table 6

Table 6. SNPs passing genome-wide significance threshold (1 × 10⁻³) in the multivariate GWAS mvLMM analysis of Chouardia litardierei phenological traits BOS and VPD, and BOF and VPD.

Figure 6

Figure 6. Manhattan plot of multivariate genome-wide association study (multi-GWAS) of (A) BOS and VPD (left) and (B) BOF and VPD traits (right). The red horizontal line indicates the genome-wide significance threshold (p = 1 × 10^-3). Each dot on the Manhattan plot signifies a SNP. The strongest associations have the smallest p-values, so their negative logarithms will be the greatest, appearing higher on the plot. Green dots indicate SNPs identified as significant in the multivariate GWAS analysis as well as in both GEMMA and GMMAT univariate analyses for each of the two plots.

GWAS candidate genes identification

The eggNOG tool provided detailed data clarifying the connection between individual SNPs/sequences and specific protein families (PFAM). To identify candidate genes potentially influencing phenological traits, we conducted eggNOG analysis on 7 SNPs that passed the genome-wide significance threshold (1 × 10⁻³) in both the single-SNP LMM and multi-SNP BSLMM analyses of C. litardierei traits, including FPD, VPD, BOF, and BOS. This analysis identified 59 queries corresponding to sequences matched to the eggNOG database for functional annotation (Supplementary File 11). Using eggNOG, we further analyzed 13 SNPs that met the same significance threshold in the multivariate GWAS analysis of BOS and VPD, uncovering 114 additional queries (Supplementary File 12). Similarly, 14 SNPs passed the same significance threshold in the multivariate GWAS analysis of VPD and BOF, resulting in 173 additional queries (Supplementary File 13). The eggNOG analysis connected sequences to protein families, which we further explored through manual inspection and a literature review to identify specific genes and PFAM domains related to the traits being studied. Some domains were common to both the univariate and multivariate GWAS results, resulting in overlaps. The most significant findings, along with their biological functions and relevant references, are summarized in Table 7.

Table 7

Table 7. List of candidate genes for regions of strong association with FPD, VPD, BOF and BOS identified by the eggNOG-mapper v2 database.

Discussion

This study aimed to advance our understanding of the genetic foundations of phenological adaptive traits in C. litardierei’s populations occupying contrasting habitats and shaped by distinct ecological pressures. To minimize the effects of phenotypic plasticity and identify heritable local adaptation traits as accurately as possible, individuals from divergent environments were grown under uniform conditions (Liu and El-Kassaby, 2019; Schwinning et al., 2022). This approach allowed the separation of genetic influences from environmental effects, revealing the heritable components driving local adaptation, where populations evolve toward optimal phenotypic and genetic configurations in response to local selective pressures (Montejo-Kovacevich et al., 2021).

Basic statistical analyses on the common garden experiment data were first performed to characterize the variations within the tested phenological traits and their potential importance for the local adaptation of studied populations in their natural habitats. Except for the duration of flowering (FPD), substantial variations in tested traits among the studied populations were revealed, highlighting their importance for adaptation to contrasting environmental pressures. However, although significant differences in phenological traits among studied groups and individual populations were present, there were many exceptions in the general pattern. For instance, although the dolomite-habitat population group began with flowering (BOF) before the remaining two groups, the Pag population from the seashore habitat was an exception, as it overlaped with all the dolomite-habitat populations. At the same time, the Pag population came into flowering significantly earlier than the Vrana Lake population, which is found in the same habitat and is even geographically closely positioned to the Pag population. Similarly, VPD was significantly shorter in the dolomite-habitat group of populations in contrast to other groups; however, the Budoške Bare population from karst poljes’ meadow habitat joined this group due to having VPD also significantly shorter than any of the remaining populations from this and the seashore habitat. Such a result supports the earlier assumption that although groups of C. litardierei population thrive in highly contrasting habitats, their differentiation into well-differentiated ecotypes remains poorly supported. This was also partially confirmed by the obtained population genetic results (Supplementary Files 3, 4). Here, only the group of populations from the dolomite habitat was substantially differentiated and formed a well-defined genetic cluster, while all the remaining populations remained clustered together, without signs of differentiation between the seashore and meadow-habitat groups. Since the ecotypes are defined as groups of populations whose differentiation is supported both genetically and phenotypically (Lowry, 2012), the studied groups do not meet these criteria. Nonetheless, some trends can be observed in the obtained results that point to certain conclusions. The dolomite-habitat population group sprouted later (BOS) and flowered earlier (BOF) but had a shorter vegetation period (VPD) than the remaining two groups. Such a shift in phenophases is likely to be significantly influenced by habitat properties. To understand how this specific habitat may affect this phenomenon, two aspects must be considered: the dolomite substrate properties and the influence of the local climate dynamics on the vegetation season. These southernmost populations of C. litardierei are usually found on bare dolomite bedrock or less frequently in dry, exposed mountainous grassland habitats developed on very shallow rendzina soils (Figure 1). Due to reduced water and nutrient capacity, accompanied by high levels of thermal conductivity and thermal capacity of the dolomite substrate (Thomas et al., 1973; Waples and Waples, 2004; Mota et al., 2021), these drought-prone habitats are known to induce heat stress in adjacent organisms and thus present a hostile environment for plant species (Mota et al., 2021). In addition, due to very sparse vegetation cover in such habitats, the substrate temperature can be expected to reach far greater values when compared to habitats covered with canopies or meadows (Oliver et al., 1987), thus further worsening already inhospitable conditions. Regarding the influence of regional climatic patterns on local vegetation, two peaks of ecosystem productivity have been observed in Mediterranean climate conditions across southern Europe – the larger one during spring and the less pronounced one during autumn. Such a modality has developed because of ecological constraints imposed by low winter temperatures on one side and summer droughts on another (Spano et al., 2013; Camarero et al., 2021), thus leaving relatively short time frames in spring and autumn suitable for development and reproduction. Consequently, it seems plausible that populations experiencing such climatic patterns, in combination with drought- and heat-stress-prone habitats, have developed short development- and reproduction-related phenophases. At the same time, the remaining C. litardierei populations inhabiting deep, moisture-retaining soils protected by dense vegetation layer which additionally reduces the increase of substrate temperature (Oliver et al., 1987), experience a less limited time frame for closing the sexual reproduction cycle. This is reflected in significant shifts in related phenophases toward later sprouting and the beginning of flowering, as well as a more extended vegetation period.

By emphasizing their heritable nature, the high PVE values observed in our study were suggested to indicate the great evolutionary importance of detected candidate loci in shaping the phenological adaptation of populations to local climatic conditions. The highest PVE value (86.95%) was exhibited by the trait VPD, suggesting that the length of the growing season in this species is predominantly determined by genetic factors. The high genetic variance observed in VPD could be reflective of adaptive mechanisms that allow C. litardierei to optimize its growth and reproductive success in response to environmental cues, such as climate and soil conditions, with strong natural selection acting on traits critical for survival in fluctuating environments. While a PVE for flowering time exceeding 95% has been reported for Arabidopsis from Cape Verde and Morocco (Neto and Hancock, 2023), highlighting the predominant genetic influence, the PVE for C. litardierei flowering period duration (FPD) was found to be 20.26%, indicating a more significant role of environmental or non-genetic factors. PVE values of 66.03% and 76.05% were exhibited by the BOF and BOS traits, respectively, indicating that genetic elements were exerting a greater influence than local environmental factors in shaping these traits. This was reinforced by the PGE values, with the highest PGE (65.72%) being observed in VPD, driven by a few major variants. In contrast, lower PGE values (25.86% and 63.19%) were found for BOF and BOS, respectively, reflecting the influence of numerous small-effect variants and a greater environmental impact. Overall, these heritability estimates and genetic findings provided evidence of the significant role played by genetic factors in shaping phenological traits in C. litardierei, emphasizing the complex interaction between genetics and environment and offering a strong foundation for future genetic, evolutionary, and adaptation studies.

In this study, multiple loci linked to phenological traits in C. litardierei were identified through univariate and multivariate GWAS approaches. The relatively low overlap of significant SNPs detected across the different GWAS models likely reflects inherent differences in their statistical assumptions and approaches to modelling genetic effects. While both frequentist methods (GLMM and LMM) applied a consistent significance threshold of < 1 × 10⁻³, the BSLMM relies on posterior inclusion probabilities, which are generally more conservative and not directly comparable to p-values. Importantly, each model is optimized for different data characteristics: LMM assumes normally distributed traits, whereas GLMM, using a Poisson distribution, is more appropriate for count-based traits with non-normal distributions. Applying trait-appropriate models increases the reliability and power of association detection, even if it results in a lower number of shared SNPs. Functional annotation of the genomic windows surrounding significant SNP loci revealed regions encoding key protein families involved in essential biological pathways related to phenological events. Among others, SNP loci were identified in regions encoding the chromo domain, which is crucial to plant chromatin-based gene regulation. In Arabidopsis, mutations in LHP1 (LIKE HETEROCHROMATIN PROTEIN 1), a gene encoding a chromo domain, have been shown to cause early flowering and reduced plant size (Gaudin et al., 2001). Overexpression of CONSTANS (CO), which activates FLOWERING LOCUS T (FT) in long-day conditions, has been found to alter chromatin at the FT locus by reducing LHP1 binding and increasing histone acetylation, suggesting LHP1 represses flowering through chromatin regulation (Adrian et al., 2010). SNP loci were also identified in regions encoding histidine phosphatase proteins, which are known to regulate plant development and stress responses, particularly through hormone signaling pathways like cytokinins that influence flowering and vegetative growth (Werner et al., 2001; Hai et al., 2020). For instance, it has been demonstrated that exogenous cytokinin application extends the vegetative phase in rice and maize by inhibiting the expression of florigen genes, such as Hd3a and ZCN8, thus delaying flowering time (Cho et al., 2022). Additionally, cytokinins have been found to interact with environmental signals like nutrient sensing (Argueso et al., 2009; Prasad, 2022), potentially aiding plant adaptation to nutrient-poor and drought-prone habitats, like those inhabited by the southern group of C. litardierei populations. Similarly, cytokinin-deficient mutants have been observed to exhibit delayed flowering on nutrient-poor substrates, underscoring cytokinin’s role in adaptation to nutrient-limited environments (Miyawaki et al., 2006). SNP loci within the genomic regions encoding aspartic proteases (APs) and CCHC-type zinc finger proteins (CCHC-ZFPs) were recognized as well. In drought-susceptible common bean cultivars, the PvAP1 gene exhibited significant upregulation under mild water stress, supporting the role of APs in drought responses (Contour-Ansel et al., 2010). CCHC-ZFPs are considered essential for growth and development, as demonstrated in Arabidopsis, where AtCSP4 has been identified as a key factor (Yang and Karlson, 2011). Additionally, SNP loci within the genomic region encoding pentatricopeptide repeat (PPR) proteins were identified. It has been reported that mutations in the Arabidopsis gene AT1G15480, encoding a P-class PPR protein, result in early flowering (Emami and Kempken, 2019). Furthermore, mutations were detected in genetic regions responsible for encoding phytochrome-interacting factor 1 (PIF1). In Arabidopsis, PIF1 has been found to play a major role in sprouting inhibition (Oh et al., 2004; Yadav, 2024), plant growth and development regulation (Soy et al., 2014), stress adaptation (Li et al., 2024), and regulation of photosynthesis initiation (Chen et al., 2013). In addition, SNP loci were identified within regions encoding the MATE domain, the protein kinase domain, and loci associated with the MLO protein family. Several MATE domain genes in O. sativa (OsMATE2, OsMATE4, OsMATE42, and OsMATE46) have been shown to regulate plant responses to abiotic stresses, such as salt and drought, through differential expression patterns (Du et al., 2021), while the protein kinase SOS2/CIPK24 has been recognized as a central regulator of salt stress response and hormonal signaling in Arabidopsis (Chen et al., 2023). Finally, the MLO protein family is considered crucial for temperature stress adaptation, as exemplified by several OsMLO proteins in O. sativa (Nguyen et al., 2016).

Here, we investigated the genetic background of phenological traits in C. litardierei, revealing significant associations between them and specific genetic variations across the genome. Our findings indicate that certain genomic regions may be instrumental in the adaptive responses of populations to contrasting environmental conditions. The genetic architecture of these phenological traits is complex, with multiple candidate loci contributing to phenotypic diversity across habitats. Using the ddRAD-seq approach and comprehensive GWAS analyses, we identified key candidate genes and multiple loci associated with phenological traits. However, the limited genome scan resolution of ddRAD-seq, particularly in large genomes like C. litardierei (3.7 Gb), leaves much genomic information unexplored. The relatively small sample size is a limitation of our study, particularly given that GWAS typically include larger cohorts to detect robust and reproducible associations. Nevertheless, our analysis revealed several biologically plausible signals, which, while requiring validation, provide a valuable foundation for future studies. These findings should be interpreted with caution, but they offer meaningful insights that can be further explored and confirmed in larger, independent populations. Functional annotation of the associated genomic regions revealed key protein families involved in vital biological pathways related to flowering time, vegetative growth, and stress adaptation. These protein families are crucial regulators of plant development, environmental responses, and abiotic stress adaptation. High narrow-sense heritability estimates indicated that genetic factors accounted for a significant portion of the phenotypic variance, with PVE ranging from 20.26% for flowering period duration (FPD) to 86.95% for vegetation period duration (VPD). This study underscores the complexity of the genetic architecture driving phenotypic diversity in plants, highlighting the critical role of genomic approaches in examining adaptive traits in non-model species exposed to diverse ecological pressures. Despite challenges in studying a wild, non-model species, this research advances our understanding of the genomic basis of adaptive divergence and ecological differentiation in C. litardierei. Expanding this research through a comprehensive Genome-Environment Association (GEA) study, incorporating more populations across the species’ distribution range, could provide deeper insights into the genomic drivers of local adaptation and phenological divergence.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA1166115.

Author contributions

SLŠ: Conceptualization, Formal Analysis, Methodology, Visualization, Writing – original draft, Writing – review & editing. NP: Conceptualization, Formal Analysis, Methodology, Visualization, Writing – original draft, Writing – review & editing. KK: Methodology, Writing – original draft, Writing – review & editing. BS: Conceptualization, Methodology, Writing – original draft, Writing – review & editing. DM: Visualization, Writing – original draft, Writing – review & editing. IR: Conceptualization, Funding acquisition, Methodology, Supervision, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This research was funded by the Croatian Science Foundation. All the research activities were financed by the “Amethyst Meadow Squill (Chouardia litardierei, Hyacinthaceae): a study system for ecological divergence” (HRZZ-IP-2020-02-8099) project.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2025.1571608/full#supplementary-material

Glossary

APs: Aspartic Proteases

BAM: Binary Alignment Map

BOF: Beginning of flowering

BOS: Beginning of sprouting

BSLMM: Bayesian Sparse Linear Mixed Model

CCHC-ZFPs: CCHC-type zinc finger proteins

Chr: Chromosome

ddRAD-seq: Double Digest Restriction Site-Associated DNA Sequencing

EGGNog: Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups

FPD: Flowering period duration

GEA: Genome-Environment Association

GEMMA: Genome-wide Efficient Mixed Model Association

GMMAT: Generalized Mixed Model Association Tests

GWAS: Genome-Wide Association Study

LMM: Linear Mixed Model

MAF: Minor Allele Frequency

mGWAS: Multivariate genome-wide association study

PGE: Proportion of Genetic Effect

PPR proteins: Pentatricopeptide repeat proteins

PFAM: Protein Family

PIF1: Phytochrome-interacting Factor 1

PIP: Posterior Inclusion Probability

PVE: Proportion of Variance Explained

SNP: Single Nucleotide Polymorphism

VPD: Vegetation period duration

References

Adole, T., Dash, J., Rodriguez-Galiano, V., Atkinson, P. M. (2019). Photoperiod controls vegetation phenology across Africa. Commun. Biol. 2, 1–13. doi: 10.1038/s42003-019-0636-7

PubMed Abstract | Crossref Full Text | Google Scholar

Adrian, J., Farrona, S., Reimer, J. J., Albani, M. C., Coupland, G., Turck, F. (2010). cis-regulatory elements and chromatin state coordinately control temporal and spatial expression of FLOWERING LOCUS T in Arabidopsis. Plant Cell 22, 1425–1440. doi: 10.1105/TPC.110.074682

PubMed Abstract | Crossref Full Text | Google Scholar

Anderson, J. T., Inouye, D. W., McKinney, A. M., Colautti, R. I., Mitchell-Olds, T. (2012). Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change. Proc. R. Soc. B: Biol. Sci. 279, 3843–3852. doi: 10.1098/RSPB.2012.1051

PubMed Abstract | Crossref Full Text | Google Scholar

Argueso, C. T., Ferreira, F. J., Kieber, J. J. (2009). Environmental perception avenues: the interaction of cytokinin and environmental response pathways. Plant Cell Environ. 32, 1147–1160. doi: 10.1111/J.1365-3040.2009.01940.X

PubMed Abstract | Crossref Full Text | Google Scholar

Bakhtiari, M., Formenti, L., Caggìa, V., Glauser, G., Rasmann, S. (2019). Variable effects on growth and defense traits for plant ecotypic differentiation and phenotypic plasticity along elevation gradients. Ecol. Evol. 9, 3740–3755. doi: 10.1002/ECE3.4999

PubMed Abstract | Crossref Full Text | Google Scholar

Bernatchez, L., Ferchaud, A. L., Berger, C. S., Venney, C. J., Xuereb, A. (2023). Genomics for monitoring and understanding species responses to global climate change. Nat. Rev. Genet. 25, 165–183. doi: 10.1038/s41576-023-00657-y

PubMed Abstract | Crossref Full Text | Google Scholar

Besnier, F., Glover, K. A. (2013). ParallelStructure: A R package to distribute parallel runs of the population genetics program STRUCTURE on multi-core computers. PloS One 8, e70651. doi: 10.1371/JOURNAL.PONE.0070651

PubMed Abstract | Crossref Full Text | Google Scholar

Bonacci, O. (2014). “Ecohydrology of karst poljes and their vulnerability,” in Dinaric karst poljes - floods for life (EuroNatur, Radolfzell am Bodensee), 25–37.

Google Scholar

Brandrud, M. K., Paun, O., Lorenzo, M. T., Nordal, I., Brysting, A. K. (2017). RADseq provides evidence for parallel ecotypic divergence in the autotetraploid Cochlearia officinalis in Northern Norway. Sci. Rep. 7, 1–13. doi: 10.1038/s41598-017-05794-z

PubMed Abstract | Crossref Full Text | Google Scholar

Bremer, B., Bremer, K., Chase, M. W., Fay, M. F., Reveal, J. L., Bailey, L. H., et al. (2009). An update of the Angiosperm Phylogeny Group classification for the orders and families of flowering plants: APG III. Bot. J. Linn. Soc. 161, 105–121. doi: 10.1111/J.1095-8339.2009.00996.X

Crossref Full Text | Google Scholar

Bui, E. N. (2013). Soil salinity: A neglected factor in plant ecology and biogeography. J. Arid. Environ. 92, 14–25. doi: 10.1016/J.JARIDENV.2012.12.014

Crossref Full Text | Google Scholar

Camarero, J. J., Valeriano, C., Gazol, A., Colangelo, M., Sánchez-Salguero, R. (2021). Climate differently impacts the growth of coexisting trees and shrubs under semi-arid Mediterranean conditions. Forests 12, 381. doi: 10.3390/F12030381

Crossref Full Text | Google Scholar

Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A., Cresko, W. A. (2013). Stacks: an analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140. doi: 10.1111/MEC.12354

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, C., He, G., Li, J., Perez-Hormaeche, J., Becker, T., Luo, M., et al. (2023). A salt stress-activated GSO1-SOS2-SOS1 module protects the Arabidopsis root stem cell niche by enhancing sodium ion extrusion. EMBO J. 42, e113004. doi: 10.15252/EMBJ.2022113004

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, H., Huffman, J. E., Brody, J. A., Wang, C., Lee, S., Li, Z., et al. (2019). Efficient variant set mixed model association tests for continuous and binary traits in large-scale whole-genome sequencing studies. Am. J. Hum. Genet. 104, 260–274. doi: 10.1016/j.ajhg.2018.12.012

PubMed Abstract | Crossref Full Text | Google Scholar

Chen, D., Xu, G., Tang, W., Jing, Y., Ji, Q., Fei, Z., et al. (2013). Antagonistic basic helix-loop-helix/bZIP transcription factors form transcriptional modules that integrate light and reactive oxygen species signaling in Arabidopsis. Plant Cell 25, 1657–1673. doi: 10.1105/TPC.112.104869

PubMed Abstract | Crossref Full Text | Google Scholar

Cho, L. H., Yoon, J., Tun, W., Baek, G., Peng, X., Hong, W. J., et al. (2022). Cytokinin increases vegetative growth period by suppressing florigen expression in rice and maize. Plant J. 110, 1619–1635. doi: 10.1111/TPJ.15760

PubMed Abstract | Crossref Full Text | Google Scholar

Contour-Ansel, D., Torres-Franklin, M. L., Zuily-Fodil, Y., de Carvalho, M. H. C. (2010). An aspartic acid protease from common bean is expressed ‘on call’ during water stress and early recovery. J. Plant Physiol. 167, 1606–1612. doi: 10.1016/J.JPLPH.2010.06.018

PubMed Abstract | Crossref Full Text | Google Scholar

Cook, B. I., Wolkovich, E. M., Davies, T. J., Ault, T. R., Betancourt, J. L., Allen, J. M., et al. (2012). Sensitivity of spring phenology to warming across temporal and spatial climate gradients in two independent databases. Ecosystems 15, 1283–1294. doi: 10.1007/S10021-012-9584-5

Crossref Full Text | Google Scholar

Cortés, A. J., Garzón, L. N., Valencia, J. B., Madriñán, S. (2018). On the causes of rapid diversification in the páramos: Isolation by ecology and genomic divergence in espeletia. Front. Plant Sci. 871. doi: 10.3389/FPLS.2018.01700/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Danecek, P., Bonfield, J. K., Liddle, J., Marshall, J., Ohan, V., Pollard, M. O., et al. (2021). Twelve years of SAMtools and BCFtools. Gigascience 10, 1–4. doi: 10.1093/GIGASCIENCE/GIAB008

PubMed Abstract | Crossref Full Text | Google Scholar

Du, Z., Su, Q., Wu, Z., Huang, Z., Bao, J., Li, J., et al. (2021). Genome-wide characterization of MATE gene family and expression profiles in response to abiotic stresses in rice (Oryza sativa). BMC Ecol. Evol. 21, 1–14. doi: 10.1186/S12862-021-01873-Y/FIGURES/7

PubMed Abstract | Crossref Full Text | Google Scholar

Duputié, A., Rutschmann, A., Ronce, O., Chuine, I. (2015). Phenological plasticity will not help all species adapt to climate change. Glob. Chang. Biol. 21, 3062–3073. doi: 10.1111/GCB.12914

PubMed Abstract | Crossref Full Text | Google Scholar

Emami, H., Kempken, F. (2019). PRECOCIOUS1 (POCO1), a mitochondrial pentatricopeptide repeat protein affects flowering time in Arabidopsis thaliana. Plant J. 100, 265–278. doi: 10.1111/TPJ.14441

PubMed Abstract | Crossref Full Text | Google Scholar

Engelmann, K., Purugganan, M. (2006). The molecular evolutionary ecology of plant development: flowering time in Arabidopsis thaliana. Adv. Bot. Res. 44, 507–526. doi: 10.1016/S0065-2296(06)44013-1

Crossref Full Text | Google Scholar

Evanno, G., Regnaut, S., Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620. doi: 10.1111/J.1365-294X.2005.02553.X

PubMed Abstract | Crossref Full Text | Google Scholar

Feder, J. L., Gejji, R., Powell, T. H. Q., Nosil, P. (2011). Adaptive chromosomal divergence driven by mixed geographic mode of evolution. Evolution 65, 2157–2170. doi: 10.1111/J.1558-5646.2011.01321.X

PubMed Abstract | Crossref Full Text | Google Scholar

Fernández-Meirama, M., Rolán-Alvarez, E., Carvajal-Rodríguez, A. (2022). A simulation study of the ecological speciation conditions in the galician marine snail Littorina saxatilis. Front. Genet. 13. doi: 10.3389/FGENE.2022.680792/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Flohr, B. M., Hunt, J. R., Kirkegaard, J. A., Evans, J. R. (2017). Water and temperature stress define the optimal flowering period for wheat in south-eastern Australia. Field Crops Res. 209, 108–119. doi: 10.1016/J.FCR.2017.04.012

Crossref Full Text | Google Scholar

Gaudin, V., Libault, M., Pouteau, S., Juul, T., Zhao, G., Lefebvre, D., et al. (2001). Mutations in LIKE HETEROCHROMATIN PROTEIN 1 affect flowering time and plant architecture in Arabidopsis. Development 128, 4847–4858. doi: 10.1242/DEV.128.23.4847

PubMed Abstract | Crossref Full Text | Google Scholar

Gaudinier, A., Blackman, B. K. (2020). Evolutionary processes from the perspective of flowering time diversity. New Phytol. 225, 1883–1898. doi: 10.1111/NPH.16205

PubMed Abstract | Crossref Full Text | Google Scholar

Gaži-Baskova, V. (1962). Geografsko raširenje livadnog procjepka ili lučike (Scilla pratensis W. et K.) Biološki Glas 15, 49–54.

Google Scholar

Goldblatt, P., Manning, J. C. (1996). Phylogeny and speciation in Lapeirousia subgenus Lapeirousia (Iridaceae: Ixioideae). Ann. Missouri Bot. Garden 83, 346–361. doi: 10.2307/2399865

Crossref Full Text | Google Scholar

Grant, V. (1981). Plant Speciation. Second edition. New York, NY: Columbia University Press 432.

Google Scholar

Hai, N. N., Chuong, N. N., Tu, N. H. C., Kisiala, A., Hoang, X. L. T., Thao, N. P. (2020). Role and regulation of cytokinins in plant response to drought stress. Plants 9, 422. doi: 10.3390/PLANTS9040422

PubMed Abstract | Crossref Full Text | Google Scholar

Hammer, D. A. T., Ryan, P. D., Hammer, Ø., Harper, D. A. T. (2001). Past: Paleontological Statistics Software Package for Education and Data Analysis. Available online at: http://palaeo-electronica.org/2001_1/past/issue1_01.html (Accessed February 3, 2025).

Google Scholar

Heslop-Harrison, J. (1964). Forty years of genecology. Adv. Ecol. Res. 2, 159–247. doi: 10.1016/S0065-2504(08)60332-3

Crossref Full Text | Google Scholar

Hill, C. B., Li, C. (2016). Genetic architecture of flowering phenology in cereals and opportunities for crop improvement. Front. Plant Sci. 7. doi: 10.3389/FPLS.2016.01906/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Hu, Z. M., Zhong, K., Weinberger, F., Duan, D. L., Draisma, S. G. A., Serrão, E. A. (2020). Linking ecology to genetics to better understand adaptation and evolution: A review in marine macrophytes. Front. Mar. Sci. 7. doi: 10.3389/FMARS.2020.545102/BIBTEX

Crossref Full Text | Google Scholar

Huerta-Cepas, J., Szklarczyk, D., Heller, D., Hernández-Plaza, A., Forslund, S. K., Cook, H., et al. (2019). eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 47, D309–D314. doi: 10.1093/NAR/GKY1085

PubMed Abstract | Crossref Full Text | Google Scholar

Johnson, S. D., Linder, H. P., Steiner, K. E. (1998). Phylogeny and radiation of pollination systems in Disa (Orchidaceae). Am. J. Bot. 85, 402–411. doi: 10.2307/2446333

Crossref Full Text | Google Scholar

Kinmonth-Schultz, H., Lewandowska-Sabat, A., Imaizumi, T., Ward, J. K., Rognli, O. A., Fjellheim, S. (2021). Flowering times of wild Arabidopsis accessions from across Norway correlate with expression levels of FT, CO, and FLC genes. Front. Plant Sci. 12. doi: 10.3389/FPLS.2021.747740/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., Mayrose, I. (2015). Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15, 1179–1191. doi: 10.1111/1755-0998.12387

PubMed Abstract | Crossref Full Text | Google Scholar

Köppen, W. (1918). Klassifikation der Klimate nach Temperatur, Niederschlag und Jahreslauf. Petermanns Geogr. Mitt 64, 193–203.

Google Scholar

Lee, Z., Kim, S., Choi, S. J., Joung, E., Kwon, M., Park, H. J., et al. (2023). Regulation of flowering time by environmental factors in plants. Plants 12, 3680. doi: 10.3390/PLANTS12213680

PubMed Abstract | Crossref Full Text | Google Scholar

Levin, D. A. (2000). The Origin, Expansion, and Demise of Plant Species. New York, NY: Oxford University Press. doi: 10.1093/OSO/9780195127287.001.0001

Crossref Full Text | Google Scholar

Levin, D. A. (2006). Flowering phenology in relation to adaptive radiation. Syst. Bot. 31, 239–246. doi: 10.1600/036364406777585928

Crossref Full Text | Google Scholar

Li, H., Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760. doi: 10.1093/BIOINFORMATICS/BTP324

PubMed Abstract | Crossref Full Text | Google Scholar

Li, Y. L., Liu, J. X. (2018). StructureSelector: A web-based software to select and visualize the optimal number of clusters using multiple methods. Mol. Ecol. Resour. 18, 176–177. doi: 10.1111/1755-0998.12719

PubMed Abstract | Crossref Full Text | Google Scholar

Li, Z. Y., Ma, N., Zhang, F. J., Li, L. Z., Li, H. J., Wang, X. F., et al. (2024). Functions of phytochrome interacting factors (PIFs) in adapting plants to biotic and abiotic stresses. Int. J. Mol. Sci. 25, 2198. doi: 10.3390/IJMS25042198

PubMed Abstract | Crossref Full Text | Google Scholar

Li, K., Wang, Y., Han, C., Zhang, W., Jia, H., Li, X. (2007). GA signaling and CO/FT regulatory module mediate salt-induced late flowering in Arabidopsis thaliana. Plant Growth Regul. 53, 195–206. doi: 10.1007/S10725-007-9218-7/FIGURES/6

Crossref Full Text | Google Scholar

Liu, Y., El-Kassaby, Y. A. (2019). Phenotypic plasticity of natural Populus trichocarpa populations in response to temporally environmental change in a common garden. BMC Evol. Biol. 19, 1–15. doi: 10.1186/S12862-019-1553-6/FIGURES/4

PubMed Abstract | Crossref Full Text | Google Scholar

Lowry, D. B. (2012). Ecotypes and the controversy over stages in the formation of new species. Biol. J. Linn. Soc. 106, 241–257. doi: 10.1111/J.1095-8312.2012.01867.X

Crossref Full Text | Google Scholar

Lowry, D. B., Modliszewski, J. L., Wright, K. M., Wu, C. A., Willis, J. H. (2008). Review. The strength and genetic basis of reproductive isolating barriers in flowering plants. Philos. Trans. R Soc. Lond. B Biol. Sci. 363, 3009–3021. doi: 10.1098/RSTB.2008.0064

PubMed Abstract | Crossref Full Text | Google Scholar

Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12. doi: 10.14806/EJ.17.1.200

Crossref Full Text | Google Scholar

Mckown, A. D., Klápště, J., Guy, R. D., Geraldes, A., Porth, I., Hannemann, J., et al. (2014). Genome-wide association implicates numerous genes underlying ecological trait variation in natural populations of Populus trichocarpa. New Phytol. 203, 535–553. doi: 10.1111/NPH.12815

PubMed Abstract | Crossref Full Text | Google Scholar

Mihevc, A., Prelovšek, M., Zupan Hajna, N. (2010). Introduction to the Dinaric Karst. Postojna: Karst Research Institute at ZRC SAZU. doi: 10.3986/9789612541989

Crossref Full Text | Google Scholar

Miyawaki, K., Tarkowski, P., Matsumoto-Kitano, M., Kato, T., Sato, S., Tarkowska, D., et al. (2006). Roles of Arabidopsis ATP/ADP isopentenyltransferases and tRNA isopentenyltransferases in cytokinin biosynthesis. Proc. Natl. Acad. Sci. U S A 103, 16598–16603. doi: 10.1073/PNAS.0603522103

PubMed Abstract | Crossref Full Text | Google Scholar

Molla, K. A. (2022). Flowering time and photoperiod sensitivity in rice: Key players and their interactions identified. Plant Cell 34, 3489–3490. doi: 10.1093/PLCELL/KOAC230

PubMed Abstract | Crossref Full Text | Google Scholar

Montejo-Kovacevich, G., Salazar, P. A., Smith, S. H., Gavilanes, K., Bacquet, C. N., Chan, Y. F., et al. (2021). Genomics of altitude-associated wing shape in two tropical butterflies. Mol. Ecol. 30, 6387–6402. doi: 10.1111/MEC.16067

PubMed Abstract | Crossref Full Text | Google Scholar

Mota, J., Merlo, E., Martínez-Hernández, F., Mendoza-Fernández, A. J., Pérez-García, F. J., Salmerón-Sánchez, E. (2021). Plants on rich-magnesium dolomite barrens: A global phenomenon. Biology 10, 38. doi: 10.3390/BIOLOGY10010038

PubMed Abstract | Crossref Full Text | Google Scholar

Neto, C., Hancock, A. (2023). Genetic architecture of flowering time differs between populations with contrasting demographic and selective histories. Mol. Biol. Evol. 40. doi: 10.1093/MOLBEV/MSAD185

PubMed Abstract | Crossref Full Text | Google Scholar

Nguyen, V. N. T., Vo, K. T. X., Park, H., Jeon, J. S., Jung, K. H. (2016). A systematic view of the MLO family in rice suggests their novel roles in morphological development, diurnal responses, the light-signaling pathway, and various stress responses. Front. Plant Sci. 7. doi: 10.3389/FPLS.2016.01413/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Oh, E., Kim, J., Park, E., Kim, J., Kang, C., Choi, G. (2004). PIL5, a phytochrome-interacting basic helix-loop-helix protein, is a key negative regulator of seed germination in Arabidopsis thaliana. Plant Cell 16, 3045–3058. doi: 10.1105/TPC.104.025163

PubMed Abstract | Crossref Full Text | Google Scholar

Oliver, S. A., Oliver, H. R., Wallace, J. S., Roberts, A. M. (1987). Soil heat flux and temperature variation with vegetation, soil type and climate. Agric. For Meteorol. 39, 257–269. doi: 10.1016/0168-1923(87)90042-6

Crossref Full Text | Google Scholar

Paris, J. R., Stevens, J. R., Catchen, J. M. (2017). Lost in parameter space: a road map for stacks. Methods Ecol. Evol. 8, 1360–1373. doi: 10.1111/2041-210X.12775

Crossref Full Text | Google Scholar

Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S., Hoekstra, H. E. (2012). Double digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PloS One 7, e37135. doi: 10.1371/JOURNAL.PONE.0037135

PubMed Abstract | Crossref Full Text | Google Scholar

Prasad, R. (2022). Cytokinin and its key role to enrich the plant nutrients and growth under adverse conditions-an update. Front. Genet. 13. doi: 10.3389/FGENE.2022.883924/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Pritchard, J. K., Stephens, M., Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959. doi: 10.1093/GENETICS/155.2.945

PubMed Abstract | Crossref Full Text | Google Scholar

Radosavljević, I., Križanović, K., Šarančić, S. L., Jakše, J. (2023). Towards the Investigation of the Adaptive Divergence in a Species of Exceptional Ecological Plasticity: Chromosome-Scale Genome Assembly of Chouardia litardierei (Hyacinthaceae). Int. J. Mol. Sci. 24, 10755. doi: 10.3390/ijms241310755

PubMed Abstract | Crossref Full Text | Google Scholar

Rathcke, B., Lacey, E. P. (1985). Phenological patterns of terrestrial plants. Annu. Rev. Ecol. Syst. 16, 179–214. doi: 10.1146/ANNUREV.ES.16.110185.001143/CITE/REFWORKS

Crossref Full Text | Google Scholar

R Core Team (2016). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing.

Google Scholar

Ritter-Studnička, H. (1954). Flora i vegetacija livada kraških polja Bosne i Hercegovine. Godišnjak Biološkog instituta u Sarajevu 7, 25–110.

Google Scholar

Rundle, H. D., Nosil, P. (2005). Ecological speciation. Ecol. Lett. 8, 336–352. doi: 10.1111/J.1461-0248.2004.00715.X

Crossref Full Text | Google Scholar

Sandring, S., Ågren, J. (2009). Pollinator-mediated selection on floral display and flowering time in the perennial herb Arabidopsis lyrata. Evolution 63, 1292–1300. doi: 10.1111/J.1558-5646.2009.00624.X

PubMed Abstract | Crossref Full Text | Google Scholar

Sargent, R. D., Roitberg, B. D. (2000). Seasonal decline in male-phase duration in a protandrous plant: a response to increased mating opportunities? Funct. Ecol. 14, 484–489. doi: 10.1046/J.1365-2435.2000.00453.X

Crossref Full Text | Google Scholar

Savolainen, O., Lascoux, M., Merilä, J. (2013). Ecological genomics of local adaptation. Nat. Rev. Genet. 14, 807–820. doi: 10.1038/NRG3522

PubMed Abstract | Crossref Full Text | Google Scholar

Schmalenbach, I., Zhang, L., Reymond, M., Jiménez-Gómez, J. M. (2014). The relationship between flowering time and growth responses to drought in the Arabidopsis Landsberg erecta x Antwerp-1 population. Front. Plant Sci. 5. doi: 10.3389/FPLS.2014.00609/ABSTRACT

PubMed Abstract | Crossref Full Text | Google Scholar

Schwartz, M. D. (2024). Phenology: An Integrative Environmental Science. 3rd Edn, ed. M. D. Schwartz. USA. Milwaukee: Springer Nature Switzerland AG 2024. Available online at: https://link.springer.com/book/10.1007/978-3-031-75027-4 (Accessed January 10, 2025).

Google Scholar

Schwinning, S., Lortie, C. J., Esque, T. C., DeFalco, L. A. (2022). What common-garden experiments tell us about climate responses in plants. J. Ecol. 110, 986–996. doi: 10.1111/1365-2745.13887

Crossref Full Text | Google Scholar

Settele, J., Bishop, J., Potts, S. G. (2016). Climate change impacts on pollination. Nat. Plants 2. doi: 10.1038/NPLANTS.2016.92

PubMed Abstract | Crossref Full Text | Google Scholar

Šilić, Č. (1990). Morfologija, horologija, ekologija i fenologija dviju grupa populacija Scilla litardierei Breistr. (Syn.: S. pratensis Waldst. & Kit. non Bergeret). Bilt. Društva Ekol. BiH Ser. B. 5, 107–116.

Google Scholar

Song, X., Liu, L., Wang, C., Zhang, W., Li, Y., Zhu, J., et al. (2023). Real-time determination of flowering period for field wheat based on improved YOLOv5s model. Front. Plant Sci. 13. doi: 10.3389/FPLS.2022.1025663/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Soy, J., Leivar, P., Monte, E. (2014). PIF1 promotes phytochrome-regulated growth under photoperiodic conditions in Arabidopsis together with PIF3, PIF4, and PIF5. J. Exp. Bot. 65, 2925–2936. doi: 10.1093/JXB/ERT465

PubMed Abstract | Crossref Full Text | Google Scholar

Spano, D., Snyder, R. L., Cesaraccio, C. (2013). “Mediterranean Phenology” in Phenology: An Integrative Environmental Science 3rd Edn, ed. M. D. Schwartz. USA. Milwaukee: Springer Nature Switzerland AG 2024. Available online at: https://link.springer.com/book/10.1007/978-3-031-75027-4. (Accessed January 10, 2025).

Google Scholar

Thomas, J., Frost, R. R., Harvey, R. D. (1973). Thermal conductivity of carbonate rocks. Eng. Geol. 7, 3–12. doi: 10.1016/0013-7952(73)90003-3

Crossref Full Text | Google Scholar

Todesco, M., Owens, G. L., Bercovich, N., Légaré, J. S., Soudi, S., Burge, D. O., et al. (2020). Massive haplotypes underlie ecotypic differentiation in sunflowers. Nature 584, 602–607. doi: 10.1038/s41586-020-2467-6

PubMed Abstract | Crossref Full Text | Google Scholar

Turesson, G. (1922). The genotypical response of the plant species to the habitat. Hereditas 3, 211–350. doi: 10.1111/J.1601-5223.1922.TB02734.X

Crossref Full Text | Google Scholar

Turner, S. D. (2018). qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. J. Open Source Softw. 3, 731. doi: 10.21105/JOSS.00731

Crossref Full Text | Google Scholar

Vicentini, G., Biancucci, M., Mineri, L., Chirivì, D., Giaume, F., Miao, Y., et al. (2023). Environmental control of rice flowering time. Plant Commun. 4, 100610. doi: 10.1016/J.XPLC.2023.100610

PubMed Abstract | Crossref Full Text | Google Scholar

Wagner, D. H. (1995). A Rule of Thumb for Botanists: The 1 in 20 Rule Oregon Flora-On-Line Newsletter. United States: Oregon State University.

Google Scholar

Walter, G. M., Clark, J., Cristaudo, A., Terranova, D., Nevado, B., Catara, S., et al. (2022). Adaptive divergence generates distinct plastic responses in two closely related Senecio species. Evolution 76, 1229. doi: 10.1111/EVO.14478

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, H., Wang, H., Ge, Q., Dai, J. (2020). The interactive effects of chilling, photoperiod, and forcing temperature on flowering phenology of temperate woody plants. Front. Plant Sci. 11. doi: 10.3389/FPLS.2020.00443/BIBTEX

PubMed Abstract | Crossref Full Text | Google Scholar

Waples, D. W., Waples, J. S. (2004). A review and evaluation of specific heat capacities of rocks, minerals, and subsurface fluids. Part 1: Minerals and nonporous rocks. Natural Resour. Res. 13, 97–122. doi: 10.1023/B:NARR.0000032647.41046.E7/METRICS

Crossref Full Text | Google Scholar

Waser, N. M. (1983). “The Adaptive Nature of Floral Traits: Ideas and Eviaence” in Pollination Biology, ed. L. Real (Academic Press), 241–285. doi: 10.1016/B978-0-12-583980-8.X5001-7

Crossref Full Text | Google Scholar

Werner, T., Motyka, V., Strnad, M., Schmülling, T. (2001). Regulation of plant growth by cytokinin. Proc. Natl. Acad. Sci. 98, 10487–10492. doi: 10.1073/PNAS.171304098

PubMed Abstract | Crossref Full Text | Google Scholar

Yadav, A. (2024). Branching paths: Unveiling functional divergence of PIF1 and PIF4 via promoter and protein variations. Plant Cell 36, 2749–2750. doi: 10.1093/PLCELL/KOAE152

PubMed Abstract | Crossref Full Text | Google Scholar

Yang, Y., Karlson, D. T. (2011). Overexpression of AtCSP4 affects late stages of embryo development in Arabidopsis. J. Exp. Bot. 62, 2079–2091. doi: 10.1093/JXB/ERQ400

PubMed Abstract | Crossref Full Text | Google Scholar

Yin, L., Zhang, H., Tang, Z., Xu, J., Yin, D., Zhang, Z., et al. (2021). rMVP: A memory-efficient, visualization-enhanced, and parallel-accelerated tool for genome-wide association study. Genomics Proteomics Bioinf. 19, 619–628. doi: 10.1016/J.GPB.2020.10.007

PubMed Abstract | Crossref Full Text | Google Scholar

Zaninović, K., Gajić-Čapka, M., Tadić, M. P., Vučetić, M., Milković, J., Bajić, A., et al. (2008). Climate Atlas of Croatia: 1961–1990: 1971–2000/Klimatski Atlas Hrvatske, Climate Atlas of Croatia: 1961–1990: 1971–2000 (Zagreb: Državni hidrometeorološki zavod).

Google Scholar

Zhou, X., Carbonetto, P., Stephens, M. (2013). Polygenic modeling with Bayesian sparse linear mixed models. PloS Genet. 9, e1003264. doi: 10.1371/JOURNAL.PGEN.1003264

PubMed Abstract | Crossref Full Text | Google Scholar

Zhou, Z., Feng, H., Ma, G., Ru, J., Wang, H., Feng, J., et al. (2024). Seasonal and vertical patterns of water availability and variability determine plant reproductive phenology. Ann. Bot. 135. doi: 10.1093/AOB/MCAE138

PubMed Abstract | Crossref Full Text | Google Scholar

Zhou, X., Stephens, M. (2012). Genome-wide efficient mixed-model analysis for association studies. Nat. Genet. 44, 821–824. doi: 10.1038/ng.2310

PubMed Abstract | Crossref Full Text | Google Scholar