Epigenetic Differentiation of Natural Populations of Lilium bosniacum Associated with Contrasting Habitat Conditions.

Epigenetic variation in natural populations with contrasting habitats might be an important element, in addition to the genetic variation, in plant adaptation to environmental stress. Here, we assessed genetic, epigenetic, and cytogenetic structure of the three Lilium bosniacum populations growing on distinct habitats. One population was growing under habitual ecological conditions for this species and the other two were growing under stress associated with high altitude and serpentine soil. Amplified fragment length polymorphism and methylation-sensitive amplification polymorphism analyses revealed that the three populations did not differentiate genetically, but were clearly separated in three distinct clusters according to DNA methylation profiles. Principal coordinate analysis showed that overall epigenetic variation was closely related to habitat conditions. A new methylation-sensitive amplification polymorphism scoring approach allowed identification of mainly unmethylated (φST = 0.190) and fully CpG methylated (φST = 0.118) subepiloci playing a role in overall population differentiation, in comparison with hemimethylated sites (φST = 0.073). In addition, unusual rDNA repatterning and the presence of B chromosomes bearing 5S rDNA loci were recorded in the population growing on serpentine soil, suggesting dynamic chromosome rearrangements probably linked to global genome demethylation, which might have reactivated some mobile elements. We discuss our results considering our earlier data on morphology and leaf anatomy of several L. bosniacum populations, and suggest a possible role of epigenetics as a key element in population differentiation associated with environmental stress in these particular lily populations.

Introduction

Plants are sessile organisms, which are inevitably and constantly exposed to environmental fluctuations. It is recognized that mutations, genetic drift, and selection are too slow processes to explain a rapid response of plants to stress factors from their habitats. Therefore, epigenetic changes and transgenerational epigenetic inheritance might play an important role in plant response to stress and ecological adaptation (Meyer 2015), especially when fast adaptation has to allow time for the permanently heritable genetic and cytogenetic processes to take over. There is growing evidence for natural variation in DNA methylation profiles between plant populations, ecotypes and accessions, and most of these variations have been associated with ecological factors and environmental stress (Labra et al. 2002; Riddle and Richards 2002; Lira-Vaughn et al. 2007; Woo and Richards 2008; Bossdorf et al. 2010; Lira-Medeiros et al. 2010; Verhoeven et al. 2010). Most of these studies were performed under controlled conditions; however, epigenetic variation and its adaptive potential in natural populations are still insufficiently explored (Angers et al. 2010; Latzel et al. 2013).

The peculiarity of DNA methylation in plant genomes is the addition of a methyl group to cytosines in three different sequence contexts (CpG, CpNpG, and CpNpN motifs), which have significant effects on both mechanisms and function of cytosine methylation. Without a change in DNA sequence, modification of cytosine by methylation can silence transcription of protein-coding genes, which is often reflected on biochemical, physiological, or morphological phenotype (Flores et al. 2013; Herrera and Bazaga 2013). Such epi-alleles can be inherited through mitosis or meiosis and can serve as multigenerational epigenetic memory to secure the survival of plant populations under changing habitat (Kalisz and Purugganan 2004; Vaughn et al. 2007; Johannes et al. 2009; Boyko et al. 2010; Kou et al. 2011; Weigel and Colot 2012; Herrera et al. 2013; Avramova 2015; Kooke et al. 2015). DNA methylation is also involved in silencing of transposable elements (TE), which are abundant in plant genomes. Due to their mobility, TE can cause small as well as large-scale genomic rearrangements, which can result in great divergence in size and organization of genomes, even among closely related species, cultivars, and populations (Lisch 2013). Preferential residence or insertion of TE is found within repetitive ribosomal RNA (rRNA) genes (Jakubczak et al. 1991; Raskina et al. 2004; Belyayev et al. 2005), where they can be involved in movement of rDNA loci and karyotype restructuring (Raskina et al. 2004; Datson and Murray 2006). Genes encoding for 18S–5.8S–26S rRNA (18S rRNA) and 5S rRNA molecules are often organized into separate clusters that are located on different chromosomes or in different loci on the same chromosome. The number and position of 18S and 5S rDNA loci are species specific and as such, these loci serve as markers in exploring karyo-evolutionary trends and/or phylogenetic relationships in plant as well as animal groups of species (Roa and Guerra 2012).

Very often, edaphic factors represent a strong selective force in diversification of plant populations or lineages, and adaptation to new soil habitat can result in changes of morphology, phenology, and physiology at the species level. Serpentine is a specific type of geological substrate (peridotite) characterized by low amounts of essential nutrients (N, P, K), low Ca/Mg ratio, high concentrations of heavy metals, as well as by low fertility, shallowness, and porosity resulting in water deficit (Proctor, 1999). Bareness of serpentine soil and low Ca/Mg ratio are identified as extremely stressful and strong selective factors in plant serpentine specialization (Murren et al. 2006; Turner et al. 2010; Cacho and Strauss 2014), leading to serpentine soil ecotypes, which is possibly the first step in the evolution of serpentine endemism (Kruckeberg 1951; Brady et al. 2005).

