Introduction

Natural environments can be subjected to modifications by anthropic interventions, natural disasters, and climate changes, among other factors. In particular, climate change will likely have an impact on plant ecophysiology, distribution, and community interactions (IPCC 2014). To face climate challenges, the presence of variability within and among natural populations is fundamental because it can provide genetic flexibility in changing environments (Nicotra et al. 2010). In sessile organisms such as plants, individual genotypes could optimise their performance through phenotypic plasticity, a phenomenon in which alternative phenotypes are expressed by the same genotype in diverse environments, and whose potential in plant evolution is being actually considered (Nicotra et al. 2010; Cooper et al. 2019; Kelly 2019). Phenotypic plasticity allows individuals to rapidly cope with environmental fluctuations (Kooke et al. 2015; Nicotra et al. 2015), which could result from climate change.

Phenotypic diversity originates from genetic variation (i.e., allele diversity as a result of differences in DNA sequence; Medrano et al. 2014) as well as from epigenetic variation (phenotypic changes unrelated to variation in DNA sequence; Richards 2006; Medrano et al. 2014). Cytosine DNA methylation, the most widely studied epigenetic mechanism, can be modulated in response to stressful conditions (Verhoeven et al. 2010; 2016) and can impact ecologically important phenotypic traits (Cubas et al. 1999; Marfil et al. 2009, 2012; Zhang et al. 2013; Cortijo et al. 2014).

In numerous plant studies carried out under controlled artificial conditions, the role of epigenetic variation has been evaluated as a source of intraspecific diversity (Tricker et al. 2012; Latzel et al. 2013; Zhang et al. 2013). The ecological consequences of this phenomenon are now being considered by performing studies under natural conditions. Epigenetic diversity in natural populations may contribute to affront diverse environmental conditions (Lira-Madeiros et al. 2010; Herrera and Bazaga 2011; Dubin et al. 2015; Richards et al. 2017) and the phenotypic variation caused by epigenetic changes could be a rapid acclimation mechanism to environmental fluctuations (Gao et al. 2010; Nicotra et al. 2015). In higher plants, however, there is little evidence on environmentally induced epigenetic changes and phenotypic plasticity during the lifespan of a single organism as well as on the dynamics of the phenomenon in natural conditions.

One of the most important non-cereal food crop worldwide (FAO UN 2018), the common potato Solanum tuberosum L., has hundreds of wild tuber-bearing relatives (Solanum spp., section Petota subsection Potatoe) distributed throughout the America, which are collectively known as ‘wild potatoes’. Wild potatoes present a great ecological diversity and can be found in even contrasting ecosystems such as deserts, humid forests, high Andean valleys or coastal regions (Spooner 2009). Phenotypic plasticity has been acknowledged as an essential attribute of this plant group (Correll 1962; Hawkes and Hjerting 1969; Masuelli et al. 2009; Camadro 2012), and numerous studies have attempted to shed light on the mechanisms that make this diversity possible. However, studies—carried out with germplasm bank accessions—have revealed that there is a lack of correlation between geographical distribution and genetic and phenotypic variability (McGregor et al. 2002; Jansky et al. 2009). Also, genotypes with different levels of pest and disease resistance coexist in the same populations (Jansky et al. 2006; Marfil and Masuelli 2014). These results would indicate that genetic variability and, possibly, other mechanisms contribute to generate the extensive natural variability of wild potatoes (Marfil et al. 2009; 2012; 2018; Cara et al. 2013).

Solanum kurtzianum Bitter and Wittm. is a diploid wild potato that grows in an arid region in western Argentina ranging in altitude from 1082 to 2490 m a.s.l. (Hijmans et al. 2002). Listed in the Global Priority Crop Wild Relative Inventory, and successfully used in potato breeding (Marfil et al. 2015), it is considered the best-adapted of all Argentinian taxonomic species to dry environments. As other species in subsection Potatoe, it has two alternative modes of reproduction with important roles in natural populations: sexual, through botanical seeds, relevant for populations establishment and as a source of genetic variability; and, asexual, through stolons and tubers, which play a role in population maintenance under several environmental conditions (Marfil and Masuelli 2014; Marfil et al. 2015).

