Genetic structure of coast redwood (Sequoia sempervirens [D. Don] Endl.) populations in and outside of the natural distribution range based on nuclear and chloroplast microsatellite markers.

Coast redwood (Sequoia sempervirens) naturally growing in southern Oregon and northern California is one of the few conifer tree species that are polyploid. Despite its unique ecological and economic importance, its population genetic structure is still insufficiently studied. To obtain additional data on its population genetic structure we genotyped 317 samples collected from populations in California (data set C) and 144 trees growing in a provenance trial in France (data set F) using 12 nuclear (five random nuclear genomic nSSRs and seven expressed sequence tag EST-SSRs) and six chloroplast (cpSSRs) microsatellite or simple sequence repeat (SSR) markers, respectively. These data sets were also used as reference to infer the origin of 147 coast redwood trees growing in Germany (data set G). Coast redwood was introduced to Europe, including Germany as an ornamental species, decades ago. Due to its fast growth and high timber quality, it could be considered as a potential commercial timber species, especially in perspective to climate warming that makes more regions in Germany suitable for its growing. The well performing trees in colder Germany could be potential frost resistant genotypes, but their genetic properties and origin are mostly unknown. Within the natural range in southern Oregon and northern California, only two relatively weak clusters were identified, one northern and one southern, separated by the San Francisco Bay. High genetic diversity, but low differentiation was found based on the 12 nuclear SSR markers for all three data sets F, C and G. We found that investigated 147 German trees represented only 37 different genotypes. They showed genetic diversity at the level less than diversity observed within the natural range in the northern or southern cluster, but more similar to the diversity observed in the southern cluster. It was difficult to assign German trees to the original single native populations using the six cpSSR markers, but rather to either the northern or southern cluster. The high number of haplotypes found in the data sets based on six cpSSR markers and low genetic differentiation based on 12 nuclear SSRs found in this study helps us study and better understand population genetic structure of this complex polyploid tree and supports the selection of potential genotypes for German forestry.

Introduction

Although climate change can negatively affect the growth of certain native species in Central Europe [1], increased precipitation and temperature in winter expected in some regions in Germany [2] can make the environment more suitable for growing certain non-native tree species in these regions. Considering climate change and the forecast for future climate, German forestry would benefit from new adaptive management strategies [3]. To secure German timber production and the growing demand for wood products in Europe, more non-native tree species should be tested for potential introduction. One of the successful examples is the North American species Douglas-fir (Pseudotsuga menziesii [Mirb.] Franco), which is widely established in European forests [4,5]. The predicted altitudinal shift in the distribution of woody species [1,6] due to climate change would make Germany more suitable for species adapted to the Mediterranean climate, such as coast redwood that has already been introduced to Europe as an exotic species [7].

Coast redwood (Sequoia sempervirens [D. Don] Endl.) is one of the four redwood species in the Cupressaceae family [8]. All four are characterized by a particular red coloured wood, an endemic narrow natural distribution range, and are listed as threatened or endangered by the International Union for Conservation of Nature (IUCN) [8,9]. Coast redwood stands out from the other three redwood species due to its valuable wood, high growth rate and natural vegetative reproduction [10-12]. In the natural range, the western USA, it is a very important timber species growing along the Pacific Coast in southern Oregon and northern California [10]. Moreover, coast redwood was introduced into several countries with similar climatic conditions for timber production, such as France [13], New Zealand [14] and China [15].

The coast redwood natural distribution range includes the three USDA plant hardiness zones 9b, 10a and 10b (minimum temperature from -4°C to 5°C) categorized as humid Mediterranean climate [16]. Therefore, coast redwood usually experiences light frost only few weeks during a year [8] and is frost sensitive, especially at young age [10]. Its needles are able to uptake water directly from the atmosphere [17], and the tree highly depends on frequent fogs that compensate the low precipitation during summer months [10,18]. Plant hardiness zones in Germany range from 7a (-17°C to -15°C) to 9a (-6°C to -4°C) [19] and currently present a high risk for planting coast redwood for timber production in these regions. The hardiness zones 8b (-12°C to -9°C) and 8a (-9°C to -6°C) along the Rhine valley in western Germany have climate conditions most similar to the natural distribution range of coast redwood, therefore most of the planted trees can be found there. To mitigate losses in German and European forests due to climate change, it is suggested to plant adapted genotypes of major timber species, native and introduced, which are able to cope with extreme weather conditions [13]. In the case of coast redwood some trees planted in the Rhine valley exhibit potential frost resistance and survived temperatures as low as -20°C [7], but their origins are unknown. Information about their origins with the particular climatic conditions would help to select additional suitable genotypes to increase the genetic diversity of German frost resistant stands. This would improve the adaptive potential and therefore decreases risks of coast redwood timber production in Germany considering future climatic uncertainties. To achieve this goal, the current genetic diversity of German stands need to be assessed and compared to diversity of populations within the natural distribution range. Identifying distinct populations within the natural range by analysing the population structure with genetic markers will be the basis to assess their different adaptive potential. Additionally, these markers could be very efficient to assign introduced plants to their origin [4,20,21], but the detection power depends much on population differentiation [21,22]. Until now, previous coast redwood studies were not able to accomplish that based on the available allozyme and microsatellite markers [23-31].

The earlier published results based on samples from the natural distribution range genotyped using one chloroplast [30] and six nuclear [31] microsatellite markers suggested high genetic diversity, but low population differentiation with a weak indication of two clusters. The weak population structure is not surprising for a long-living, wind-pollinated tree species with relatively small area and no physical barriers for gene flow. We anticipated that additional seven expressed sequence tag SSR (EST-SSR) and six chloroplast SSR (cpSSR) markers combined with the already available five nSSRs can enhance the population structure resolution power and help us identify divergent populations. This is the first study of coast redwood where a combination of EST-SSR, nSSR and cpSSR markers was used. In contrast to the random nuclear nSSRs, EST-SSR markers could be potentially under selection [32] and, although cpSSRs are more conservative and have a lower mutation rate than random nuclear nSSRs, they could be also very polymorphic and informative [33], especially for tracing long distance migration because of their strict paternal inheritance via pollen in coast redwood [34]. In our study we used already published five microsatellite or nuclear simple sequence repeat (nSSR) markers [29,31,35] together with newly developed six cpSSR and seven EST-SSR [36] markers. This is also the first study presenting data on SSR markers genotyped in both Californian and German coast redwood populations.