Lilium bosniacum (Beck) Beck ex Fritsch is a rare species endemic to the Balkan peninsula (Beck von Mannagetta 1887; Muratović et al. 2005; Muratović, Robin, et al. 2010). This flowering plant is a perennial species, ranging in height from 60 to 80 cm. When growing under optimal conditions, Bosnian lilies flower from late May to early July and form 3–10 Turk’s-cap flowers per stem. The latest source (Euro+Med 2006) confirmed the existence of 20 species of the genus Lilium L. in Europe. The list was extended with mostly regional endemics, as was the case with L. bosniacum (Beck) Beck ex Fritsch (Muratović et al. 2005; Muratović, Robin, et al. 2010). Earlier ecological studies on L. bosniacum showed that its typical environmental conditions are limestone or sometimes dolomite geological substrates, deep soils rich in nutrients, full exposure to sunlight, and an altitude from 1,200 to 1,300 m (Lakušić and Kutleša 1971; Muratović et al. 2005; Muratović, Bogunić, et al. 2010; Muratović, Hidalgo, et al. 2010; Muratović, Robin, et al. 2010). No Lilium species, ecotypes or populations were described as a serpentinophyte taxon. However, in the year 2005, we found a unique small population of L. bosniacum at the microsite Kladanj, growing on serpentine soil, therefore exposed to stress associated with this type of substrate.

In the present study, we investigated genetic, epigenetic, and cytogenetic structure of three natural populations of L. bosniacum from contrasting habitats. We have chosen the following populations: a unique population growing on serpentine soil, at the microsite Kladanj (P01); an alpine population growing under stresses associated with high altitude at the site Jahorina (P02); and the population growing under habitual ecological conditions for this species at the site Crepoljsko (P03). We used amplified fragment-length polymorphism (AFLP) and methylation-sensitive amplification polymorphism (MSAP) in order to estimate the degree of genetic and epigenetic diversification in the three populations. Fluorescence in situ hybridization (FISH) was used to map rDNA loci in karyotypes of the three populations, since preliminary results indicated variability in the number/position of 18S rDNA loci in the Bosnian lily population from the microsite Kladanj, growing on serpentine soil (Muratović et al. 2005).

Materials and Methods

Population Sites and Plant Material

Lilium bosniacum is an endemic and protected species in Bosnia and Herzegovina, so we were careful in collecting plant material from the natural sites. The study was conducted at three natural sites (hereafter called populations) located in southwest Bosnia and Herzegovina (fig. 1, table 1). The populations were small; therefore, the maximal number of plants that were collected at the sites was 15 per population. In general, the ecological conditions for Bosnian lily are altitudes higher than 1,000 m above the see level (m a.s.l.) and a limestone substrate. We have chosen a population growing at the habitual conditions (Crepoljsko; P03), and two populations growing under unfavorable conditions concerning altitude, temperature range, light availability, and type of soil. The population at the site Jahorina (P02) was alpine, growing at a very high altitude of 1,850 m a.s.l. and on limestone substrate. The population Kladanj (P01) was growing at an unusually low altitude of 853 m a.s.l. and on serpentine bare substrate, which is unusual for this species (table 1, fig. 1). This woodland population was very small (therefore we name it a microsite)—plants were growing as an enclave groups of four to five individuals between trees, or were growing as solitaires on clearings, probably due to the extremely steep terrain (Muratović et al. 2005; Muratović, Bogunić, et al. 2010). In this woodland population, the plants had the lowest light availability compared with other two populations. Also, a harsh continental climate conditions—below freezing temperature values during 3 months, severe and late frosts and abundant precipitation—characterized the site of Kladanj. The density of P01 and P02 populations counted about two to three individuals per m2. The straight-line distance between the populations are as follows: P01 and P02—57 km; P01 and P03—36 km; P02 and P03—32 km (fig. 1).

. 1.

—Geographical location and photos of the three studied sites of L. bosniacum. Section of alpine geographic region in central Bosnia and Herzegovina with locations of three sampled populations: P01 microsite Kladanj; P02 site Jahorina, and P03 site Crepoljsko. Symbol colors indicate the habitat type (red = woodland, serpentine soil, limited amount of light; green = meadow, limestone soil, high exposure to sun; and blue = mountain meadow, limestone soil with rocks, high exposure to sun). The straight-line distances between the populations are as follows: P01 and P02—57 km; P01 and P03—36 km; P02 and P03–23 km.

Table 1

Overview of the Sampled Populations of Lilium bosniacum

Population ID	Region	Latitude	Longitude	Habitat Type	Geological Substrate	Altitude (m)	Exposure
P01	Muške Vode near Kladanj	44°14' 11.21″N	18°34' 04.99″E	Coniferous woodland	Serpentine	853	W-SW
P02	Mt. Jahorina	43°43' 08.66″N	18°34' 02.65″E	Subalpine meadow	Limestone	1851	W-SW
P03	Mt. Crepoljsko	43°54' 59.24″N	18°29' 11.10″E	Wet upland meadow	Limestone	1200	E-SE