One approach to address whether natural populations are able to respond to climate change is to perform common garden experiments under natural conditions (Berend et al. 2019). Along this line, it is advisable to use genetically uniform backgrounds to obtain evidence about the effect of the environment on both epigenetic variability and phenotypic plasticity. Hence, clonal plants established in experimental gardens (EGs) represent an excellent experimental model for this type of studies. In accordance with these considerations, the present work was performed with three genotypes of S. kurtzianum growing in its natural habitats: the Villavicencio Natural Reserve (VNR). This reserve is an Andean Protected Area in Western Argentina with contrasting environmental conditions, and it has been the centre of previous studies by our research group (Ibañez et al. 2017). The aim of the present work was to evaluate the phenotypic plasticity and the possible establishment of differential DNA methylation patterns and their dynamics in biological replicates (clones) growing at contrasting altitudes along three growing seasons. We aspire to answer: what is the extent at which environmental factors imposed by in situ growing conditions can induce phenotypic and epigenetic changes and, moreover, the stability of these changes if they actually occur.

Materials and methods

Plant material

We evaluated morphological, biochemical, and phenological phenotypes and analysed the dynamics of DNA methylation patterns in clones of three S. kurtzianum genotypes. These clones were cultivated during 3 successive years, both permanently and as reciprocally transplanted, in two Andean EGs, located at 1141 and 2113 m a.s.l. in the VNR. To establish the EGs, and according to the number of available tubers, one genotype from each of the three populations collected in the VNR and studied by Marfil and Masuelli (2014)—1228M, QH5 and 2166M—was selected and labelled as G-1, G-2, and G-3, respectively. Twelve tubers (clones) per genotype were cultivated in pots at the campus of Facultad de Ciencias Agrarias, UNCuyo (33° 00′24″ S, 68° 52′19″ W, 940 m a.s.l.) from November 2013 to January 2014.

EGs and reciprocal transplant approach

When the pot-grown plants presented three fully expanded leaves, they were transplanted into 10-L pots with a mixture of autoclaved local soil and organic compost (3:1), to prevent seed contamination. These pots were used to establish the field experiment in January 2014, in two EGs located in the VNR. One EG (hereinafter called ‘EG-1100’) was located at 1141 m a.s.l. (S32° 34′39.57″ S, 68° 56′45.65″ W), and the other (hereinafter called ‘EG-2100’) at 2113 m a.s.l. (32° 35′06.47″ S, 69° 05′56.85″ W). Half of the clones of each genotype were placed in each EG following a completely random design (Ibañez et al. 2017). Overall, six clones × three genotypes × two EG were evaluated. Plants were weekly irrigated with KSC II PHYT-actyl fertiliser (Timac Agro, USA). In September 2014, tubers formed during the growing season (November 2013 to April–May 2014) were collected and conserved in the laboratory at 4 °C until dormancy breakdown. The sprouted tubers were used as starting material for the next growing season. The cultivation conditions were repeated with the modifications detailed below (Fig. 1).

Fig. 1: Experimental garden approach.
figure 1

Three Solanum kurtzianum genotypes were cultivated in two experimental gardens located within an Andean Protected Area at 1100 and 2100 m a.s.l. (EG-1100 and EG-2100, respectively) during three successive seasons (2015, 2016 and 2017). As an example, only one genotype is shown. For each growing season, tubers (clonal propagules) obtained from the previous growing season were used as starting plant material, and two groups of clones were generated: permanent and transplanted clones.

Full size image

The cultivation conditions in 2015, 2016 and 2017 were modified by incorporating reciprocal transplants between EGs (Fig. 1). In each season, clones were divided into two groups. Three clones per genotype were cultivated repeatedly in the same EG as in the previous season (permanent clones). For example, clones of G-1, G-2 and G-3 cultivated in EG-1100 in 2014 were cultivated again in EG-1100 in 2015. The other three clones were transplanted (transplanted clones) into the other EG (reciprocal transplants). For example, clones of G-1, G-2 and G-3 cultivated in EG-1100 in 2014 were transplanted into EG-2100 in 2015, whereas those grown in EG-2100 in 2014 were transplanted into EG-1100 in 2015 (Fig. 1). In summary, in 2015, 2016 and 2017, three biological replicates of each of the three genotypes and condition (permanent or transplanted) were cultivated in each EG, that is, three permanently grown clones and three transplanted clones of each of the three genotypes were cultivated in each EG, totalling 18 plants per EG.

Environmental characterisation

Two data loggers (iButton Hygrochron, Maxim Integrated, USA) were placed in each EG at 30 cm from ground to register atmospheric temperature (T°) and relative humidity (RH). One additional data logger (iButton Thermochron, Maxim Integrated, USA) was placed in a pot in each EG to register soil temperature. Data were registered every hour from the beginning of the 2016 season until the end of the 2017 season. Also, to characterise each season (from 2015 to 2017), T°, RH and precipitation data from the nearest weather station at “El Plumerillo”, located ~40 km away from the VNR, were used (data provided by the Meteorological Information Centre of Argentina).