We expected that this new marker combination would improve the resolution of population structure within the natural distribution range and help us discriminate populations. In case of highly differentiated population structure we should be able to assign German frost resistant trees to their origin. In addition to the samples representing the complete natural distribution range, we collected samples from sites with suboptimal growing conditions in California and included them in our analysis hoping that they could help us identify the origin of German trees that are supposedly adapted to drier and colder climate and could likely originate from sites with colder and drier environment. The obtained data can help to select the most appropriate genotypes for in-situ conservation management strategies and timber production in and outside the natural distribution range.

The main objectives of this study were to 1) compare genetic population structure resolution based on nuclear (EST-SSRs together with nSSRs) vs. chloroplast (cpSSRs) microsatellite markers, 2) study genetic diversity based on EST-SSRs, nSSRs and cpSSRs within and outside the natural distribution range, and 3) assign trees growing in Germany to the original populations.

Materials and methods

Plant material

Three data sets were used in this study: French (F), Californian (C) and German (G). In 1984 an international provenance test was established by Kuser et al. [37] (so-called “Kuser provenance test”) using seedlings representing 90 provenances in the natural distribution range of coast redwood (S1 Table). The 180 individuals from this test were propagated by cuttings and used to establish hedge orchards in California, Oregon, Spain, Britain, France and New Zealand, including a planting site established in 1990 in St. Fargeau, Central France. For data set F, depending on accessibility, either needle or cambium (when braches with needles were unreachable) material was collected in 2014 from 153 surviving trees in St. Fargeau and frozen until DNA extraction. The data set F was compared with another “Kuser provenance test” experiment established in 1984 at a site close to Berkeley, California, named the Russell Reserve. Douhovnikoff and Dodd [31] partitioned 135 Russell Reserve trees into 17 watersheds according to their original geographic location and the GPS data (Fig 1) for population genetic structure analysis based on six nSSR markers [31]. To compare our data with their results and population based analyses [31], the St. Fargeau samples in the data set F were also partitioned into the same 17 watersheds also based on their original geographic location (S2 Table). Original locations of the ‘Kuser’ samples placed on the map with the mean monthly temperature pattern for the time period 1979–2013 are presented in S1 Fig.

Fig 1

Map of 17 watersheds (A-Q) along a latitudinal range from 42° 12' to 35° 55' N within the natural distribution range of Sequoia sempervirens for the French data set F (Fig 1 in Douhovnikoff and Dodd [31] reproduced with permission from The American Midland Naturalist).

The Californian data set C consisted of 309 tree samples collected from 16 locations within the natural distribution range in central and northern California in 2017. Number of samples, altitude, historic precipitation and mean January temperature per location are listed in S3 Table. The locations were chosen according to their geographic position and climatic conditions. The optimal climatic conditions for coast redwood are along the coast, with high oceanic influence and humidity and high ground water level [38]. The locations selected for sampling were at higher altitude, more interior, with colder and less foggy conditions, and therefore represented drier and more extreme habitats [18,39]. The original locations of 12 frost resistant trees in the St. Fargeau provenance site were also considered while selecting sources for the potentially frost resistant trees. These 12 frost resistant trees were identified by the owner in 1992, after a late frost occurrence in 1991 (unpublished data). Among the selected sites it was possible to visit 16 locations (Fig 2), and trees with a minimum diameter at breast height (DBH) of 15 cm were randomly sampled at each location. It should be noticed that there is no information proving that the sampled California stands represented natural populations. Some or many of them could be at least managed, if not even planted. The collected needle material was dried and stored in silica gel until DNA extraction.

Fig 2