The bulbs of lily plants were taken from their natural sites and then were grown in pots in the greenhouse of the University Paris-Sud XI (Orsay, France) until the roots were developed enough for karyotype analysis (i.e., until meristematic cells started to proliferate, as could be seen as white meristematic tips). The karyotypes were recorded for 10–15 lily plants per each of the three populations (table 2). For AFLP and MSAP analyses, young undamaged leaves were collected from lily plants at the natural sites and then frozen in liquid nitrogen in situ. The frozen leaf material was kept at −80 °C during transportation and DNA was isolated immediately upon arrival to the laboratory. Leaves were collected at the three population sites during the same season (from May to late June) and at the same phenological stage (fully expanded leaves), in order to exclude possible developmental variation in DNA methylation. The same size of the population sample (15 individuals per population) was taken for genetic (AFLP) and epigenetic (MSAP) analyses.

Table 2

Number of Individual Plants Possessing rDNA Karyotypes A, B, C, and D

Population ID	Region	Total no. of Plants	No. of Plants Karyotype A	No. of Plants Karyotype B	No. of Plants Karyotype C	No. of Plants Karyotype D
P01	Kladanj	10	3	3	2	2
P02	Jahorina	15	15	0	0	0
P03	Crepoljsko	15	15	0	0	0

Chromosome Preparations, Fluorochrome Banding, and FISH

Root tips were pretreated in 0.05% colchicine aqueous solution at room temperature for 5 h and subsequently at 4 °C for 12 h. Roots were then fixed in 3:1 (v/v) ethanol:acetic acid. Meristems were macerated in an enzymatic mixture composed of 3% cellulase R10 (Onozuka Yakult Honsha Co., Tokio, Japan), 1% pectolyase Y23 (Seishin Pharmaceutical Co., Tokyo, Japan), and 4% hemicellulase (Sigma-Aldrich, St. Louis, MO, USA) in 0.01 M citrate buffer (pH 4.6) for 50 min at 37 °C. Chromosome spreads were prepared by gentle squash in a drop of 45% acetic acid. Preparation of nuclear suspension was given in Liu and Whittier (1994), and protocol for releasing DNA fibers from nuclei was given in Parra and Windle (1993).

Chromomycin A3 (CMA) staining was done as described in Siljak-Yakovlev et al. (2002). CMA was used at a concentration of 0.125 mg/ml and chromosomes were stained for 13 min. Simultaneous FISH on metaphase chromosomes and DNA fibers using 18S-5.8S-26S (18S) and 5S rDNA probes was performed according to the protocol already described (Muratović et al. 2005). The 18S and 5S rDNA probes were labeled with Cy3 (Amersham, Saclay, France) and digoxigenin-11-dUTP (Roche Diagnostics, Meylan, France), respectively. For fiber FISH, the probes and DNA fibers were denatured simultaneously at 80 °C for 2 min. Slides were analyzed using an epi-fluorescence Zeiss Axiophot microscope coupled with Princeton CCD camera and Metaview image analyzer.

AFLP and MSAP Analyses

AFLP was used for analysis of genetic structure, while MSAP was used to analyze epigenetic diversity in the three populations of L. bosniacum. We analyzed a total of 45 individuals (15 individuals per population). Liquid nitrogen-frozen leaf tissue was manually ground to a fine powder and then genomic DNA was isolated using the DNeasy 96 Plant extraction Kit (Qiagen, Hilden, Germany), according to manufacturer’s instructions. DNA concentration was measured using a NanoVueTM instrument (GE Healthcare, Little Chalfont, UK). The AFLP analysis was carried out following the protocol described previously with some modifications (Carović-Stanko et al. 2011). Restriction of the genomic DNA was performed using a high concentration restriction enzymes EcoRI (100,000 U/ml, NEB), and Tru1I/MseI (50,000 U/ml, NEB). Specific adapters were ligated using T4 DNA Ligase (Fermentas). Restriction digestion and adapter ligation, pre-amplification, and selective amplifications were performed in a GeneAmp PCR System 9700 (Applied Biosystems, Waltham, MA, USA). Four primer combinations were chosen for selective amplification: VIC-EcoRI-ACG + TruI-CGA, PET-EcoRI-ACC + TruI-CGA, FAM-EcoRI-ACA + TruI-CAC, and NED-EcoRI-AGA + TruI-CAC. In the MSAP reaction, two separate AFLP reactions were performed on each sample, using frequent cutter enzymes with differential methylation sensitivity. The HpaII/MspI isoschizomer pair recognizes the same DNA sequence (5′-CCGG-3′); however, unlike MspI, HpaII cannot cut the internally methylated 5′-CCGG-3′ sequence. Comparison between DNA fragments, observed after EcoRI/HpaII and EcoRI/MspI restrictions, allows identification of the methylation status of the 5′-CCGG-3′ sequence. The MSAP protocol was almost the same as the AFLP protocol except 5 U restriction enzymes HpaII or MspI (50 U/100 U, New England Biolabs, Ipswich, MA, USA) were used (instead of Tru1I) together with analogous HpaII/MspI adaptors and primers. Due to the large genome size of L. bosniacum, which contains ∼67 pg/2C of DNA (Muratović et al. 2005), and consequently a large number of restriction fragments, the two primer combinations with six additional bases in the HpaII/MspI primer for selective amplification were used: VIC-EcoRI-AGCG + HpaII/MspI-TCACAG and FAM-EcoRI-ACA + HpaII/MspI-TCACAG.