Phenotypic characterisation

All variables described below were measured for all the plants cultivated in 2015, 2016 and 2017 at both EGs.

Leaves

The 6th and 7th fully expanded leaves were removed from plants in anthesis, wrapped with aluminium foil and conserved on ice. Once in the laboratory, leaves were photographed and the total leaf area, terminal leaflet area and rachis length were measured with AxioVison software (Carl Zeiss MicroImaging GmbH, Jena, Deutschland); leaflet number was also registered. Leaflet thickness was estimated through the area/weight ratio in two 2 cm2 leaflet discs dried at 40 °C until reaching a constant weight. Photoprotective pigments as ultraviolet absorbing compounds (UVAC) and anthocyanins were determined in the leaflet discs according to Ibañez et al. (2017).

Tubers

At harvest, the total tuber number (TN) and total and mean tuber weight were recorded for each plant, except in the 2017 season in which data from EG-2100 had to be excluded due to the harmful action of stray animals on both substrate and tubers. Harvested tubers were stored at 4 °C and periodically examined for sprouting (dormancy breakdown). Sprouted tubers were maintained at room temperature in darkness, and the sprouting percentage (S%) and sprout length (SL) were recorded every 2 weeks for 60 days.

DNA extraction

Leaves from each genotype were collected when anthesis began. Leaf collection was done at the same time from plants growing in both EGs and repeated in each of the three growing seasons. Sample preparation and DNA extraction and quantification were done according to Marfil et al. (2015).

MSAP analysis in EG approach

MSAP markers were used to estimate epigenetic variability according to Cara et al. (2013). Three selective primer combinations without fluorescence were assayed: EcoRI + AGC/HpaII-MspI + ATC, EcoRI + AGA/HpaII-MspI + ATC and EcoRI + ACA/HpaII-MspI + ATC. The amplification products were separated by electrophoresis on a 6% denaturing polyacrylamide gel, silver stained and scanned for manual scoring. To calculate the error rate, a subset of 12 samples were analysed in duplicate for each season, using the same extracted DNA.

From the EcoRI/HpaII and EcoRI/MspI raw data, each locus could present multistate information. To construct binary epigenetic matrices, each locus in the multistate matrix was split into three subloci (or subepiloci) following the codification ‘Mixed Scoring 2’ described by Schulz et al. (2013): fragments present in both the HpaII and MspI lanes were codified as Type I, fragments present in the HpaII lane but absent in MspI lane were codified as Type II and fragments present in the MspI lane but absent in HpaII lane were codified as Type III. Then, two types of methylation states were considered: non-methylated (Type I) and methylated (Types II + III). In the methylated state, both amplification patterns (II + III) were clustered in one category according to the findings of Fulneček and Kovařík (2014), who indicated that the hemi-methylated state (fragment presence in the HpaII lane) could not be determined with accuracy. The absence of fragments on both lanes was not used for calculating the methylation state because it could be the result of genetic changes.

Data analysis

Environmental data

In the 2016 and 2017 seasons, each EG was characterised according to the monthly mean values of RH, atmospheric T° and soil T°. Then, a comparison between EGs was carried out by ANOVA and Fisher’s LSD with InfoStat (version 2017; Grupo InfoStat, Córdoba, Argentina). Moreover, in order to understand the seasonal variability between EGs, the monthly environmental data were compared between EGs. To this end, the difference in EG-2100 respect to EG-1100 was calculated for each environmental parameter, with negative or positive values indicating, respectively, lower or higher values in EG-2100. The comparison among seasons was made by ANOVA and Fisher’s LSD with InfoStat. To explore the relationships among environmental factors, EGs and seasons, a principal component analysis (PCA) was performed with variance-scaled and centred data with ‘FactoMineR’ and ‘factoextra’ packages in R (Le et al. 2008; Kassambara and Mundt 2017; R Core Team 2018). Two PCA were constructed: one with data from the data loggers placed in each EG, and the other from data provided by the nearest weather station.

Phenotypic data

