Introduction
Silene L. (Caryophyllaceae) is a key plant genus for studying crucial questions interrelating ecology and evolution (Bernasconi et al., Reference Bernasconi, Antonovics, Biere, Charlesworth, Delph, Filatov, Giraud, Hood, Marais, McCauley, Pannell, Shykoff, Vyskot, Wolfe and Widmer2009). Silene ciliata Pourret is a Mediterranean alpine plant that occurs in alpine pastures (Escudero et al., Reference Escudero, Gimenez-Benavides, Iriondo and Rubio2005) protected by the European Habitats Directive (Council Directive 92/43/EEC, 1992). The species is threatened by global warming (Giménez-Benavides et al., Reference Giménez-Benavides, Escudero and Iriondo2007, Reference Giménez-Benavides, Escudero, García-Camacho, García-Fernández, Iriondo, Lara-Romero and Morente-López2018) and has been included in catalogues of threatened species from different countries and regions (Dray, Reference Dray1985; Fernández et al., Reference Fernández, Prieto, Altuna, Arregui, Gutiérrez, Sánchez and Janices2007; Sanz et al., Reference Sanz, Carreño and Albacete2010; Légifrance, 2019). The seedling stage of this species experiences great mortality and it is, thus, subjected to strong selective pressure. This pressure may be qualitatively different between environmental conditions that are most commonly found at the species populations (optimal conditions) and those found only at the extreme of the species ecological range (marginal conditions). Previous studies have shown local adaptation patterns in optimal and marginal populations (Giménez-Benavides et al., Reference Giménez-Benavides, Escudero and Iriondo2007, Reference Giménez-Benavides, Escudero, García-Camacho, García-Fernández, Iriondo, Lara-Romero and Morente-López2018; García-Fernández et al., Reference García-Fernández, Escudero, Lara-Romero and Iriondo2015). Thus, S. ciliata is a key exponent to evaluate Mediterranean alpine species responses to oncoming global warming and corroborate evolutionary hypotheses (e.g. García-Fernández et al., Reference García-Fernández, Escudero, Lara-Romero and Iriondo2015; Kyrkou et al., Reference Kyrkou, Iriondo and García-Fernández2015; Lara-Romero et al., Reference Lara-Romero, García-Fernández, Robledo-Arnuncio, Roumet, Morente- López, López-Gil and Iriondo2016). Accordingly, there is great interest in developing genomic resources for the species, which would improve the understanding of the genetics of adaptation in Mediterranean alpine environments.
In this context, we carried out a transcriptomic study of seedlings of S. ciliata, with the following objectives: (1) to identify Single Nucleotide Polymorphisms (SNPs) in the transcriptome of the species, (2) to describe the biological function of the polymorphic genes expressed, (3) to examine diversity patterns of potential adaptive value in optimal and marginal populations of the species and identify loci that may be involved in local adaptation processes.
Experimental
We used RNeasy Plant Mini-Kit (QIAGEN) to extract and isolate RNA from six seedlings grown under controlled conditions, one for each of the six studied populations located in Central Spain. Three populations were located at the high edge and the other three at the low edge of the species elevational range (Table 1). The high and low-elevation edges represent optimum and marginal (warmer and drier) environmental conditions for the species, respectively. The quality of RNA was evaluated with a Qubit (Invitrogen, Carlsbad, CA, USA). One sequencing run was carried out in an Illumina platform through 100 bp paired-end reads. Trimming was carried out with software Trimmomatic (Bolger et al., Reference Bolger, Lohse and Usadel2014). Then, S. ciliata transcriptome was aligned with the genome of Silene latifolia (GenBank reference: GCA_900095335.1) using BWA software (Li and Durbin, Reference Li and Durbin2010).
Geographic and climatic characteristics of the six studied populations of Silene ciliata in Central Spain and genetic assessment of the sampled seedlings
For each plant, the population of origin, environmental classification, minimum annual temperature in (Min T), snowpack accumulation in thaw months (February, March and April), geographical coordinates, percentage of mapping with the reference genome, total number of sequences retained after trimming and estimators of genetic diversity are provided. Trimming configuration after the optimization process was: leading: 5, trailing: 5, sliding window: 4:15, minlen: 50. See trimmomatic manual (http://www.usadellab.org/cms/?page=trimmomatic) for further details on the selection of trimming steps and their associated parameters.
H O, observed heterozygosity; F i, coefficient of inbreeding.
a Nucleotide diversity per site ×10−3. Climatic variables were obtained from the ACPI (Atlas Climatic de la Peninsula Ibérica, http://opengis.uab.es/wms/iberia/index.htm). Monthly snowpack was calculated following the methodology proposed by López-Moreno and co-workers (Reference López-Moreno, Vicente-Serrano and Lanjeri2007) and it ranged from 0 to 1.
SNPs were identified using Reads2snp (Gayral et al., Reference Gayral, Melo-Ferreira, Glemin, Bierne, Carneiro, Nabholz, Lourenco, Alves, Ballenghien, Faivre, Belkhir, Cahais, Loire, Bernard and Galtier2013) and filtered with VCFtools 4.1 (Danecek et al., Reference Danecek, Auton, Abecasis, Albers, Banks, DePristo, Handsaker, Lunter, Marth, Sherry and McVean2011). Only, biallelic SNPs with no missing data and at least seven reads per genotype were retained to prevent the inclusion of false positive SNPs (Swarts et al., Reference Swarts, Li, Romero Navarro, An, Romay, Hearne, Acharya, Glaubitz and Mitchell2014; Marano et al., Reference Marano, Marcorin, Castelli and Mendes-Junior2017). Paralogous and singleton SNPs were further deleted. VCFtools was also used to estimate the genetic variation in the whole genome. Blastx software (Altschul et al., Reference Altschul, Gish, Miller, Myers and Lipman1990) and the database of SWISS-PROT (Bairoch and Apweiler, Reference Bairoch and Apweiler2000) were used to annotate the biological function (i.e. gene ontology terms) of the sequences carrying SNPs. We applied two different measures to detect candidate SNPs with unusually high-allele frequency differentiation between elevations. We first calculated allele frequency differences (AFDs) between low and high elevations at the individual allele level (Turner et al., Reference Turner, Bourne, Von Wettberg, Hu and Nuzhdin2010; Stölting et al., Reference Stölting, Paris, Meier, Heinze, Castiglione, Bartha and Lexer2015). SNPs were considered unusually divergent if AFDs were ≥3 SDs higher than the genome-wide average. Second, following Muller et al., Reference Muller, Latreille and Tollon2011, we computed the dispersion of each allele (m), which is the average elevational distance of an allele copy to the average elevation of all copies of that allele (β). Then, 1000 permutations of allele copies among all studied geographical locations were performed, to subsequently estimate m and β in each of the permutations. SNPs were considered unusually clustered if m were ≥2 SDs higher than the permutation average. We performed a GO term enrichment analyses using Fisher's exact tests to assess whether sequences containing selected SNPs were enriched in a biological function.
Discussion
Eighty percent of RNA sequences were conserved after trimming from an initial average of 30,000,698 ± 3,715,368 sequences. The percentage of sequences mapped against the reference genome ranged between 37.9 and 45 (Table 1). After filtering, we identified 147,118 SNPs distributed throughout 12,688 complete and partial sequences (SNPs per sequence: mean ± SD = 11.6 ± 14.51). Posterior probability per SNP was higher than 0.985 for all SNPs (mean ± SD = 0.9995 ± 0.0011). In total, 8023 polymorphic sequences were annotated. Their most common function was related to cellular processes, metabolic processes and biological regulation (Fig. 1, Table S1). Annotated sequences were deposited in the GenBank (BioProject ID: PRJNA528948). This extensive dataset provides a novel genomic resource for S. ciliata, and a significant step towards a better understanding of its genetics.
Most common biological functions found in polymorphic genes. The graph was created using the web Gene Ontology Annotation Plotting WEGO 2.0 (Ye et al., Reference Ye, Zhang, Cui, Liu, Wu, Cheng, Xu, Huang, Li, Zhou, Zhang, Bolund, Chen, Wang, Yang, Fang and Shi2018).
According to identified SNPs, individuals from high and low elevation presented similar values of genetic diversity (Table 1, all Wilcoxon rank tests: P ≥ 0.2). Fi was positive in five out of six plants (Table 1) and did not differ between elevations (Wilcoxon rank test: P = 0.7). Previous studies on S. ciliata using neutral markers also found similar estimates of genetic diversity across elevations (García-Fernández et al., Reference García-Fernández, Segarra-Moragues, Widmer, Escudero and Iriondo2012; Lara-Romero et al., Reference Lara-Romero, García-Fernández, Robledo-Arnuncio, Roumet, Morente- López, López-Gil and Iriondo2016). Overall, 775 sequences carrying SNPs associated with elevation were selected by at least one of the implemented approaches, but both shared only 130 of them. Almost 90% (n = 118) of them were successfully annotated (Table S1), but they were not significantly enriched in any biological process after FDR correction. However, about 15% of these sequences were associated with responses to stress (n = 19) and abiotic stimulus (n = 16) (Table S2). This is particularly interesting for the identification of loci related to adaptation to climate change (Giménez-Benavides et al., Reference Giménez-Benavides, Escudero and Iriondo2007, Reference Giménez-Benavides, Escudero, García-Camacho, García-Fernández, Iriondo, Lara-Romero and Morente-López2018), and improving the understanding of adaptation processes in the species.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S1479262119000157
Acknowledgements
Samples were sequenced in the ETH Zurich Genetic Diversity Centre (GDC). Niklaus Zemp from the GDC for providing access to the reference transcriptome of Silene latifolia. Alex Widmer and the Plant Ecological Genetics group (ETH Zurich) for their helpful recommendations and discussions throughout the study. This work was supported by projects AdAptA (CGL2012-33528) and EVA (CGL2016-77377-R) of the Spanish Ministry of Economy and Competitiveness (MINECO). CLR was supported by a Juan de la Cierva post-doctoral fellowship (MINECO: FJCI-2015-24712) and by a European Science Foundation ESF Exchange Grant (Reference No. 4794) within the ESF activity entitled ‘Conservation Genomics: Amalgamation of Conservation Genetics and Ecological and Evolutionary Genomics'. SSB was supported by a FPI pre-doctoral fellowship (MINECO: BES-2017-082317).