Both AFLP and MSAP amplified fragments were separated using capillary electrophoresis in an ABI3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). The amplified DNA fragments in the size range of 50–125 bp were scored using GeneMapper 4.0 software (Applied Biosystems, Foster City, CA, USA) as being present (1) or absent (0) in order to create a binary matrix. For AFLPs, the estimation of error rates was performed as suggested in Herrmann et al. (2010). Three individual DNAs (6.67% of total sample) were duplicated and the output from GeneMapper (fragment sizes and peak heights) was analyzed using the R script “scanAFLPv1-3.r” (Herrmann et al. 2010).

Data Analysis

In MSAP analyses, the banding patterns of EcoRI/HpaII and EcoRI/MspI were classified as follows: (i) fragments present in both profiles (1/1): unmethylated state, (ii) fragments present only in the EcoRI/MspI profiles (0/1): hemi- or full-methylation of internal cytosine (HMeCG or MeCG sites), (iii) fragments present only in the EcoRI/HpaII profiles (1/0): hemi-methylation of external cytosine (HMeCCG sites), and (iv) fragments absent in both profiles (0/0): uninformative state. We used the “Mixed-Scoring 2” approach as implemented in R script “MSAP_calc.r” (Schulz et al. 2013) to generate a matrix of binary data for each of the three subepiloci: u-subepilocus (coding the unmethylated sites), m-subepilocus (coding the methylated sites: HmeCG or MeCG sites), and h-subepilocus (coding the hemimethylated sites: HMeCCG sites).

“MSAP_calc.r” was used to assess within-population diversity based on both AFLP and MSAP markers in terms of the number of polymorphic loci (Np), number of private loci (Npr), and Shannon’s information index H (Lewontin 1995). The frequency down-weighted marker values (Schönswetter and Tribsch 2005) were calculated using AFLPdat (Ehrich 2006). To facilitate interpretation, the values were range standardized to obtain the rarity index (RI) (Winkler et al. 2010).

Pairwise distances among individual plants were calculated using Dice’s coefficient. Mantel’s test was used to compute and test the linear correlation between the matrices based on AFLP and MSAP markers. Significance levels were obtained after 10,000 permutations using the NTSYS-pc ver. 2.2 (Rohlf 2008). PCoAs based on Dice’s distance matrices were performed using PAST version 2.01 (Hammer et al. 2001) to visualize the relationships between populations.

The analysis of molecular variance (AMOVA) (Excoffier et al. 1992) was performed using ARLEQUIN ver. 3.0 (Excoffier et al. 2005). AMOVA was used to partition the total genetic (based on AFLPs) and epigenetic (based on MSAPs) diversity among populations. The variance components were tested statistically by nonparametric randomization test using 10,000 permutations.

Results

Genetic and Epigenetic Structure of the Three Populations of L. bosniacum from Different Habitats

The AFLP analysis yielded a total of 663 polymorphic loci. The estimated error rate per primer combination ranged from 3.57% to 6.40% with a mean of 4.73%, which fell into the typical range of 2–5% usually found in AFLP data set (Bonin et al. 2004). Assessment of genetic diversity across the populations (table 3) revealed that the percentage of the polymorphic loci (%P) as well as Shannon's information index (H) was almost equal in all three populations. The percentage of private loci (%Npr) ranged from 0.60 (P02) to 1.96 (P03), while the RI ranged from −0.05 (P02) to 0.03 (P03). MSAP analysis resulted in 1,000 polymorphic loci that yielded 2,605 polymorphic subepiloci, consisting of 670 u-subepiloci, 969 h-subepiloci, and 966 m-subepiloci. By considering all MSAP subepiloci, both % P and H were comparable across populations with slightly lower values than those for AFLP markers. On the other hand, % Npr ranging from 5.91 (P02) to 8.41 (P03) was considerably higher on the epigenetic than on the genetic level. RI calculated from MSAP data showed a wider range of values across populations (from −0.11 for P01, to 0.11 for P03) when compared to that based on AFLP data.

Table 3

Genetic (AFLP) and Epigenetic (MSAP) Diversity within the Three Populations of Lilium bosniacum

Parameter	Population	AFLP	MSAP
			All	u-Subepiloci	h-Subepiloci	m-Subepiloci
Np
	All	663	2605	670	969	966
%P
	P01	88.08	71.82	57.61	80.60	72.88
	P02	88.08	72.13	60.45	70.28	82.09
	P03	88.54	74.24	65.07	79.26	75.57
%Npr
	P01	1.36	7.75	14.78	14.78	5.59
	P02	0.60	5.91	9.25	9.25	6.21
	P03	1.96	8.41	17.01	17.01	4.14
H
	P01	0.61	0.43	0.33	0.50	0.43
	P02	0.63	0.45	0.34	0.42	0.57
	P03	0.63	0.46	0.40	0.50	0.48
RI
	P01	0.02	−0.11	−0.06	0.01	−0.25
	P02	−0.05	−0.01	−0.44	−0.17	0.43
	P03	0.03	0.11	0.50	0.15	−0.19