S% and SL were analysed with a generalised linear mixed model in InfoStat. Days after harvest and growing season were considered random factors, whereas EG and transplanting were considered fixed factors, with logit used as link function for binomial distribution for S% and identity for SL in InfoStat. TN was analysed by using a generalised model in InfoStat, with EG and transplanting as fixed effects and log as link function for the Poisson distribution. The effects of EG and transplanting in the total leaf area, terminal leaflet area, rachis length, leaflet thickness, UVAC and anthocyanin compounds and total and mean tuber weight were determined by ANOVA and Fisher’s LSD comparison in InfoStat. To determine associations between plants from each EG and phenotypic traits, a PCA was performed with centred and variance-scaled data using ‘FactoMineR’ and ‘factoextra’ packages in R. Phenotypic plasticity was quantified for exploring genotype responses to the environment, and the relative distances plasticity index (RDPI) was calculated according to Valladares et al. (2006). In brief, the distance among each trait value for all pair of genotypes (grown in the two different micro-environments) and defined as the absolute value, was calculated as the difference between the values measured in each EG, and the relative distance was then obtained by dividing this difference by the sum of both values. The RDPI could take values ranging from 0 to 1, corresponding to absence of plasticity to maximal plasticity, respectively. To test the consistency of each evaluated trait over time, Spearman rank correlation coefficients were calculated (Garnier et al. 2001) with the package ‘stats’ in R, which compares whether the ranking for each trait remained similar both temporally (between growing season) and spatially (between EG). The rank was independently established for each evaluated genotype and by considering all replicates cultivated in each EG.

MSAP

Three different binary epigenetic matrices were constructed to analyse DNA cytosine methylation: (i) a combined matrix, with the full set of data (three types of subepiloci) coding the presence (1)/absence (0) of fragment types I, II and III, (ii) non-methylated matrix, in which the non-methylated loci were coded as fragments type I = non-methylated state = 1/fragments type II and III = methylated state = 0 and (iii) methylated matrix, in which the methylated loci were coded as fragments type II and III = methylated state = 1/fragments type I = non-methylated state = 0. The combined matrix contains the full set of data, while the non-methylated and methylated matrices were constructed to test the potential functionality of these specific methylation states (Schulz et al. 2013; Preite et al. 2015). In the three codified matrices, ambiguous type-IV loci were coded as zeros. Monomorphic loci and singletons (i.e., fragments present only in one sample) were excluded, thus, all the analyses were performed with polymorphic markers. From these binary matrices, the epigenetic distance matrices were calculated by using Sorensen–Dice coefficient (Sneath and Sokal 1973). In addition, and in order to test the consistency of the performed analysis, the Jaccard coefficient was used to calculate the epigenetic distance from the combined matrix (Supplementary information, Table S1 and Fig. S1). To determine relationships among plants from the EG based on MSAP polymorphism, a principal coordinates analysis was performed with the ‘vegan’ package (Oksanen 2015) in R. An analysis of variance through distance matrices was made with the adonis function in the vegan package to assess the effects of EG, genotype and transplanting for each evaluated season. This function studies differences in the mean epigenetic distance among plants within each source of variability (EG, genotype and transplanting) and allows to examine the significance of those partitions (Oksanen 2015). DNA methylation levels of plants cultivated in the EGs were computed as the observed proportion of methylated and non-methylated loci in the multistate matrix. Differences in the proportion of each methylation state (methylated or non-methylated) were evaluated by a generalised linear model, with log as link function for binomial distribution with InfoStat Software.

Correlations

A simple Mantel test was performed in order to detect linear correlations between distance matrices from epigenetic (combined matrix) and phenotypic variability in the 2015, 2016 and 2017 seasons by using InfoStat. All phenotypic traits were considered in the analysis of the first two seasons, but only leaf traits were considered for the analysis of the three seasons because the tuber phenotypic data of the 2017 season were incomplete.

Results

Environmental differences in EGs among seasons

The RH, soil and atmospheric T° were registered in 2016 and 2017 with data loggers in each EG (Table S2). The 2016 season was colder and wetter than the 2017 one. In fact, on average and when comparing the 2017 vs. the 2016 seasons, the former had lower RH (66.2% vs. 74.1%) and higher atmospheric (13.1 °C vs. 12.0 °C) and soil T° (17.8 vs. 13.9 °C) than the latter (P value < 0.0001).

In both years, the weather in EG-2100 was colder and wetter than in EG-1100 (Table S2). On average, in EG-2100 the mean annual T° was 36.3% and 33% lower, and RH 12.8% and 17% higher than in EG-1100 for 2016 and 2017, respectively. The differences in mean soil T° between EGs were similar to the observed for the atmospheric T° (Table S2).