Collection sites for the Californian data set C depicted by stars in the California map generated using the SimpleMappr online software (https://www.simplemappr.net).

The German data set G included in total 147 samples representing six sites: “Sequoiafarm Kaldenkirchen”, “Aboretum Burgholz”, “Weltwald Bad Grund” and three German Botanic Gardens in Bayreuth, Chemnitz and Göttingen, respectively (S2 Fig and S4 Table). The “Sequoiafarm Kaldenkirchen” was established by the Martin family in the 1950’s using coast redwood seedlings from probably one single tree in the Californian “Schenck-Grove”, Prairie Creek Redwoods State Park (Humboldt County, California, USA). Cuttings from these “best performing” clones were planted all over Germany and are known as the “Martin-Clones” (M. Geller, pers. communication). More than 100 individuals were planted in the “Aboretum Burgholz” representing the “Martin-Clones” and other sources in unknown numbers and combinations. The “Sequoiafarm Kaldenkirchen” and “Aboretum Burgholz” are both located in the Rhine valley in the plant hardiness zones 8a and 8b [19], respectively. The “Weltwald Bad Grund” is an arboretum located in the western Harz Mountains near Bad Grund with relatively cold winters and dry summers, typical for the plant hardiness zone 7b. Coast redwood cuttings and possibly seedlings from trees of the “Aboretum Burgholz”were planted in the Weltwald, but also in unknown number and combinations. Coast redwood trees collected from botanic gardens in Chemnitz, Bayreuth and Göttingen represented few individuals with unknown origins that survived in the coldest areas in Germany. Chemnitz and Bayreuth are in the plant hardiness zone 7a, and Göttingen in zone 7b. Fresh needle material was collected from all these individual trees in Germany and frozen at -20°C until DNA extraction.

For allele scoring validation, needles from 30 approximately three-year-old clones, including 2–3 ramets per clone, respectively, were collected in the Allerweltsgrün nursery (Köln) and frozen at -20°C (S3 Fig).

DNA isolation and PCR amplification of the microsatellite (SSR) markers

DNA was isolated from needles or cambium using the DNeasy Plant Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. The isolated DNA was diluted in ddH2O 1:10 for PCR amplification and stored at -20°C.

Five nSSR, seven EST-SSR and six chloroplast (cpSSR) microsatellite markers were used in this study to genotype samples in all three data sets F, C and G, respectively (S5 Table). The same touch-down PCR program was used for all 18 PCR primer pairs following the protocol described in Breidenbach et al. [36]. The PCR products were separated and visualized using the ABI Genetic Analyser 3130xl with GENSCAN ROX 500 as an internal size standard.

Verification of the nuclear microsatellite markers

The PCR primer nucleotide sequences for the 12 nuclear SSR markers were mapped to the coast redwood draft nuclear and chloroplast genome assemblies (made recently publicly available) to verify annealing sites for the microsatellite markers used in this study. The primer sequences were mapped against scaffolds downloaded from the NCBI GenBank (accession number VDFB00000000) using the CLC Genomic Workbench v11.0.1 software (Qiagen, Hilden, Germany) with the following parameters: complete match for the last 15 nucleotides at 3’ end, maximum mismatch for 3 nucleotides per annealing site in total and a maximum of 400 nucleotides between annealing sites.

Genotyping of the nuclear microsatellite markers

GeneMapper 4.1. (Thermo Fisher Scientific, USA) was used for visualization and fragment size calling of the PCR products. Coast redwood is a hexaploid species [40], which complicates microsatellite scoring [35,41]. In this study, we used the genotyping method of Pfeiffer et al. [42] to identify consistent and reproducible alleles. Following this method, all fluorescent peaks detected in 2–3 ramets per each of 30 different clones and representing different PCR amplified fragments were assigned in each sample to one of the 10 rank categories according to the shape and intensity of the peak and the motif of the primer. To determine which ranking has the highest probability, input files with different ranking combinations from category 7 to 10 with all 12 primers were generated. For example, for the ranking combination 7, all peaks of each sample assigned to the category 7 and above (8, 9 and 10) were included in the analysis, for the ranking combination 8, all peaks of each sample assigned to the category 8 and above (9 and 10) were included into the analysis. The R-package “polysat” [43,44] specifically developed to study genetic data of polyploid species was used for each ranking combination data set to calculate the pairwise ‘bruvo genetic distance’ between individuals [45]. The supposedly compound RW56 marker, in fact, demonstrated alleles which sizes followed regular tetranucleotide repeat differences in our study. Therefore, it was treated as a regular SSR marker with a tetranucleotide motif.

Using the R-package “ape” [46] and the pairwise distance matrix, a neighbour-joining (NJ) consensus tree was generated based on 1000 bootstraps. Since each marker differed in the allele pattern and peak shape, 25 different ranking combinations for the markers were tested. The ranking combination with the highest bootstrap value and correct grouping of the ramets representing the same clone was used for all further genotype scoring. The NJ results were confirmed with the function “assignClones” of the R-package “polysat” [43] using the miss-matching threshold of 0.2. The markers final ranking categories were from 7 to 10 (S5 Table). The NJ tree (NJT) based on the final ranking combination of the markers and 1000 bootstraps is presented in S3 Fig. In addition, for the data sets F and C, scored genotypes were converted also into a presence-absence matrix, and the Nei’s genetic distances ([47] after [48]) were calculated between watersheds F and populations C using the AFLP-SURV software with 1000 permutations [49]. The genetic distance matrix was used to generate a NJT with the PHYLIP v.3.69 software [50], respectively. The consensus tree was visualized using the FIGTREE software [51]. NJTs, using cpSSRs, were based on Nei’s genetic distance ([47] after [48]) and 1000 bootstraps for the reference data sets C and F and generated using the R-packages “adegenet” [52,53] and “poppr” [54,55]. To facilitate comparison between the phylogenetic trees based on the watersheds in data set F and Douhovnikoff and Dodd [31], watershed I (Mendocino County) was also used as a root of the NJT.

To compare genetic diversity levels within and outside the natural distribution range, the data sets F, C and G were also combined into one set and partitioned into the following three subgroups: NORTH and SOUTH representing the watersheds and locations located either North or South of the San Francisco Bay, respectively, and GER representing the German trees in the data set G. The R-package “polysat” was also used for the calculation of number of alleles, Shannon index as a parameter of genetic diversity and number of private alleles for each population in C, watershed in F, and group, respectively. Further, pairwise Jost’s D and FST parameters of genetic differentiation were calculated and bootstrapped 1000 times to identify significant values using the R-package “polysat”. Additionally, analysis of molecular variance (AMOVA) was conducted for data sets C, F and G with 999 permutations, respectively using the converted presence-absence matrix and the software GenAlEx [56].

STRUCTURE analysis

Based either on 12 nuclear SSR or six cpSSR markers, input files with the original allele sizes for the data sets F and C were used by the STRUCTURE v2.3.4. program [57-60] to infer the most likely potential number of distinct genetic clusters (K) by testing for different number of K from 1 to 24 with 20 iterations per run using the MCMC with 10 000 burn in’ and 100 000 final iterations assuming the admixture model [61].

The STRUCTURE HARVESTER [62] and the ClumPPAK [63] programs were used to visualize the STRUCTURE results and to help determining the most likely number of clusters using the ΔK approach. The STRUCTURE analyses were repeated also with the “locprior” optional parameter engaged (which is used to assist the clustering by specifying populations a priori), because low population differentiation was expected. Following the STRUCTURE HARVESTER results, the data sets F and C were divided each into two separately analysed groups—12 northern populations and four southern populations located to the north or south of the San Francisco Bay, respectively, to test for further clustering within these two groups using the same STRUCTURE settings.

In addition, the German sample (G) together with French (F) and Californian (C) samples was also analyzed with STRUCTURE for K = 17 (number of watersheds) and K = 30 (total number of subpopulation samples or groups in the F and C data sets) by using USEPOP flag (known origin) for F and C samples, but not for G (unknown origin) samples to assign individual German trees to either a particular watershed or a particular sample, respectively. The Q-values in two Q-matrices obtained for K = 17 and K = 30, respectively, were used as genetic traits to generate two pairwise genetic distance matrices using geometric index of genetic similarity developed by Zhivotovsky [64]. Briefly, using Q-matrix , where n is a total number of trees, K–number of clusters used by STRUCTURE to generate the Q-matrix, and each tree is presented by Q-values, so that their sum in each row equals 1: P+P++P=1 for any, i (i=1,2,…n), then, index of similarity (r) between trees i and j is calculated as . This index is always positive (non-negative) and does not exceed 1. The advantage of this index is that it takes into account even small Q-values. When vectors (trees in our case) are completely identical, then, it is equal to 1, and when the vectors are completely orthogonal, it is 0. The distance is simply calculated as 1-r.

Then, the obtained distance matrices were used for cluster analysis by generating NJTs with the R-package “ape” [46] and Principle Coordinate Analysis (PCoA) plots with GenAlEx [56].

Genotyping of the chloroplast microsatellite (cpSSR) markers and haplotype network

Due to the haploid nature of chloroplasts, genotyping of the cpSSR markers was easier and unambiguous. Based on all three data sets, a haplotype network was built using the Goldstein distance [65] in the program EDENetwork v2.28 [66]. The program calculates a weighted network based on the Goldstein distance between haplotypes with an automatically calculated percolation threshold of 2.67 [67].

Genetic assignment using the chloroplast microsatellite (cpSSR) markers

The individuals from the German data set G were assigned to the watersheds and populations in data sets F and C using the cpSSR markers and the GeneClass2 software [68]. For data set C, eight closely located sampling sites (< 1.0 km) were combined into four: PUC and LPF were combined into LPF (watershed M), MtMa1 and MtMa2 into MtMa (watershed O), HW1 and HW2 into HW (watershed G), and JDF and CCJDF into JDF (watershed J). A quality assessment (QI) was based on the self-assignment test following Hintsteiner et al. [4] using the computation criteria of Rannala and Mountain [69], the simulation algorithm of Paetkau et al [70] with 10 000 resampled individuals, and Type 1 error of 0.01 [69,71].

Results

Mapping of the PCR primer nucleotide sequences for the SSR markers to the draft coast redwood nuclear and chloroplast genome assemblies found annealing sites facing each other in a correct configuration allowing amplification for all six chloroplast markers and eight out of 12 nuclear markers (S7 Table). Three different single scaffolds included annealing sites for three markers ss36782, ss73361 and ss114481, respectively, two different scaffolds included annealing sites for each of two markers RW48 and ss73307, respectively, three, four and six different scaffolds included annealing sites for ss73978, ss91170 and ss74800 markers, respectively. The chloroplast markers all had only a single annealing site pair per each marker in the chloroplast genome, but multiple scaffolds also contained annealing sites. The latter could be due to either incomplete removal of chloroplast sequences from the nuclear genome assembly or presence of chloroplast sequences in the nuclear genome as a result of misassembling or migration of chloroplast genes to the nuclear genome [72-74].

The NJTs for the data sets F and C had low bootstrap values below 60% indicating weak differentiation between watersheds and populations based on 12 nuclear SSR markers (Figs 3 and 4), except for two watershed pairs P-Q and E-G in the data set F (Fig 3) and populations AET and CCJDF in the data set C (Fig 4). The NJTs based on the cpSSR markers showed higher bootstrap values for most clusters, but they did not reflect their geographic relationship (S4 and S5 Figs).

Fig 3

The neighbour-joining tree (NJT) for the French data set F partitioned into 17 watersheds (A-Q) following Douhovnikoff and Dodd [31] based on Nei’s genetic distance ([47] after [48]) calculated using 12 nuclear SSR markers.

Watersheds from A to Q are distributed from north to south in central California (see also Fig 1 and S1 Table). Bootstrap values are presented as percentage.

Fig 4

The neighbour-joining tree (NJT) for populations in the Californian data set C based on Nei’s genetic distance ([47] after [48]) calculated using 12 nuclear SSR markers.

County names of the sampled locations are in brackets. Bootstrap values are presented as percentage.

Highest genetic diversity values were found in northern and southern populations HW2 and MtMa1 belonging to data set C, and in northern watersheds G and L belonging to data set F (Table 1). The highest Shannon index (5.58) was found in the group of northern populations (NORTH) and the lowest (3.76) in the group of German samples (GER). Populations and watersheds with high genetic diversity contained also high number of alleles (Table 1). Highest numbers of private alleles were found in populations AET and RERI and watersheds J and N. AMOVA results revealed high variation within populations for all three data sets C, F and G. Highest variation among populations was found for data set C (Fig 5). Pairwise Jost's D- and FST-values were low and insignificant (S6 Table). Within data set C, population AET was the most and HW2 the least differentiated in comparison to all other populations, except RERI and MtMa1 (S6A Table). Within data set F the watersheds O, P and Q located to the south of the San Francisco Bay were differentiated from the watersheds B, C, D and E located to the north of the San Francisco Bay (S6B Table).

Table 1

Diversity measures for populations in Californian data set (C), 16 watersheds (A-Q) in French data set (F) and three groups consisting of the trees pooled from populations located above (NORTH) and below (SOUTH) San Francisco Bay in both C and F sets, and all German genotypes together (GER) based on 12 nuclear SSR markers.

Watershed	Population	Data set	N	Total number of alleles	Shannon Index	Number of private alleles
North
A		F	8	83	2.08	3
B		F	8	79	2.08	4
C		F	7	70	1.95	2
D		F	5	62	1.61	1
E		F	5	62	1.61	1
F		F	8	80	2.08	1
G		F	13	95	2.56	1
G	HW1	C	29	158	3.14	4
G	HW2	C	47	189	3.78	5
G	WEO	C	11	117	2.16	0
I		F	7	79	1.95	1
J		F	8	85	2.08	6
J	CCJDF	C	16	128	2.52	3
J	JDF	C	13	122	2.30	3
H		F	6	64	1.79	1
K		F	8	81	2.08	3
L		F	9	86	2.20	5
M		F	5	67	1.61	3
M	AET	C	6	65	1.39	5
M	ANG	C	6	83	1.61	0
M	DL	C	15	117	2.56	3
M	LPF	C	12	126	2.48	2
M	PUC	C	26	151	3.03	5
M	EN	C	17	138	2.77	2
M	TC	C	17	136	2.77	1
South
N		F	6	77	1.79	6
N	RERE	C	11	109	2.20	3
O		F	7	79	1.95	5
O	RERI	C	26	158	3.14	6
O	MTMA1	C	35	149	3.21	3
O	MTMA2	C	22	148	3.00	4
P		F	4	56	1.39	3
Q		F	4	62	1.39	1
Groups
	NORTH	C & F	267	240	5.58	60
	SOUTH	C & F	86	192	4.45	21
	GER	G	44	152	3.76	18

Fig 5

AMOVA results based on 999 permutations for data sets C and F and groups NORTH, SOUTH and GER (G) using the 12 nuclear SSR marker genotypes transformed into binary data.

The STRUCTURE analysis based on the six cpSSR markers suggested very low differentiation with maximum two clusters (K = 2) for the data sets C and F (Figs 6 and 7). The “locprior” function did not affect much the STRUCTURE results, therefore only results obtained with this function engaged are presented here. Additional STRUCTURE analysis performed separately within the northern and southern subpopulations did not find additional clusters and confirmed that the differentiation observed in the data sets F and C was mainly due to the differences between northern and southern populations. Results based on the 12 nuclear SSR markers suggested K = 2 and 6 with the “locprior” function engaged and disengaged for data set C, respectively, and K = 2 for both “locprior” function options for data set F (S7 Fig). However, the clustering was not strongly expressed (S8 Fig).

Fig 6

STRUCTURE analysis demonstrating the most likely number of clusters (K) using the “locprior” function for 16 populations in the Californian data set C based on 6 cpSSR markers and two separately analysed groups; 12 northern populations (K = 3) and four southern populations located to the north or south of the San Francisco Bay, respectively (K = 4).

Fig 7

STRUCTURE analysis demonstrating the most likely number of clusters (K) using the “locprior” function for the French data set F based on 6 cpSSR markers.

Using “assignClones” of the R-package “polysat” [43] 95 individual trees from the Burgholz stand in the data set G were identified as likely clones of the same SF72 tree from the “Sequoiafarm Kaldenkirchen”. Other individuals from Burgholz and Bad Grund were very similar to the SF76 tree, but with insufficient statistical support to assign them all unambiguously to one clone. Each clone was represented by only one individual for the diversity estimate of German samples (Table 1).

Based on six cpSSR markers, 63 haplotypes were found in the data set F, 109 in C and 22 in G. In total, all three data sets together contained 150 different haplotypes (Fig 8). The haplotype network and its percolated clusters generated by the EDENetwork software indicated no geographic structure in relationships between the haplotypes (Fig 8). More than half of all haplotypes were unique—89 haplotypes occurred only once within the three data sets. The haplotypes GG and KK were the most frequent and occurred with 17% and 14% frequency, respectively.

Fig 8

Haplotype network based on six cpSSR markers genotyped in three data sets F, C and G.

Size of the nodes reflects number of individuals assigned to the respective haplotype. The percolation threshold is 2.67 Red lines present the closest links, blue to green, yellow and no colour lines indicate links following decreasing relationship between haplotypes.

Individuals with identical genotype representing the same clone were analysed as a single entry (S4 and S8 Tables). Accuracy of the assignments was evaluated by the quality index (QI), which was low for both F (QI = 7%) and C (QI = 16%). Most of the assignments of German trees to the southern or northern populations divided by the San Francisco Bay in the reference data set C with scores above 48% were in agreement with the assignments in the reference data set F (Tables 1 and S8). Samples with scores below 48% were assigned inconsistently between the two references, to geographically very different watersheds. For example, B11 and SF77 were assigned to RERI south of the San Francisco Bay (watershed O) based on the reference data set C, but to Del Norte close to the Oregon border (watershed B) based on reference F (Fig 1 and S8 Table). The individuals from the Sequoiafarm, which were half- or full siblings according to the owner, were assigned to various watersheds (F) and locations (C).

We used also pairwise individual tree genetic distances between all 396 trees in German (G), French (F) and Californian (C) data sets based on similarity in composition of their Q-values for cluster analysis to assign individual German trees to a particular watershed or a subpopulation by generating NJTs and PCoA plots, which are presented in S9–S12 Figs, respectively. The individual Q-values for K = 17 (number of watersheds) and K = 30 (total number of subpopulation samples or groups except German sample) and genetic distance based on them for all 396 trees in German (G), French (F) and Californian (C) data sets are presented in S9 Table. Q-plots based on the STRUCTURE results for K = 17 and K = 30 for all 396 trees in German (G), French (F) and Californian (C) data sets are presented in S13 Fig. The German trees closely clustered together with genetically similar trees in French and Californian data sets F and C representing different watersheds (S9 and S11 Figs) and groups (S10 and S12 Figs) are presented in Table 2.

Table 2

German trees in the data set G clustered together with genetically similar trees representing different watersheds (S9 and S11 Figs) and groups (S10 and S12 Figs) in French and Californian data sets F and C, respectively, based on the pairwise similarity in composition of Q-values obtained by STRUCTURE for K = 17 and K = 30.

K = 17		K = 30
German trees	F and C	German trees	F and C
SF75	N-G_HW_111	SF75	N-G_HW_226
BA030549	S-N_F76	BA030549	N-G_HW_111, N-G_HW_207
CH1	N-H_F136	CH1	N-H_F136
SF76	N-H_F14, S-O_RERI_11	SF76	N-J_CCJDF_09, N-J_CCJDF_15
GOEB	S-O_RERI_18	GOEB	S-O_RERI_18
GOEP	S-O_MTMA_1B5	GOEP	S-Q_F93 | S-O_MTMA_208, S-N_RERE_01
SF3, SF86	N-F_F171, N-F_F174	SF86, SF88	N-F_F171 | N-F_F172, N-D_F25
SF88	N-G_F32, N-K_F40
SF63 | SF64 | SF69, SF70 | SF91, SF93	N-K_F129	SF64, SF93 | SF3 | SF63, SF91 | SF69, SF70 | SF82, SF90	N-I_F51 | N-K_F130, N-M_F44 | N-G_F29, N-K_F129
SF80 | SF82, SF90	N-F_F178	SF80	N-F_F178
SF67, B42	S-N_F52	B42	S-N_F52
SF74	S-Q_F83	SF67, SF74	N-M_EN_15, N-M_EN_17 | N-I_F63 | N-M_AET_ANGST03, N-M_DL_D15
SF66, SF77	N-I_F37, N-G_F35 | N-F_F167	SF77	N-C_F17, N-G_HW_101 | N-G_HW_218, N-M_AET_ANGST06
		SF66	N-B_F11, N-G_F35
SF81 | SF79, BG120	N-B_F70 | N-J_F38, N-G_F160	SF81SF79, BG120	N-B_F70, S-Q_F83 N-J_F38
B30	S-O_MTMA_1A8	B30	N-G_HW_219, N-M_PUC_301
SF71	N-M_PUC_303	SF71	N-G_HW_227
SF73	N-M_TC_01	SF73	N-B_F10
B4	N-M_EN_08	B4	N-L_F57
SF72	N-M_LPF_09	SF72	N-M_EN_08
BG121	N-H_F137	BG121	N-M_DL_D04 | N-M_LPF_04 | N-M_DL_D02, N-M_DL_D03 | S-O_RERI_17 | N-G_HW_230, N-M_DL_D05
GOEK	N-G_HW_223
SF84	N-G_HW_245	SF84	N-G_HW_245
BA10067	S-O_RERI_22	BA10067	S-O_MTMA_1C8
BG124 | BG123 | BG125, BG128 | B102 | BG127 | BG122, BG126	S-O_MTMA_1C8, S-O_MTMA_213 | S-O_MTMA_217 | S-O_MTMA_214	B102, BG127 | BG122 | BG126, GOEK | BG124 | BG128 | BG123, BG125	N-M_EN_13 | N-G_WEO_HRED | S-O_MTMA_205, S-O_MTMA_210

The first letter N or S in the tree ID name in the F and C data sets means North or South location, respectively. The second letter (A-Q) means watershed (see also Table 1). The trees that have F preceding the number in the name belongs to the F data set (for instance, N-H_F136 or S-Q_F83), all other trees belong to the C data (for instance, N-G_HW_111 or S-O_RERI_18). Highlighted by bold are the same area part (N or S), watershed or group in the tree names for the tress that were associated for both K.

Discussion

The presented study is the first study investigating coast redwood populations within and outside the natural distribution range with a high number of nuclear SSR markers and chloroplast SSRs. The new and significantly expanded marker set confirmed results of previous studies [27,28,31]. Surprisingly, neither the novel nuclear nor chloroplast microsatellite markers could improve resolution of population differentiation within the natural distribution range of coast redwood despite high polymorphism of the used microsatellite markers. It is interesting that supposedly frost resistant German trees were assigned to both of the two clusters found within the natural range, which shows a high genetic diversity and adaptive potential in the investigated German trees.

The NJT for the data set F based on 12 nuclear SSR markers in our study was in consensus with the NJT calculated for the same 17 watersheds (A-Q) by Douhovnikoff and Dodd [31] using clones from the Russell Reserve. The discrepancies between the two trees could be explained by very low bootstrap support for most clusters in both trees and by using the pairwise FST values for clustering in Douhovnikoff and Dodd [31] instead of the Nei’s genetic distance used in our study. We did not use the regular pairwise FST for generating NJT because their values could be greatly underestimated for highly polymorphic microsatellite markers [75]. Moreover, the calculations in Douhovnikoff and Dodd were based on less than half of nuclear SSR markers (6 vs. 12), and only two were common in both studies.

The high variation within populations and low differentiation between populations are expected in wind pollinated tree species with no apparent physical or environmental barriers between populations, such as coast redwood [76]. The German samples in the data set G were closer to the northern rather than to southern populations based on both Jost’s D- and FST-values (S6C Table), which makes sense, although it was not statistically significant. As expected, the level of genetic diversity and number of private alleles in GER are less than in the original populations based on nuclear SSRs (Table 1), but not as low as detected typically for recently introduced tree species [77,78].

The STRUCTURE analyses based on the cpSSR markers were able to identify the San Francisco Bay as a border or geographic barrier between two main clusters in both data sets C and F. The STRUCTURE analyses within each of the two areas, north and south of the San Francisco Bay, did not reveal any additional clusters, neither in the dataset F nor C. Both Jost’s D- and FST-values are in consensus with the STRUCTURE and NJT results suggesting that CCJDF and AET in data set C and O, P and Q in data set F were the most differentiated groups compared to the remaining populations. AET represents the most interior location in data set C and surprisingly had higher genetic distance to the geographically closer populations in Napa Valley rather than to the southern populations (S6A Table).

Our data confirmed the San Francisco Bay as a border suggested already earlier by Brinegar [30] based on a single chloroplast marker. It is also identical with one of two borders identified by Sawyer et al. [79] based on the soil conditions and water availability provided by precipitation and fog. The low population structure might also be explained by transferring potentially non-local genotypes within the distribution range due to the long and intensive use of coast redwood as a timber species [80]. Trees from areas north of San Francisco Bay might have been planted in the south and vice versa. This uncertainty about the origin of genotypes and their adaptation to the local conditions should be considered in natural regeneration and conservation management strategies of coast redwood. Unfortunately, this is hard to detect and could only be verified in a comprehensive and detailed population genetic study of natural coast redwood populations in California and Oregon using genome-wide markers.

Unlike coast redwood, haplotype differentiation in other conifer species usually reflects the geographic origin of populations [81-83]. The lack of strong population structure and low genetic differentiation can explain the inconsistent assignment of German trees to populations in both reference data sets C and F and the low QI for them. The chloroplast genome is inherited paternally in coast redwood [34]. It does not prevent the cpSSRs markers to be used for the assignment, but it can be affected much by gene flow promoted by pollen and may not be as accurate as assignment based on nuclear markers. Combining the individual haplotypes closely connected together in the haplotype network into several clusters would not help for the individual assignments because individual trees composing the same cluster had different geographic origin. This conflict between geographic origin of an individual tree and position of its haplotype in the haplotype network, as well as discrepancy between data sets C and F in the assignments of German trees can be explained by the possibility that not the all sampled California stands in the data set C represented necessarily natural populations. The potential errors in the records regarding tree origin in the reference populations need to be also taken into account when considering the reliability of the assignment of individuals to an origin. Individuals with wrongly identified origin in the reference population can decrease the assignment quality [71]. However, for reliable assignment a stronger differentiation between populations in a reference data set and sufficient sample size of each reference population are needed [4,84,85]. Sample sizes were possibly insufficient for some populations in the reference C and F ranging from 4 to 47 individuals per population or watershed. However, all but five German trees were assigned correctly to the northern and southern clusters identified in the STRUCTURE analyses (Table 1, Figs 6 and 7 and S7 and S8 and S8 Table).

We also tried to find those trees in the data sets F and C which would be genetically similar to the German trees using pairwise individual similarity in their Q-value compositions to calculate pairwise genetic index of similarity and genetic distance between all trees and to visualize their clustering by generating NJTs and PCoA plots (S9–S12 Figs and S9 Table). Different German trees were similar to trees representing different watersheds (S9 and S11 Figs) and groups (S10 and S12 Figs). Some German trees were associated with the same trees or the same watershed or group in the F or C data sets for both K (for instance, SF75, CH1, etc.; Table 2). There is a trend of German trees to cluster more with trees from North watersheds and populations. These data can be used to trace the origin of the German trees, but cautiously. Additional verification using genome-wide markers or sequencing is needed to infer the true origin of the German trees.

The differentiation between populations based on the cpSSR markers is usually less than based on nuclear SSRs in conifers [86]; particularly, due to long distance migration by pollen and because the mutation rate in chloroplast microsatellite loci is lower [87]. In case of coast redwood, the paternal inheritance of chloroplasts [34] and the long distance gene flow via wind-dispersed pollen and seeds could explain the similarly low differentiation results for both marker sets. Ribeiro et al. [88] found similar results when comparing population differentiation based on AFLP markers with the one based on the paternally inherited cpSSR markers in the wind pollinated conifer Pinus pinaster Aiton. In addition, Petit et al. [89] showed that for various conifer species genetic differentiation based on bi-parentally inherited markers correlated with differentiation based on paternally inherited markers due to the similar gene flow vectors, but the latter was lower in general [90].

The reliable clone identification based on microsatellites in data set G (S6 Fig) confirmed results of previous studies based on allozyme and AFLP markers [27,91]. General difficulties to accurately genotype microsatellite markers in polyploid organisms excludes analyses based on their allele frequencies [35,92-95]. It concerns also coast redwood, but genotyping problems could be even more aggravated due to a high probability of somatic mutations in basal sprouting shoots in these extremely long living trees, which can result in different genotypes of different tissues and clones originated from the same tree [35]. The somatic mutations can help to maintain a high genetic diversity level in coast redwood populations, despite the low seed germination rate and high clonal growth. The assumption of Hardy-Weinberg equilibrium is also tricky due to common clonal growth in coast redwood populations [92], where trees within 40 m radius can belong to the same clone [91].

The correct estimation of the null allele frequencies and allelic dosage are two major difficulties associated with genotyping polyploid organisms using microsatellite markers [32,37]. In our study, the risk of null alleles should be reduced since null alleles are in general less frequent in EST-SSRs due to their location in more conserved regions [32], and the Bruvo distance used as the genetic distance measure in this study does not require same allelic dosages between individuals [43,45]. The fragment scoring following Pfeiffer et al. [42] and confirmation of allele scoring using multiple ramets of the same clones reduced the scoring errors due to the stutter effects in this study. Narayan et al. [35] considered only those alleles for the multilocus lineages (MLL) that were consistent in at least two different tissue types of the same tree, and they found a low error rate in allelic dosage in coast redwood based on the Bruvo distance. The possible number of genotypes in polyploid species, such as coast redwood, can be very high [42], hence, just a few microsatellite markers can be sufficient for clonal identification and population structure analysis.

The number of chloroplast haplotypes found in this study was exceptionally high (150) based on six cpSSR markers genotyped in 579 samples in total in all data sets. This extremely high number of chloroplast haplotypes was not observed earlier in any tree species. Similar high genetic diversity in chloroplasts was observed in another conifer species, Abies nordmanniana (Steven) Spach., with 111 haplotypes in 361 individuals, although the sampling range included a several times larger area than for coast redwood [96]. Although speculative, there could be several reasonable explanations such as plastid heteroplasmy and accumulation of somatic mutations during the long lifespan of coast redwood, which need additional research for verification. Plastid heteroplasmy has been observed in several plant species and can occur also due to bi-parental inheritance, hybridization, introgression, recombination between chloroplast genomes and plastid leakage [97,98]. According to Wolfe and Randle [99], heteroplasmy due to hybridization can be a stable condition. Additionally, long-living perennials and vegetatively reproducing plants can maintain a high genetic variation level in chloroplasts [99]. Scott et al. [100] did not rule out an allopolyploidization event within the “Sequoia-clade” as a source of the polyploidy of coast redwood. Possible hybridization events, longevity of trees and clonal propagation as a source and maintenance of chloroplast heteroplasmy apply to coast redwood. One can expect multiple fragments representing different chloroplast alleles amplified due to chloroplast heteroplasmy, but we observed exclusively single fragments in this study. However, it might be due to the uneven allelic dosage within a sample, and as a result only one allele (fragment) is amplified during PCR [99]. To our knowledge, we are unaware of any studies specifically demonstrating correlation between polyploidy and chloroplast heteroplasmy in plants; although, in the genus Medicago both polyploidy and chloroplast heteroplasmy are common [101]. However, finding the reason for the exceptionally high number of haplotypes in coast redwood requires further investigations.

Conclusions and future directions

Coast redwood forests used to have a continuous distribution along the Pacific Coast in California before intensive logging started at the beginning of the nineteenth century [102]. Therefore, current coast redwood populations are considered as remnant and fragmented populations. However, being long living clonal trees with a high somatic mutation rate, coast redwood maintained its high genetic diversity despite multiple bottlenecks [15,34]. Although a limited number of German trees were introduced only decades ago, they have relatively high genetic diversity.

The combination of low sexual reproduction and local adaptation could cause insufficient ability to meet the challenges of climate change and will increase the pressure on coast redwood [9]. California has become drier in the last 2000 years [16], and the very important fog has been declining in its frequency during the last century [39,103]. Considering these threats O’Hara et al. [15] emphasized the necessity to find drought tolerant genotypes and to increase genetic diversity, especially for southern populations. Genetic and physiological mechanisms behind drought resistance are similar to those that are behind frost tolerance, therefore geographic variation in drought tolerance in tree species often overlaps with variation in frost tolerance [104-106]. The identification of water stress resistant genotypes would benefit both Californian and German forestry. Tolerant genotypes would not only provide Germany with a valuable timber species considering climate change, but also presents suitable resources for ex-situ conservation programs for coast redwood, which was already suggested for the sister species Sequoiadendron giganteum (Lindl.) J. Buchholz [8].

Further studies of the genetic structure of coast redwood populations based on functional markers that are potentially under selection, such as non-synonymous SNPs or SNPs in regulatory genes, can help to identify local adaptation to particular environmental stress factors of interest [107]. More studies should be also done using mitochondrial DNA markers, which can help to infer maternally based gene flow and lineages due to their maternal inheritance. They are distributed by seeds and should be less affected by long-distance dispersal than pollen [108]. The results obtained in these prospective studies will provide additional important data for sustainable timber production and conservation management, in situ and ex situ, for coast redwood, especially in the context of predicted climate change.

Map of the original “Kuser’s” samples in the data set F.

Original locations of the “Kuser’s” samples in the data set F are presented on the map together with the mean monthly temperature pattern for the time period 1979–2013 indicating colder and warmer temperatures by darker and lighter shades of grey (http://chelsa-climate.org).

(PDF)

Click here for additional data file.

Map of the German sampling locations (data set G).

Darker and lighter shades of grey highlight colder and warmer mean monthly temperatures for the time period 1979–2013 (http://chelsa-climate.org).

(PDF)

Click here for additional data file.

Neighbour-joining tree of 84 ramets representing 30 different clones (with 2–3 ramets per clone) provided by the Allerweltsgrün nursery (Köln).

It is based on 12 nSSR markers with the final ranking combination obtained according to Pfeiffer et al. [45] and the ‘bruvo genetic distance’ [48]. Numbers indicate bootstrap values (not percentage). The sample ID represents the clone ID followed by the ramet consecutive number.

(PDF)

Click here for additional data file.

Neighbour-joining tree of the “Kuser’s” samples in the French (F) data set (St. Fargeau) grouped into 17 watersheds (A-Q).

It is based on six cpSSR markers and Nei’s genetic distance ([50] after [51]) with 1000 bootstraps. Numbers indicate bootstrap values in terms of percentage.

(PDF)

Click here for additional data file.

Neighbour-joining tree of the 16 Californian reference populations represented by samples collected in 2017 (data set C).

It is based on six cpSSR markers and Nei’s genetic distance ([50] after [51]) with 1000 bootstraps. Numbers indicate bootstrap values in terms of percentage. County names of sampling sites are in brackets.

(PDF)

Click here for additional data file.

Neighbour-joining tree of the 143 individual trees collected in Germany (data set G).

It is based on 12 nSSR markers with the final ranking combination obtained according to Pfeiffer et al. [45] and the ‘bruvo genetic distance’ [48] with 1000 bootstraps. Numbers indicate bootstrap values (not percentage).

(PDF)

Click here for additional data file.

STRUCTURE results for the reference data sets C and F based on 12 nSSR markers.

Plots for data sets C and F are presented for the most likely number of clusters (K) inferred without (a) and with (b) the “locprior” function engaged (K = 6 and K = 2 for C and K = 2 for F, respectively). L(K) and ΔK statistics generated by the ClumPAK software are also presented.

(PDF)

Click here for additional data file.

STRUCTURE statistics for the reference data sets C and F based on 6 cpSSR markers.

Presented are the most likely number of clusters (K) and results of STRUCTURE runs for K from 1 to 20 or 25 with 20 iterations each for the complete data set C (a), northern subset C (b), southern subset C (c), complete data set F (d), northern subset F (e), and southern subset F (f).

(PDF)

Click here for additional data file.

Neighbour-joining trees (NJTs) based on pairwise individual tree genetic distances between all 396 trees in German (G), French (F) and Californian (C) data sets generated using Q-values for K = 17.

(PDF)

Click here for additional data file.

Neighbour-joining trees (NJTs) based on pairwise individual tree genetic distances between all 396 trees in German (G), French (F) and Californian (C) data sets generated using Q-values for K = 30.

(PDF)

Click here for additional data file.

Principal coordinate analysis (PCoA) plots based on pairwise individual tree genetic distances between all 396 trees in German (G), French (F) and Californian (C) data sets generated using Q-values for K = 17.

(PDF)

Click here for additional data file.

Principal coordinate analysis (PCoA) plots based on pairwise individual tree genetic distances between all 396 trees in German (G), French (F) and Californian (C) data sets generated using Q-values for K = 30.

(PDF)

Click here for additional data file.

Q-plots based on the STRUCTURE results for K = 17 and K = 30 for all 396 trees in German (G), French (F) and Californian (C) data sets.

(PDF)

Click here for additional data file.

ID, geographic and watershed data for the French (F) set of the “Kuser’s” samples.

(DOCX)

Click here for additional data file.

Original latitude and number of samples representing 17 watersheds in the French (F) data set (St. Fargeau) and the Russell Reserve according to Douhovnikoff and Dodd [34].

(DOCX)

Click here for additional data file.

Geographic and climatic data for sampled locations in California (data set C).

(DOCX)

Click here for additional data file.

Sampling site, clone ID, and geographic location for samples collected in Germany (data set G).

(XLSX)

Click here for additional data file.

Data for 18 microsatellite (SSR) markers used in the study.

(DOCX)

Click here for additional data file.

Pairwise Jost's D- (in the left corner lower triangle) and FST- (in the right upper corner triangle) values based on 12 nSSR marker for 16 populations of data set C (A), 17 watersheds of data set F (B), and among three groups NORTH, SOUTH, and GER, respectively (C). All values were insignificant.

(XLSX)

Click here for additional data file.

Verification of the 18 PCR primer pairs of the 18 microsatellite (SSR) markers by mapping their primer nucleotide sequences to the draft coast redwood nuclear and chloroplast genome assemblies.

(XLSX)

Click here for additional data file.

The German samples in the data set G assigned to the populations and watersheds in the reference data sets representing samples from California (C) and the “Kuser Provenance test” in St. Fargeau, central France (F), respectively, and assignment quality scores (%) for two ranks.

(DOCX)

Click here for additional data file.

Q-matrices for K = 17 and K = 30 and genetic distance matrices based on them.

(XLSX)

Click here for additional data file.