Note.—Np, total number of polymorphic loci; %P, percentage of private loci; %Npr, percentage of privated loci; H, Shannon’s information index; RI, rarity index; u-subepilocus (coding of unmethylated sites); m-subepilocus (coding of HmeCG or MeCG sites); h-subepilocus (coding of HMeCCG sites); P01 = Kladanj (serpentine population), P02 = Jahorina (nonserpentine population), P03 = Crepoljsko (nonserpentine population).

Pairwise distances between individuals in a population, as calculated using Dice's coefficient based on MSAP data (average 0.773, ranging from 0.626 to 0.869), tended to be higher than those based on AFLP data (average 0.522, ranging from 0.379 to 0.636). A weak correlation was found between matrices based on AFLP and MSAP markers. Similarly, correlations between matrices based on u-, h- and m-subepiloci were weak (ru/m = 0.121; ru/h = 0.322) to moderate (rh/m = 0.474).

PCoAs based on genetic distances tended to separate the populations along the first axis but a high dispersion of the samples and overlap among populations along the second axis indicated high levels of intra-population diversity and absence of a strong population structure (fig. 3). On the contrary, PCoA of epigenetic distances (combined data) showed clear separation of all three populations along both axes, suggesting high levels of epigenetic divergence among populations (fig. 3). The separation of populations was slightly less clear when PCoAs were carried out based on distance matrices using solely u-, h-, or m-subepiloci (fig. 3).

. 3.

—PCoA of Dice’s distances among individuals belonging to three L. bosniacum populations (P01, P02, and P03) based on genetic (AFLP) and epigenetic (MSAP) data. (A) Genetic distances based on AFLP. (B) Epigenetic distances based on MSAP including all subepiloci. MSAP results for individual subepiloci (u-, h- and m-subepiloci) are shown in (C), (D) and (E), respectively.

One-way AMOVA based on genetic data showed that 4.95% of the total variation was attributable to the differences among populations (φST = 0.050; P < 0.0001; table 4). In concordance with the results of PCoA, AMOVA based on combined epigenetic data indicated much higher levels of the population differentiation when compared with genetic data—11.72% of the epigenetic variance actually resides among populations (φST = 0.117; P < 0.0001). Considering u-, h-, and m-subepiloci separately, the highest φST was observed for u-subepiloci (0.190) followed by m- (0.118) and h-subepiloci (0.073).

Table 4

AMOVA for the Partitioning of AFLP and MSAP Diversity among Lillium bosniacum Populations

Loci	Source of Variation	Variance Components	% Total Variation	φST	P
AFLP	Among populations	5.21	4.95	0.050	<0.0001
	Within populations	99.96	95.05
MSAP All	Among populations	36.20	11.72	0.117	<0.0001
	Within populations	272.60	88.28
MSAP u-subepiloci	Among populations	12.94	19.00	0.190	<0.0001
	Within populations	55.17	81.00
MSAP h-subepiloci	Among populations	8.40	7.32	0.073	<0.0001
	Within populations	106.30	92.68
MSAP m-subepiloci	Among populations	14.86	11.80	0.118	<0.0001
	Within populations	111.13	88.20

To gain more insight into the sources as well as robustness of the observed genetic and epigenetic differences between populations, we conducted additional pairwise comparison (table 5). Again, the results were in concordance with those of PCoA (fig. 3) and AMOVA for all three populations (table 4), giving similar φST for all pairs of AFLP (φST 0.028–0.074) and total MSAP loci (φST 0.114–0.123), supporting the observed much greater epigenetic than genetic differentiation. When MSAP results were broken down according to individual subepiloci, the findings from molecular variance partitioning were still supported, although the φST values were less homogeneous. The h-subepiloci uniformly showed the least amount of divergence, while the u- and m-subepiloci showed divergence consistent with the overall results, except for P1/P3 (m-subepiloci) and P1/P2 (u-subepiloci) population pairs, where the difference was less pronounced.

Table 5

Population Pairwise ϕST Values (below diagonal) and Their Significance (above diagonal) Based on AFLP and MSAP Markers

Loci	Population	Population
		P1	P2	P3
AFLP	P1		***	***
	P2	0.028		***
	P3	0.074	0.047
MSAP All	P1		***	***
	P2	0.114		***
	P3	0.115	0.123
MSAP u-subepiloci	P1		***	***
	P2	0.033		***
	P3	0.283	0.163
MSAP h-subepiloci	P1		***	***
	P2	0.074		***
	P3	0.047	0.099
MSAP m-subepiloci	P1		***	***
	P2	0.156		***
	P3	0.066	0.122

Note.—ns, nonsignificant value.

Significant at *P < 0.05, **P < 0.01, ***P < 0.001.

rDNA Karyotypes in the Three L. bosniacum Populations from Different Habitats