Monthly differences between EG in the environmental parameters were also compared between the 2016 and 2017 seasons (Table S3). In February (summer in the southern hemisphere), when samples were collected for DNA extraction, differences between the EG-2100 respect to the EG-1100 were lower in 2016 than in 2017 for RH (9.04% vs. 14.14%) and for both atmospheric T° (−6.1 °C vs. −7.1 °C) and soil T° (−4.1 °C vs. −6.3 °C).

The environmental data registered in ‘El Plumerillo’ weather station indicated that the 2015 season was drier and warm (54.9% RH and 17.3 °C) than the 2016 (62.0% RH and 16 °C) and 2017 (55.6% RH and 16.8 °C) seasons. Also, in the 2016 season, the highest accumulated precipitation (338.4 mm) was register in comparison with that of the 2015 and 2017 seasons (116.2 and 180.8 mm, respectively).

To summarise, a PCA was performed with the environmental information (Fig. 2). In 2016, EG-1100 was associated with the highest values of atmospheric T° and the minimum RH, and EG-2100 was associated with the highest mean RH (Fig. 2a). In 2017, EG-1100 was associated with the highest soil temperature and EG-2100 with the maximum RH (Fig. 2a). According to the information of the closest weather station (Fig. 2b), the 2016 season was associated with highest values of RH, 2017 with the highest values of atmospheric T°, and 2015 with intermediate values of both RH and T°, and the lowest value of accumulated precipitation.

Fig. 2: Environmental differences among seasons.
figure 2

Biplot representation of principal component analysis for the environmental variables registered in: a two mountain experimental gardens (EG-1100, turquoise dots and EG-2100, pink dots) with an elevation difference of about 1000 m within an Andean Protected Area during two seasons (2016 and 2017); and in b the Protected Area nearest weather station ‘El Plumerillo’ for 2015, 2016 and 2017 seasons. Analysis based on mean (Mean_RH), minimum (Min_RH) and maximum (Max_RH) relative humidity, mean (Mean_T), minimum (Min_T) and maximum (Max_T) atmospheric temperature, mean (S_Mean_T), minimum (S_Min_T) and maximum (S_Max_T) soil temperature and accumulated precipitation (Acc_PP) variables.

Full size image

Phenotypic variability

The PCA summarises all phenotypic traits measured for each genotype and EG for the three evaluated seasons (Fig. 3). Leaf traits (Fig. 3a) showed differences among seasons and EGs. The 2015 season was associated with the highest values of leaflet number and UVAC compounds for both EGs, with a small difference between them. In 2016, both EGs were associated with the highest values of total and terminal leaflet area, rachis length and anthocyanin compounds. In 2017, both EGs were associated with the highest values of leaflet thickness. For the tuber traits (Fig. 3b) measured in 2016, EG-2100 was associated with the highest values of mean and total tuber weight and TN. In 2017, EG-1100 was associated with the highest values of S% and SL.

Fig. 3: Phenotypic variability in three Solanum kurtzianum genotypes growing within an Andean Protected Area.
figure 3

Principal component representation based on phenotypic variability in three S. kurtzianum genotypes cultivated during three successive seasons in two mountain experimental gardens (EG) with an elevation difference of about 1000 m within an Andean Protected Area. a Biplot representation of principal component analysis for leaf variables: leaflet thickness, 6th leaf area (LA6l), 7th leaf area (LA7l), terminal leaflet area (TLA), 6th leaf rachis length (RL6l), 7th leaf rachis length (RL7l), 6th lateral leaflet numbers (LlN6l), 7th lateral leaflet numbers (LlN7l), anthocyanins and UVACs. b Biplot representation of principal component analysis for tuber traits: sprout length (SL), sprouting percentage (S%), mean tuber weight (MTW), total tuber weight (TTW) and tuber numbers (TN). The phenotypes were characterised in plants cultivated in EG-1100 (turquoise diamonds) and EG-2100 (pink diamonds) in 2015, 2016 and 2017 seasons. In 2017, tuber variables were analysed only in EG-1100 because the EG-2100 data were incomplete. Small symbols are individual plants and large symbols indicate EG mean.

Full size image

A differential response was observed for the same genotypes according to the EG in which they were grown. For example, in the 2015 season, G-1 presented thicker leaves, G-2 presented more tubers and G-3 presented shorter rachis in EG-2100 than in EG-1100 (Table S4). In 2016, G-1 and G-2 presented larger sixth terminal leaflet area, and G-1 and G-3 more anthocyanin compounds in EG-2100 than in EG-1100 (Table S5). In 2017, only G-2 and G-3 presented differences between EGs (Table S6).