General chromosome features (number, size, and morphology) recorded in plants from the three L. bosniacum populations corresponded to the already described karyotype for L. bosniacum (Muratović et al. 2005). However, striking variations in number and position of 18S rDNA loci were revealed in the population from the microsite Kladanj (P01), which was growing on serpentine soil. While the standard karyotype of L. bosniacum has four 18S rDNA loci (located on chromosomes 2, 6, and 12) and one 5S rDNA locus located on chromosome 7 (karyotype 1, fig. 2), we found plants from serpentine population having five or six 18S rDNA loci: the new loci appeared on chromosome 7, collocated with the 5 S rDNA locus (karyotypes 2, 3, and 4; fig. 2, vertical arrows), and on chromosomes 1 and 3 (fig. 2, horizontal arrows). FISH performed on DNA fibers confirmed partial overlapping of 5S and 18S rDNA clusters (fig. 2). In the rDNA karyotype 4, the intercalary 18S rDNA locus disappeared from the long arm of the chromosome 12 and new loci appeared on chromosomes 1 and 3 (fig. 2). The rDNA karyotype 3 differs from karyotypes 2 and 4 by intensity of hybridization signals on chromosomes 2 and 12 (fig. 2, asterisks). On the other hand, the number and position of the 5S rDNA loci was the same in all plants analyzed from the three populations. The number of plants possessing specific rDNA karyotypes is given in table 2.

. 2.

—Karyotypes of serpentine L. bosniacum after simultaneous FISH (A–F) using 18S (Cy3-red) and 5S (FITC-green) rDNA probes, and partial karyotype after CMA staining (G). A–D: Idiograms showing standard (karyotype 1) and three “deviant” rDNA karyotypes B, C, and D (karyotypes 2, 3, and 4, respectively). Vertical arrows show the new 18S rDNA locus juxtaposed to 5S rDNA locus, and horizontal arrows show transfer of rRNA genes from the locus on the chromosome 12 to new chromosomal loci on the chromosomes 1 and 3. Asterisks indicate amplification of rRNA genes seen as more intense hybridization signal. (E) Metaphase plate corresponding to the rDNA karyotype 2, with one supplementary 18S locus on the chromosome 7 (arrow). (F) Karyotype with standard number/position of 18S and 5S rDNA loci (karyotype 1), but containing four B-chromosomes which bear two 5S rRNA loci (arrows); (G) Junction between CMA positive NORs (18S rDNA loci) of terminal and paracentromeric positions suggesting mobility of rRNA genes. (H) Double FISH on DNA fibers showing co-localization of the signals corresponding to 5S and 18S rRNA genes.

CMA banding revealed that distribution of GC-rich heterochromatin was the same in rDNA karyotypes 1, 2, and 3; nevertheless, there was a difference in the intensity of CMA bands at paracentromeric positions of chromosomes 2 and 12 in the rDNA karyotype 3, reflecting the difference in quantity of GC-rich heterochromatin (fig. 2, asterisks). In the rDNA karyotype 4, two additional CMA bands appeared at the position of the new 18S rDNA loci on chromosomes 1 and 3 (fig. 2, horizontal arrows).

In addition, one to four B chromosomes (Bs) were found in almost all plants from the population P01, bearing two 5S rDNA loci positioned at the chromosome termini (fig. 2, arrows).

Discussion

In the present study, we assessed genetic, epigenetic, and cytogenetic structure of three natural populations of L. bosniacum growing on different habitats, two of them associated with various types of environmental stress. We found that, despite of genetic similarity, the three populations showed clear epigenetic differentiation, probably affecting karyotype rearrangements observed in the P01 population growing on serpentine soil.

The AFLP and MSAP analyses showed little genetic (φST = 0.050), but clear epigenetic (φST = 0.117) differentiation of the three L. bosniacum populations growing on contrasting habitats. The correlation between the genetic and epigenetic distances (r = 0.318) between the three L. bosniacum populations was lower than reported in maize (r = 0.419) (Roy et al. 2015), rice (r = 0.520) (Choi et al. 2016), or potato (r = 0.780) (Cara et al. 2013), suggesting that only 10% (R2 = 0.101) of the epigenetic divergence between individuals could be explained by genetic distance. PCoA offered more insight by revealing that the three populations group into three distinct clusters based on the DNA methylation data, while failing to show a clear separation based on the genetic distances (fig. 3). Centers of clusters form an almost equilateral triangle, indicating similar epigenetic distance between all pairs of populations. Interestingly, when comparing genetic and epigenetic inter-individual variability (fig. 3) we found less epigenetic diversity and higher epigenetic than genetic inter-individual similarity. This is apparent when examining the scattering radius within each population (same color group) for AFLP and MSAP analyses (fig. 3), with clearly larger spread for AFLP than MSAP, even when MSAP data were broken down according to individual subepiloci. This finding also suggests that habitat conditions strongly affect epigenetic component of the genome in a convergent manner. Very similar situation was observed in potato natural populations and mangrove populations from contrasting environments (Marfil et al. 2009; Lira-Medeiros et al. 2010).

We assumed that the results of multivariate analyses (population diversity assessment; AMOVA, PCoA) have not been affected substantially by the inherent error rates of AFLP and MSAP techniques. In case of AFLPs, the estimated error rate fell into the expected range, while in the case of MSAPs no technical duplicates were included in the analysis (assuming that similar error rates could be expected) since we used essentially the same experimental procedure for MSAPs as for AFLPs. Analyzing Lavandula latifolia populations, others (Herrera and Bazaga 2016) reported that the average error rate of MSAPs (3.1%) was slightly higher than that of AFLPs (2.2%). However, in Viola elatior populations, the average error rate of MSAPs (2.1%) was lower than that of AFLPs (3.3%) (Schulz et al. 2014).

When the epigenetic distance was partitioned according to the type of epigenetic modifications (fig. 3), same clusters emerged as with the combined distance, suggesting that the contribution from different types of epigenetic modification was additive, that is, going in the same direction when separating populations. The MSAP scoring approach that we used in this work (Schulz et al. 2013) revealed that the relative contribution of the hemimethylated sites (h-subepiloci) to the overall population differentiation (φST) was much lower (φST = 0.073), and probably insignificant for population differentiation, than the contribution of other epigenetic states—the fully methylated, m-subepiloci (φST = 0.118), and the unmethylated, u-subepiloci (φST = 0.190). These results suggest that a functional difference in the methylation types of subepiloci contributes to the differentiation between populations. Most of the total variation was attributed to u-subepiloci indicating that demethylated state of the epiloci was associated with unfavorable environmental conditions present on the habitats. Similar results were shown in a number of studies, that is, biotic and abiotic stressors and environmental stress-induced demethylation and transcriptional activation of stress related genes (Wada et al. 2004; Choi and Sano 2007) or global genome demethylation (Lira-Medeiros et al. 2010; Wu et al. 2013; Schulz et al. 2014).

Number/position of rDNA loci is a species-specific karyotypic feature. Surprisingly, in the population P01, we found variation in rDNA karyotype regarding 18S rDNA loci. The three karyotypic events occurred in this population: 1) Emergence of the new 18S rDNA cluster on the chromosome 7 (karyotypes 2, 3, and 4) at the same site populated previously by 5S rRNA genes (karyotype 1); 2) amplification of rRNA genes within rDNA clusters on chromosomes 2 and 12 (karyotype 3); and 3) transfer of 18S rRNA genes from amplified rDNA clusters of chromosomes 2 and 12 (karyotype 3), as well as from intercalary 18S rDNA cluster of chromosome 12 (which completely disappeared in the karyotype 4), to the new paracentromeric locations on the chromosomes 1 and 3 (karyotype 4). We assume that the karyotypes 2 and 3 represent intermediate karyotypes from the karyotypes 1 and 4. These intermediate karyotypes possess a 18S rDNA locus on chromosome 7 and amplified 18S rRNA genes at loci on chromosomes 2 and 12, but still not the new 18S rDNA loci on chromosomes 1 and 3 (karyotype 4). It is important to emphasize that all four karyotypes are found within the same population; nevertheless, we do not know if the karyotype 4 is the “final” karyotype in this obvious karyotype diversification within the serpentine population of L. bosniacum.

The rDNA loci movements can be explained by several mechanisms including structural chromosome rearrangements (Rousselet et al. 2000), ectopic recombination between terminal rDNA loci (Pedrosa-Harand et al. 2006), and rDNA transposition mediated by TEs (Raskina et al. 2004; Datson and Murray 2006). Since 18S rDNA loci in the standard L. bosniacum karypotype have an intercalary or a paracentromeric position, the ectopic recombination between terminal rDNA sites is excluded as a possibility for rDNA movements. As no changes that would suggest extensive structural rearrangements were observed in chromosome size or morphology, the most plausible mechanism by which rDNA loci moved within L. bosniacum karyotype would be transposition mediated by TEs possibly reactivated due to stress-related global DNA demethylation that we observed in the P01 population (Slotkin and Martienssen 2007). The observed preferential insertion of TEs within rRNA genes had been described in several other species (Jakubczak et al. 1991; Raskina et al. 2004; Belyayev et al. 2005). A comparison of CMA banding and FISH patterns of the serpentine L. bosniacum suggests that the emergence of the new 18S rDNA loci on the chromosomes 1 and 3 in the rDNA karyotype 4 could be associated with amplification of GC-rich heterochromatin (observed in the karyotype 3). It is likely that amplified GC-rich heterochromatin, associated with the new 18S rDNA loci, contains some uncharacterized TE element. Reactivated by stress-related DNA demethylation, these mobile elements can be involved in transfer of rRNA genes to the new chromosomal locations in karyotype of the plants from the P01 population. The two ribosomal gene families were located on different chromosome pairs in standard karyotype of L. bosniacum. This is known as an “S-type arrangement,” which is the most frequent amongst angiosperms. In three deviant karyotypes (2, 3, and 4) found in the population P01, the 5S rRNA genes are linked to 18S rRNA genes, which represent the less common L-type arrangement (Garcia et al. 2012). Two different rDNA clusters overlap in the contact area, which was clearly visible on extended DNA fibers (fig. 2). To the best of our knowledge, this is the first report of coexistence of S- and L-type of rDNA arrangements within the same species, which is also suggestive of chromosomal rearrangements.