Plasticity index values, measured as RDPI, ranged from 0.05 to 0.49, with the only statistical differences among seasons being for SL and S% whereas genotypes also presented different RDPI, ranging from 0.05 to 0.46, with statistical differences detected only for S% and sprout tuber length (Table 1). Leaf traits were less plastic than tuber traits. On average, RDPI for leaf traits was 0.19, 0.16 and 0.21 for the 2015, 2016 and 2017 seasons, respectively. On the other hand, RDPI for tuber traits was 0.19, and 0.33 for the 2015 and 2016 seasons, respectively. The same trend was observed when RDPI was compared among genotypes. RDPI in leaf traits was 0.18, 0.19 and 0.19 for G-1, G-2 and G-3, respectively, whereas RDPI for tuber traits was 0.32, 0.24 and 0.23 for the same genotypes, respectively.

Table 1 Relative distance plasticity index (RDPI) for three Solanum kurtzianum genotypes growing during three successive seasons (2015, 2016 and 2017) in two mountainous experimental gardens with an elevation difference of about 1000 m within an Andean Protected Area.
Full size table

The most frequent result was the lack of consistency of the evaluated traits between seasons, although some of them showed consistency for at least two seasons (e.g., 6th leaf rachis length, leaflet thickness, mean tuber weight and S% for G-1; Table S7).

Epigenetic variability

By using three primer combinations, a total of 182 loci were evaluated, with 0.5% of error rate calculated from 12 duplicated samples. According to the alternative methylation states at each locus (multistate matrix), a total of 306 binary epialleles or subepiloci (i.e., alternative methylation states: non-methylated condition I, and methylated condition II + III) were computed. The number of polymorphic markers varied for each genotype and season (Table S8). In the multivariate analysis using the combined matrix to calculate epigenetic distances, genotypes explained the greatest proportion of the variability in the MSAP profiles (except for 2017). In addition, significant effects of the EGs were detected in two out of three seasons analysed (Table 2). Due to the largest influence of genotypes in the observed epigenetic variability, the three evaluated genotypes were independently analysed, and significant effects of the EGs were also detected in two out of three seasons analysed (Table 3). The three genotypes presented differential methylation patterns (i.e., the combination of multiple methylated and non-methylated states at different loci) according to the EG where they were cultivated in the 2015 and 2017 seasons (Fig. 4). In 2016, however, no epigenetic differentiation between EGs was observed (Fig. 4). These results were further confirmed with the multivariate analysis of the distance matrix (Table 3). The analysis of reciprocal transplants between EGs showed no methylation differences between transplanted and permanent clones for the three genotypes and the three analysed seasons (Table 3). When these analyses were performed separately for the methylated and non-methylated loci, the clustering patterns were similar to those obtained previously (Fig. S2 and Table S9). The significant effects of the EGs on the DNA methylation variability were maintained for the 2015 and 2017 seasons, except for the methylated loci of genotypes 2 and 3 in 2015 (Table S9).

Table 2 Multivariate analysis of variances in three Solanum kurtzianum genotypes cultivated during three seasons in two mountain experimental gardens (EG) with an elevation difference of about 1000 m within an Andean Protected Area.
Full size table
Table 3 Multivariate analysis of variances in three Solanum kurtzianum genotypes cultivated during three seasons in two mountain experimental gardens (EG) with an elevation difference of about 1000 m within an Andean Protected Area.
Full size table
Fig. 4: Epigenetic variability in three Solanum kurtzianum genotypes growing during three successive seasons within an Andean Protected Area.
figure 4

Principal coordinate representation based on a combined matrix of methylated and non-methylated loci derived from methylation-sensitive amplified polymorphism (MSAP) markers evaluated in clones from genotypes 1, 2 and 3 cultivated in two mountain experimental gardens (EG) with an elevation difference of about 1000 m: EG-1100 (turquoise) and EG-2100 (pink) for the 2015, 2016 and 2017 seasons. Each group presented a 0.90 confidence interval.

Full size image

When the DNA methylation levels were analysed, the non-methylated state was the most common in the three assayed genotypes and during the three seasons (Table S10). Moreover, no differences in DNA methylation levels induced by the EG were observed (Table S11).

Epigenetic and phenotypic variability correlations

The simple Mantel test indicated that the epigenetic variability was associated with the phenotypic variability for all the considered traits in the 2015 and 2016 seasons (R = 0.36, P value = 0.0001; Fig. S3a). Also, for the 2015, 2016 and 2017 seasons, the epigenetic variability was associated with the leaf phenotypic variability (R = 0.37, P value = 0.0020; Fig. S3b).