One to four Bs, possessing two terminal 5S rDNA loci, were found in all plants of the population P01. The previous study revealed one B-chromosome in karyotype of lily plants of atypical alpine population Bjelašnica (Muratović et al. 2005), however, containing no rRNA genes. Others (Belyayev et al. 2010) have demonstrated that increased numbers of Bs correlate with significant increases in the copy number of some TEs in small and marginal populations of Aegilops speltoides. Concentration of TEs around and inside ribosomal genes can increase fragility of rDNA sites (Raskina et al. 2008). One possible scenario for the genesis of Bs, bearing terminal 5S rDNA loci, in the PO1 population of L. bosniacum would be that they arose through A-chromosome breaks at rDNA sites associated with some unknown mobile elements. Indeed, we observed some junctions between NORs of terminal and paracentromeric positions, indicative of rRNA gene mobility (fig. 2). Some authors (Raskina et al. 2004, 2008) have stressed that micro-evolutionary genomic changes can be estimated indirectly from the mobility of rDNA clusters and that the dynamics of rDNA loci should be regarded as a strong indicator of intra-genomic changes.

One could object that we did not record any rDNA rearrangements in another population subjected to environmental stress (P02 population, Jahorina), even though this population was epigenetically differentiated as well, mostly by demethylated u-subepiloci. At this stage, we do not have any comprehensive explanation, although both populations were growing under unfavorable ecological conditions. Mechanisms by which TEs are involved in adaptation to local environments including high altitude, light availability etc., have been described in plants (for a review, see Casacuberta and Gonzalez 2013 and references therein). However, the links between TE activation and stress response are still by far more complex and exceed our present understanding of that process. First, TE reactivation depends on the type of mobile elements and the host genome (Lisch 2013; Chuong et al. 2017); second, specific TE sequences respond to different stressors; and third, TE reactivation can be mediated by other cellular mechanisms and not simply be a result of a direct response to a specific stressor (for a review, see Casacuberta and Gonzalez 2013). Here, we analyzed genome-wide methylation using isoschizomers MspI and HpaII, which limited our investigation to CpG dinucleotides. Even though substantial number of intragenic TEs in plant genomes are controlled by CG methylation (Dowen et al. 2012; Le et al. 2015), many of them are controlled by methylation of CpNpN and CpNpG motifs as well, which may respond differently to a particular stressor (Dowen et al. 2012).

According to the existing literature and our previously published data on L. bosniacum, it is possible to distinguish several basic morphological traits and phenology for lily plants growing under varying ecological conditions and substrates (table 6). Harsh environmental conditions, such as unusually high altitude and serpentine soil, have impact on morphology and phenology of Bosnian lilies—plants are of lower height and possess a smaller number of flowers per stem (1–4) than lily plants growing under habitual conditions (table 6, fig. 4). In addition, leaves of Bosnian lily plants growing on serpentine soil possess a lower density and larger stomata size (length and width) in comparison to leaves of plants from any other L. bosniacum population studied before (see table 6; and Muratović, Bogunić, et al. 2010). These traits are known to be adaptive for plants which regulate transpiration in conditions of water deficit (Brady et al. 2005; Pavlova 2009), which is the case of the Bosnian lily plants coping with osmotic stress associated with the serpentine substrate.

Table 6

Overview of Morphological Characters and Phenology of Bosnian Lily Populations Growing in Diverse Environmental Conditions

Environmental Conditions	Substrate	Altitude (m a.s.l.)	Bulbs Position in the Soil	Height (cm)	Flowering Time	Number of Flowers	Leaf Anatomy
							Stomatal Densitya	Length × Width of Stomata (µm)
Optimalb	Limestone (dolomite)	1,200–1,500	15–20 (25)	60–80	V–VI	3–10	13/6	72 × 50
Extremely high altitudesc	Limestone	1,800–2,000	<10 (15)	<30 (40)	VI–VII	1–2	12/6	80 × 49
Serpentinec	Serpentine	850–1,000	<10 (15)	30–50 (60)	IV–V	1–4	10/5	96 × 59

On abaxial epidermis/on adaxial epidermis of apical half of the leaf.

According to Matthews (1980), Muratović et al. (2005), Muratović, Bogunić, et al. (2010), Ritter-Studnička (1964).

According to Muratović et al. (2005), Muratović, Bogunić, et al. (2010), Ritter-Studnička (1968).

. 4.

—Bosnian lily plants from: (A) meadow population at site Crepoljsko (P03), and (B) woodland serpentine population at microsite Kladanj (P01). Two plants differ in their height and number of flowers.

In conclusion, our results suggest that epigenetic variation might be the first, immediate source of variability playing a role in stress response of L. bosniacum plants, which would provide the necessary time frame for genetic processes to take over and contribute to long-term stable heritability of changes reflecting adaptation to a new environment (Becker et al. 2011; Schmitz et al. 2011). Prominent rearrangements of 18S rDNA loci and presence of Bs bearing rRNA genes in L. bosniacum from the population P01 (microsite Kladanj) are suggestive of a micro-evolutionary process, probably governed by TE mobility as a result of global genome demethylation induced by stress present on specific habitat, i.e. serpentine soil (Raskina et al. 2004, 2008). In addition, some morphological traits and leaf anatomy are found on plants of this unique population growing on serpentine substrate. Therefore, this population certainly deserves further focus and more profound study in order to identify if the adaptive processes are going on within this population (Via 2009).