Discussion

The understanding of the molecular basis of environmentally induced changes in the phenotype is of great interest in breeding programs. In addition, the study of these changes in wild species could improve our understanding and management of the consequences of climate change. Our work is an initial cross-disciplinary study between ecology and molecular biology, which aims to understand the effects of a real-world context on morphological, biochemical and phenological phenotypes and epigenetics. Epigenetic modifications can contribute to phenotypic plasticity with particular implications in fluctuating environments (Kooke et al. 2015). Species with asexual reproduction (mitotic inheritance) have been proposed as ideal experimental models to study the interplay of phenotypic plasticity, epigenetic variability and adaptation (Castonguay and Angers 2012; González et al. 2018). However, the temporal and spatial phenotypic plasticity of a given genotype and the epigenetic dynamics during successive cycles of agamic reproduction have been overlooked. The present study was set out in this conceptual framework. First, we worked with a species with, presumably but unmeasured, high phenotypic variability, with agamy as an alternative mode of reproduction (Marfil et al. 2015). Second, we performed the assay in a mountainous region, a scenario that offers pronounced environmental gradients. This in situ experiment enabled us to describe the spatial and temporal dynamics of phenotypic variability and DNA methylation patterns during three cycles of clonal reproduction.

In order to obtain evidence about the environmentally induced epigenetic and phenotypic variation in a context of ecological realism, two EGs with differences of ~1000 m in altitude were used. Depending on the season, these two EGs presented differences between 33% and 36.6% in atmospheric mean T° and between 12.8% and 17% in atmospheric mean RH. With this experiment, we showed that clones from different genotypes underwent significant DNA methylation changes and, also, the expression of significant phenotypic plasticity in responses to the two different natural micro-environments. Also, we found that the methylation patterns were greatly influenced by the genotype, indicating that the environmentally induced methylation changes are genotype-dependent, in agreement with previous results (Dubin et al. 2015).

The phenotypic evaluation showed the key influence of the growing seasons. The 2015 season was characterised by highest levels of UVACs. These photoprotective pigments are a common plant response to increased levels of ultraviolet-B radiation (UV-B; Berli et al. 2013). However, Ibañez et al. (2017) observed that the high levels of UVACs in S. kurtzianum could not only be explained by high UV-B levels. In addition to UV-B, UVACs accumulate in response to cold and are thought to be important for the expression of tolerance against this stress (Schulz et al. 2015). In February 2015, the mean temperature was 10% lower than in the 2016 and 2017 seasons (23.7, 26.3 and 26.4 °C, respectively; P value = 0.0014), likely determining the higher UVAC levels in 2015.

Leaf area of higher plants has been widely studied in relation to diverse environmental factors (Peppe et al. 2011). Wright et al. (2017) found that the combination of climatic factors better explained differences in leaf area than individual climatic factors. These authors observed that in the driest sites, leaf area decreased with temperature during the warmest months. Also, Peppe et al. (2011) observed that leaf size varied in a sensitive way according to water availability more than with temperature differences. In the present work, the wetter season (2016, accumulated precipitation ~350 mm) and the wetter EG (EG-2100) were associated with bigger leaf size, whereas the 2017 season, which was the warmest and driest season, was associated with thicker leaves, a trait that has been correlated with dry ecosystems (Coneva and Chitwood 2018). In addition, the unique and unpredictable effects of the growing season on the phenotype were supported by the lack of consistency in the evaluated traits throughout seasons for each EG.

The differential epigenetic patterns between clones cultivated in the two EGs were observed in two out of the three evaluated seasons. There were several environmental differences among seasons indicating that each year had a distinct effect on the DNA methylation variability. The high precipitation regime in 2016 could likely have determined the lack of epigenetic differentiation between EGs for the three evaluated genotypes. The 2016 season registered the highest precipitation in the previous 4 years (412-mm annual accumulated precipitation). In addition, we found that the differences in RH between EGs were lower in the 2016 than in the 2017 sampling collection month (in February, the differences in RH between EGs were 36% lower in 2016 than in 2017). In the desert ecosystem in which the experiment was done, it is likely that water availability could be the principal factor inducing differential environmental conditions, which in turn, would determine differential epigenetic patterns in the assayed genotypes.

As asexual reproduction circumvents meiosis, a reduced epigenetic resetting would be expected in each cycle (Verhoeven and Preite 2014). However, we found an extensive mitotic reprogramming of methylation patterns principally associated with the growing season. A key question that emerges from this result is whether the documented epigenetic changes were environmentally directed or stochastically induced (see Verhoeven and Preite 2014). If the stochastic changes would have been predominant, a high dispersion among replicates and genotypes cultivated in each season would have been expected. The lack of discrimination between transplanted vs. permanent clones indicates that a programmed epigenetic response was activated in each year. In addition, the close epigenetic similarity between transplanted and permanent clones cultivated in the same EG indicates that the previously established methylation patterns were mostly reprogrammed according to the current environmental conditions. It is important to mention that the methylation reprogramming was induced in a short-time scale. Gao et al. (2010) working with Alternanthera philoxeroidesals also found a rapid restructuring in DNA methylation states and reversible phenotypic responses to particular environmental factors. This lack of stress-memory could be an adaptive response as the stressful conditions are mainly transient and unpredictable (Crisp et al. 2016). The results obtained in the present work suggest that the epigenetic mitotic reprogramming contributes in the potato asexual propagation mode. The rapid response to environmental changes through phenotypic plasticity and methylation reprogramming may possibly allow persistence of potato populations over time (Erazzú et al. 2009; Marfil et al. 2015). This could explain why potatoes were considered very difficult to eradicate in most regions (Hawkes 1954).

We found significant correlations between environmentally induced DNA methylation changes and phenotypic variation. This result indicates that methylation changes in wild potatoes play a role in phenotypic plasticity induced by specific environmental differences. Even though phenotypic plasticity in potatoes is a quality that has been recognised for decades by multiple authors, systematic studies have not addressed this characteristic (see Hawkes and Hjerting 1969). By comparing the RDPI, a significant difference between seasons and among genotypes was associated with tuber traits. In addition, higher values of phenotypic plasticity were found in sprouting behaviour, which is a crucial life-history stage influencing the survival and maintenance of wild S. kurtzianum populations (Marfil et al. 2015). Thus, this high tuber plasticity could lead to the successful establishment and persistence of natural populations under fluctuating environmental conditions. Notwithstanding, the evaluated genotypes also presented differences in tuber plasticity that indicated a possible local selection pressure on this trait. Using RDPI, phenotypic plasticity has been characterised in traits of different species (e.g., Godoy et al. 2012; Wei et al. 2019, Hiatt and Flory 2020). The significance of phenotypic plasticity to different environmental conditions, have often yielded species‐specific conclusions. However, under certain conditions, higher trait plasticity could contribute to the successful spreading of invasive plants (Hiatt and Flory 2020). In our study, because the clustering differentiation was comparable regardless the epigenetic distance matrix used (combined, methylated or non-methylated), the environmentally induced methylation patterns associated to phenotypic variation could not be attributed to a specific methylation state. Thus, it is not advisable to make inferences of functionality based on the different observed methylation states.

Medrano et al. (2014), who worked with the perennial herb Helleborus foetidus, proposed that the epigenetic variation possibly confers this species with the capacity to exploit a broad range of ecological conditions. Similarly, in wild-potato populations, it has been reported that new epigenetic patterns, established after interspecific hybridisation, would enable the hybrids to cope with new environments (Cara et al. 2013). Our results suggest that the epigenetic response to environment experiencing short-term fluctuations allows genotype persistence and promote the maintenance of the genetic diversity of populations. The evolution conducted by selective pressures tends to reduce the within population genetic diversity. However, one of the fundamental characteristics of wild-potato populations is their great genetic variability within populations (Jansky et al. 2006, 2009; Spooner et al. 2009; Marfil and Masuelli 2014; Marfil et al. 2015). The environmentally induced phenotypic and epigenetic variation could be contributing a buffering effect against selective processes in nature (O’Dea et al. 2016; Donelson et al. 2017). Also, our study contributes to understand how environmental exposure induces reversible epigenetic modifications that alter plant phenotypes, which not necessarily lead to genomic fixation of environmentally induced traits in later generations. However, since the most common methylation state found in this work was non-methylated, it is possible that transposable elements are activated, inducing genotypic changes that are inherited through the germ line with an increment of genetic variability (see Pimpinelli and Piacentini 2019).

In summary, we demonstrated that in an herbaceous perennial species with asexual reproduction and in an Andean desert with stochastic environmental changes principally associated with precipitation levels, the phenotype and methylation patterns are rewritten every year. The results obtained over three seasons did not show any evident memory effect of the previous environment, but a reprogramming of methylation patterns and high phenotypic variation according to the growing conditions. These reprogramming in methylation patterns and phenotype, in turn, could contribute to genotype persistence, preventing the loss of genetic variability in potato natural populations.