Genetic evidence of human mediated, historical seed transfer from the Tyrolean Alps to the Romanian Carpathians in Larix decidua (Mill.) forests

Abstract

Key message

Historic transfer of larch from Alpine sources to Southern and Eastern Carpathians has been verified by means of nuclear genetic markers. Tyrolean populations can be differentiated into a north-western and south-eastern group, while Romanian populations are separated according to the Southern and Eastern Carpathians. Low-level introgression from Alpine sources is found in autochthonous Carpathian populations.

Context

Large scale human mediated transfer of forest reproductive material may have strongly modified the gene pool of European forests. Particularly in European larch, large quantities of seeds from Central Europe were used for plantations in Southern and Eastern Europe starting in the mid nineteenth century.

Aims

Our main objective was to provide DNA marker based evidence for the anthropogenic transfer of Alpine larch reproductive material to native Carpathian populations.

Methods

We studied and compared 12 populations (N = 771) of Larix decidua in the Alps (Austria, Italy) and in the Southern and Eastern Carpathians (Romania) using 13 microsatellites.

Results

High genetic diversity (H_e = 0.752; R_S = 9.4) and a moderate genetic differentiation (F_ST = 0.13; G′_ST = 0.28) among populations were found; Alpine and Carpathian populations were clearly separated by clustering methods. A Tyrolean origin of plant material was evident for one out of four adult Romanian populations. In the transferred population, a genetic influence from Carpathian sources was found neither in adults nor in juveniles, while the natural regeneration of two Romanian populations was genetically affected by Alpine sources to a minor degree (2.2 and 2.9% allochthonous individuals according to GeneClass and Structure, respectively).

Conclusion

Tracing back of plant transfer by means of genetic tools is straightforward, and we propose further studies to investigate gene flow between natural and transferred populations.

1 Introduction

European larch (Larix decidua Mill.) is a monoecious, anemophilous pioneer tree (Pinaceae) with a disjunct natural distribution in Europe. This deciduous conifer naturally occurs in the montane to the subalpine altitudinal belt, but can also be found in the Polish lowlands; in addition to pure stands, it occurs also in mixed stands associated with stone pine (Pinus cembra), Norway spruce (Picea abies), Swiss pine (Pinus mugo), Silver fir (Abies alba), and European beech (Fagus sylvatica) (Mayer 1992). From the forelands of the Alps and Carpathians, which represented the main refuge areas of the species during the last glacial period, it has re-colonised Central Europe (Wagner et al. 2015b). The natural distribution is split into five regions: (1) Alps, (2) Eastern Sudetes, (3) Polish lowlands, (4) Western Carpathians with scattered populations in the High and Low Tatra Mountains as well as Beskids, Fatra, and Ore Mountains, and finally (5) scattered occurrences in the Eastern and Southern Carpathians and in the Apuseni Mountains of Romania (Wagner et al. 2015a). Differences in seed weight (Simak 1967), stomata rows (Maier 1992), and female flower colour (Geburek et al. 2007) exist among populations and provide an indication about their geographical origin.

European larch is an economically important forest tree species mainly used for its timber, while its turpentine as well as larch shingles are locally also of importance (Øyen 2006). In the second half of the nineteenth century, the rapid establishment of the railway system in Central and Northern Europe enhanced the trade with plant material and triggered large-scale cultivation outside and inside of the native range (Pardé 1957; Jansen and Geburek 2016). In consequence, the current distribution area of L. decidua has more than doubled to 1.5 million ha in Europe (Pâques et al. 2013); the species has been extensively cultivated in temperate European forests especially in Germany, Poland, France, Denmark, Great Britain, Sweden, and Norway.

In this study, we focused on L. decidua in the Southern and Eastern Carpathians, where it was affected by an intense seed transfer from other areas of the native range, especially from the Tyrolean Alps. In particular, forests in Transylvania and along the Prahova Valley were affected by these transfers starting in the middle of the nineteenth century (Rubţov 1965). Cultivations of unknown origin were reported in Transylvania (Sighişoara), the Southern (Bucegi Mountains, region of Azuga) and Eastern Carpathians (Dofteana), as well as south of the Carpathian arc (Câmpina, Piteşti) (Gava 1963; Rubţov and Mocanu 1958). Although the origin of this material was not recorded, the use of Alpine material is likely. During that time, Transylvania was part of the Austro-Hungarian Empire and cultivations with Alpine seed sources were promoted (Rubţov and Mocanu 1958). Also, the Austrian seed trading companies “Jenewein” (Innsbruck) and “Steiner” (Wiener Neustadt) advertised the use of larch material especially from Tyrol and the Vienna Woods in the Southern Carpathians (Rubţov 1965). Since 1890, the trade of Alpine plant material has been documented and Austria was probably the main supplier (Gava 1963). However, which Alpine sources were used and where exactly the material was planted in Romania remains unknown. Historic records indicate, that the transferred material originated mainly not only from the Mieminger Plateau (North Tyrol, Austria) and the Etsch and Eisack valleys (South Tyrol, Italy) (Gothe 1961) but also from South (e.g. St. Johann, the Fiemme valley) and North Tyrol (e.g. Götzens, Seefeld) (Jansen and Geburek 2016).

Larch plantings with transferred Alpine provenances, especially when seeds from the Western Alps were used outside of the Alpine range, were reported to perform poorly (Weisgerber and Šindelár 1992). Through such a plant transfer the gene pool of recipient autochthonous larch populations may have been adversely affected. Outbreeding depression may be one of the effects, i.e. the reduction in progeny fitness after populations have been mixed compared to the fitness of progenies from crosses between trees that are more closely related (Frankham et al. 2017). Theoretically such effects resulting from chromosomal incompatibilities causing (partly) sterility, deleterious epistatic interactions between diverged genes, and/or disruption of co-adapted gene complexes cannot be ruled out in transferred larch populations, but these effects are not very likely. Instead, local adaptation losses are more probable due to the introduction of maladapted genes from Alpine larch sources. Theory predicts that the genetic impact on the native population depends on the degree of maladaptation of the transferred plant material and on the size of the recipient population (Kopp and Matuszewski 2014; Kremer et al. 2012). In any case, a transfer can cause genetic swamping, i.e. loss of the integrity of the local gene pool as a result of the introgression of genes of the transferred population through extensive hybridization with the local source (Hufford and Mazer 2003).

Romanian larch populations are scattered and hence are often genetically isolated. Located at the southern fringe of the native range they represent rear-edge populations and are therefore important from a gene conservation point of view (Fady et al. 2016). Therefore, we addressed the following research questions based on data from nuclear microsatellite markers: Can a potential Alpine origin be verified and—if so—can the Alpine sources be narrowed down to certain regions? Are there detectable effects on the gene pool of the Carpathian populations from Alpine seed sources in adults and/or juveniles? And if this is true, how large is the genetic impact on a local scale? Moreover, our data are of relevance for natural population differentiation based on in situ samples in the studied regions.

2 Materials and methods

2.1 Population sampling

Needle or cambium samples were collected from 12 populations (in situ) from the Alps and Carpathians. In the Alps material was collected in North Tyrol (Austria) and South Tyrol (Italy) (hereafter, the denotation “Tyrol” stands for both South and North Tyrol); Carpathian material was sampled in the Vâlcea, Prahova, and Braşov counties of Romania representing the native population groups, except for R3, which presumably originated from a Tyrolean source (Fig. 1, Table 1). The sample sites were selected based on documented evidence (Gothe 1961; Rubţov and Mocanu 1958) and the help of local experts. Additionally, we sampled individual trees in three locations (spot checks, SC; overall N = 10; Table 1) in the vicinity of the sampled native stands due to their Tyrolean-like habitus (Geburek et al. 2007; Rubner and Svoboda 1944). Plant tissue was collected from 50 adult dominant trees in forest stands, considering at least 30 m of inter-tree distance to avoid sampling of closely related individuals. To determine the potential recent influence of the Tyrolean variety on the Carpathian gene pool, we furthermore sampled 50 juveniles (diameter at breast height < 10 cm or height < 1.3 m) from the natural regeneration in each of these stands. Geographical coordinates were recorded with a 60CSx GPS device (Garmin International, Inc.). Approximately 3-cm² sized cambium samples at the stem base were collected by using a hollow punch, while mainly needle samples were taken from juveniles. Collected materials were dried and stored in zip lock plastic bags containing silica gel.

Fig. 1

Sampled populations of Larix decidua within the natural distribution range (dashed line). Colours indicate geographical regions (blue—Austria (North Tyrol), green—Italy (South Tyrol), red—Romania). Map from Wagner (2013) was modified based on data from Stănescu et al. (1997)

Table 1 Geographical origin of 12 Larix decidua populations used in this study

2.2 DNA extraction and microsatellite genotyping

Total DNA was extracted from approximately 100 mg of cambium (dried in silica gel) or needle tissue, after freezing in liquid nitrogen and homogenisation using a Qiagen Tissue Lyser device (Qiagen Inc.). Genomic DNA was extracted by using the DNeasy 96 Plant Kit (Qiagen Inc.) following the manufacturer’s protocol. Quality and concentration of obtained DNA were measured using an ND-1000 spectrophotometer (NanoDrop Inc.). DNA was stored at 4 °C. All samples were genotyped with 13 highly polymorphic nuclear microsatellite loci (Wagner et al. 2012). Multiplex PCR amplification was optimised to be performed in a 10-μl reaction volume containing 2–10 ng of genomic DNA, 5-μl HotStarTaq Master Mix (Qiagen Inc.), double distilled water, and 0.3 μM of forward and reverse primers each. We used the following cycling protocol on a TC-412 Programmable Thermal Controller (Techne): 35 cycles with 94 °C for 30 s, 56 °C for 90 s, and 72 °C for 60 s. Before the first cycle, a prolonged denaturation step (95 °C for 15 min) was included and the last cycle was followed by a 30-min extension at 72 °C. Genotyping was outsourced to a commercial provider (ecogenics, Balgach, Switzerland) using Applied Biosystems (Foster City, CA, USA) chemistry and allele calling tools with manual checking of scores.

2.3 Data analyses

Genotypes with missing data for more than two loci were excluded from the analysis. Individuals from SC1, SC2, and SC3 were only used for Structure and GeneClass analyses (see below) (Table 1).

2.3.1 Genetic diversity within populations

Microsatellite data were analysed for allele dropout, null alleles, and slight changes in allele sizes during PCR amplification using Micro-Checker 2.2.3 (van Oosterhout et al. 2004). Furthermore, the Bayesian approach of individual inbreeding coefficient F_i implemented in INEst 2.2 (Chybicki and Burczyk 2009) was used to check raw data for null alleles, to estimate unbiased inbreeding coefficients F_i, and unbiased observed and expected heterozygosity within populations (settings: 500,000 Markov chain Monte Carlo iterations, burn-in of 50,000, thinning parameter of 50). Deviations from Hardy-Weinberg expectation (HWE) expressed as heterozygote excess or deficiency and the genotypic disequilibrium among pairs of loci were assessed using Genepop 4.7.0 implementing Monte Carlo Markov chain simulations of Fischer’s test with 10,000 dememorisations, 100 batches, and 10,000 iterations (Raymond and Rousset 1995; Rousset 2008).

Genetic diversity indices, including deviations from HWE, number of alleles, effective number of alleles (N_e), observed (H_o), expected heterozygosity (H_e), and number of private alleles (A_p) were calculated per locus and as means over all loci with corresponding standard errors (SE) using GenAlEx 6.502 (Peakall and Smouse 2006, 2012). The number of alleles corrected for equal sample size (allelic richness R_S) with a rarefaction to 39 individuals (minimum sample size; Table 1) was calculated using Fstat 2.9.3.2 (Goudet 1995). Descriptive summarised statistics (e.g. mean, SE) were performed using R 3.3.2 (R Development Core Team 2016).

Additionally, we characterised the spatial genetic structure (SGS) within the (sub)populations using the Nason estimator of multilocus kinship coefficient F_ij (Loiselle et al. 1995) implemented in SPAGeDi 1.5 (Hardy and Vekemans 2002). Ten distance intervals were defined by the program so that the number of pairwise comparisons within each distance interval was approximately constant. Individual sampling locations and gene copies were each permuted 10,000 times. SE of mean F_ij values per distance class was generated by jack-knifing over loci. We assessed the significance of estimated values by comparing them with the confidence interval of the coefficients under the assumption of no SGS. Natural and non-natural (e.g. afforestation) populations should be different as in planted populations SGS assessed with selectively neutral markers should be random. Evidence of recent bottlenecks within populations was assessed using the graphical test of Luikart et al. (1998) implemented in Bottleneck 1.0.02 (Piry et al. 1999), which in most cases is not sensitive to deviations from Hardy-Weinberg proportions (Luikart et al. 1998).

2.3.2 Genetic diversity among populations

To determine population differentiation and the relationships among populations, we first performed a principle co-ordinate analysis (PCoA) for visualising population clusters due to relationships among inter-individual genetic distances and to identify a set of reduced dimension traits (eigenvectors) according to the number of populations using GenAlEx (Peakall and Smouse 2006, 2012). In order to take into account the unequal sample sizes, we used Nei’s unbiased genetic distance (Nei 1978) as algorithm for the PCoA. Secondly, we computed an unweighted pair group method arithmetic average dendrogram (UPGMA) based on Nei’s standard genetic distance (Nei 1972) using 1000 bootstrapped matrices created by a microsatellite analyser (Dieringer and Schlötterer 2003). For consensus tree construction, we used the programs Neighbour and Consense implemented in the Phylip 3.63 package (Felsenstein 1989).

Furthermore, we used Structure 2.3.4 (Pritchard et al. 2000) to infer the most likely number of population clusters and attempt to assign individuals to populations by reference to their genotypes. The software uses a Bayesian clustering algorithm to pool individuals into a predefined number of clusters (K) by minimising deviations from HWE and gametic-phase disequilibrium within the clusters. We used the admixture-model, where each individual does not have any information about population affiliation, with K values ranging from 2 to 10 and a run length of 800,000 iterations with a burn-in period of 200,000. Four runs per K were performed for reasons of iteration. In these computations, we included all available samples, also groups (spot checks) with a sample size < 10 (cf. Table 1). To predict the appropriate number of clusters for the whole data set, we used ΔK, an ad hoc quantity related to the rate of change of the log probability of data with respect to the number of clusters (Evanno et al. 2005), implemented in Structure Harvester 0.6.94 (Earl and von Holdt 2012). We defined individuals with more than 30% probability to belong to a different genetic cluster as introgressed. In addition, following the analysis of the whole data set, also putative autochthonous Romanian and Tyrolean adult populations were run separately with the same settings to analyse unbiased population structure on the regional level.

Information about population differentiation at various levels of population aggregations was obtained using an analysis of molecular variance (AMOVA) implemented in Arlequin 3.5 (Excoffier and Lischer 2010), with which global F_ST (Weir and Cockerham 1984) and R_ST (Slatkin 1995) were calculated separately with 10,000 permutations. Thereby, we defined the population aggregations according to the two major groups (Tyrolean vs. Romanian populations), as well as to their demographic grouping (adults vs. juveniles). To take allele size and stepwise mutations into account, we compared the observed differentiation F_ST with R_ST. If stepwise mutations have contributed to differentiation, R_ST is expected to be significantly higher than F_ST. In addition, we calculated G′_ST (Hedrick 2005), which unlike the G_ST value (Nei 1972), also regards the restricted amount of genetic variation given by overlapping sets of alleles among subpopulations (Hedrick 1999); this was done using GenAlEx (Peakall and Smouse 2006, 2012) with 10,000 permutations and bootstraps, respectively. To assess the effect of geographical proximity (and perhaps population history) on genetic diversity of the two regions, we compared diversity indices (R_S, H_o, H_e, F_is; using Fstat [Goudet 1995]) and differentiation (F_ST, G′_ST; using a specific R script [R Development Core Team 2016], C. Dobeš, unpublished) among groups. Populations were grouped according to the Structure analysis, and significance was assessed based on 10,000 permutations.

Finally, we analysed the existence of first-generation migrants in the Romanian juvenile populations using GeneClass 2.0 (Piry et al. 2004). Thereby, we used the Bayesian criterion of Baudouin and Lebrun (2000) for computing the likelihood of occurrence of the individual genotype within the population where the individual has been sampled. L_home as the test statistic was used in order to prevent a potential bias in likelihood ratios caused by missing source populations (ghost populations) (Paetkau et al. 2004). For probability computation and better performance of calculation in case of ghost populations (Leblois 2011), we chose the Monte-Carlo re-sampling algorithm of Cornuet et al. (1999) with 10,000 simulated individuals and a type I error (probability of detecting a resident as immigrant) of < 0.05.

3 Results

3.1 Genetic diversity within populations

Average null allele frequency across all loci per population ranged between 1.9 and 8.5% (mean of 0.045; SE = 0.004; data not shown). These values fall well within the range between 5 and 8%, which according to a simulation study (Chapuis and Estoup 2007) have only minor effects on classical estimates of population differentiation. Neither large allele dropouts nor scoring of stutter peaks were identified in any of the loci. We found no evidence for genotypic disequilibrium between loci within populations (P < 0.05). However, all populations showed a significant deficit of heterozygotes, except population T9 and T12 in North Tyrol, thus in 11 populations highly significant (P < 0.001) deviations from HWE were detected (Table 2). In total, 113 alleles, ranging from six (locus Ld101_565) to a maximum of 38 (locus bcLK263_550) with an average of 18.46 alleles per locus were detected. The lowest mean number of observed alleles (N_a) was calculated in population R2 (6.77), whereas the highest value was found in ST5 (11.23). Expected heterozygosity (H_e) values were mainly below values of observed heterozygosity (H_o) also indicated by positive fixation indices. Generally, lower diversity values were found in Romanian compared to Tyrolean samples (Table 3). We detected a total of 32 private alleles at low frequency in 13 populations with a minimum of one private allele both in R3 and T10 and a maximum of five private alleles in ST7. In most cases, SGS did not deviate from random expectation and significant F_ij values were only obtained in first distance classes, allowing no estimation of gene dispersal distances. The mean pairwise kinship coefficient F_ij (Loiselle et al. 1995) over all loci was significant with 0.015 in the first distance class in Romanian and with 0.010 in Tyrolean populations in each first distance class. The highest, significant F_ij detected was 0.025 in the first distance class (up to 100 m) in R1 (juveniles). Moreover, R1 (adults, 81 m), R2 (juveniles, 72 m), R3 (adults, 58 m; the planted stand), R4 (juveniles, 11 m), as well as ST7 (70 m), ST8 (133 m), and T9 (67 m) showed significant positive deviations from the permuted mean values in the first distance class. Recent bottlenecks in each population were not detected.

Table 2 Genetic diversity of all Larix decidua populations and age cohorts based on 13 SSR markers studied

Table 3 Differences in genetic diversity between Tyrolean and Romanian populations

3.2 Genetic diversity among populations

The PCoA based on Nei’s unbiased genetic distance (Nei 1978) revealed clusters for Romanian and Tyrolean populations. However, a strong affiliation of the Romanian population R3 to Tyrolean populations was evident (Fig. 2). The first two axes of the PCoA explained 75.7 and 11.2% of the observed variation, respectively. Similar results of relationships among populations according to geographic proximity were provided by the UPGMA dendrogram, where each region, i.e. Romania, North and South Tyrol formed its own clade (Fig. 3).

Fig. 2

Principal co-ordinate analysis showing the multivariate relationships of 12 Larix decidua populations. The first and second axes explain 75.7 and 11.2% (86.9%) of the observed variation, respectively. Circles indicate possible clusters of related populations

Fig. 3

UPGMA consensus tree showing relationships between 12 Larix decidua populations. Numbers beside the branches indicate bootstrap support for nodes. Circles indicate possible grouping of related populations

Results of Structure consistently showed high assignment coefficients for Tyrolean and Romanian clusters, resulting from strong differentiation between these population groups (Fig. 4a). Results interpreted using the method of Evanno et al. (2005) revealed two clusters (K = 2) as the most likely group supporting the results from PCoA as well as from the UPGMA dendrogram. The Tyrolean cluster can be further differentiated into two groups, i.e. a north-western (T9, T10, T11, T12, ST8) and a south-eastern group (ST5, ST6, ST7) (Fig. 4b). Also, in the pure Romanian populations (excluding R3), a strong clustering was present, especially when these were analysed separately (Fig. 4c), corroborating the findings from PCoA and UPGMA. Population R3 and spot check SC3 (in close vicinity to R2) growing in Romania were assigned to the north-western Tyrolean cluster, while SC1 and SC2 (in close vicinity to R1) were assigned to the south-eastern cluster. According to our threshold, we found three juveniles in R1 and one individual both in the juveniles and adults, respectively, in R4 that were assigned to the Tyrolean clusters (Fig. 4a).

Fig. 4

a Bar plots of individual population assignment performed with Structure for all sampled Larix decidua individuals for two assumed ancestral clusters (K = 2). Immigrants were identified both by Structure and GeneClass are marked by a black arrow; immigrants identified by Structure alone by a blue arrow. b Bar plots of individual population assignment performed with Structure for the Tyrolean populations with K = 2. c Bar plots of individual population assignment performed with Structure for the native Romanian populations with K value of 3

The AMOVA revealed an F_ST value of 13.3% and an R_ST value of 8.01% within the adults when groups were defined according to the Structure analysis (Tyrolean vs. Romanian populations; Fig. 4). Within the juveniles, an F_ST value of 13.6% and an R_ST value of 9.22% were calculated. G′_ST was estimated within adults as 27.7 and 32.6% for juveniles, respectively (Table 4).

Table 4 Results of analysis of molecular variance separated for adults and juveniles

For the GeneClass assignments, we checked whether first-generation migrants of Tyrolean origin could be identified in the Romanian putatively autochthonous juveniles (natural regeneration). Here, we found two individuals in R1 (juveniles) to be linked to T11 and R3, respectively, and one individual in R4 (juveniles) linked to ST7. Of particular interest is that in R3 (which is obviously from north-eastern Tyrolian origin), no first-generation migrants from autochthonous Romanian populations were detected.

4 Discussion

4.1 Genetic diversity within populations

Our results of nuclear microsatellite diversity are similar to other studies on L. decidua (King et al. 2013; Pluess 2011; Wagner et al. 2012) and to other conifer species with similar habitat demands in boreal or temperate forest associations such as P. abies (King et al. 2013; Tollefsrud et al. 2009), P. cembra (Dzialuk et al. 2014; Lendvay et al. 2014), or Pinus sylvestris (Pazouki et al. 2016). In general, high levels of diversity within populations in long-lived, outcrossing, and late successional coniferous taxa were obtained. However, such comparisons are always limited as samples sizes and markers differ among studies. Nevertheless, comparisons with the results of Wagner (2013) obtained by the same marker-set showed similar results of F_ST values, but slightly higher inbreeding coefficients than ours, which might be the effect of the presence of null alleles. Similar to Wagner (2013), in our study marker bcLK263 and bcLK211 had the highest and Ld42 had the lowest number of alleles. Comparisons with species within the genus Larix reveal higher microsatellite diversity for L. decidua (H_e = 0.75) than for L. lyallii and L. occidentalis (0.42 and 0.58; Khasa et al. 2006) as well as for L. gmelinii (0.41–0.60; Zhang et al. 2015), L. sibirica, and L. cajanderi (0.63 and 0.56; Oreshkova et al. 2013). H_e estimates in L. kaempferi (0.72–0.76; Nishimura and Setoguchi 2011) are similar to our results. In principle, differences can be due to historical bottlenecks affecting populations besides other evolutionary processes.

Almost all populations show a significant homozygote excess compared to HWE. Similar results are also reported in other studies of Larix spp. (Larionova et al. 2004; Nishimura and Setoguchi 2011; Oreshkova et al. 2013). However, when null alleles in our study were taken in account as proposed by Chybicki and Burczyk (2009), only a very small homozygote excess was observed. Unexpectedly, we did not find a general trend that in juveniles an excess of homozygotes at neutral gene loci was more pronounced than in adults, as to be expected because inbred individuals are predominantly selected against. Juvenile cohorts often differ genetically from adults, as diversity (e.g. loss of rare alleles; Kettle et al. 2007) or homozygote excess (Fujio and von Brand 1991; Stoeckel et al. 2006) might decrease during growth depending on the different selection rate against homozygotes (deleterious homozygote recessive alleles). However, it remains open whether early postzygotic selection due to embryonic lethals (recessive lethal genes that eliminate most selfed embryos during seed development; Savolainen et al. 1992) or during later ontogenetic stages other forms of viability selection occurred as in our study we did not analyse seed material but rather established young trees.

Genetic diversity was lower in the Carpathians compared to the Tyrolean populations, especially in the south-eastern cluster. Larch forests were present in the Carpathians approximately from 11,500 years ago (Magyari et al. 2012). At the same time also in Tyrol, larch forests became established as macrofossil-pollen records indicate (Heiss et al. 2005; Oeggl and Wahlmüller 1994). Lower genetic diversity of the Romanian autochthonous populations (R1, R2, R4) might be explained by the origin of these populations from small, isolated, and genetically depauperate refugia in the western part of the Southern Carpathians (Wagner et al. 2015b). On the other hand, the gene pool of Tyrolean larch might have become enriched by the merging of populations derived from several large Alpine refugia, as the data of Wagner et al. (2015b) suggest, but this would need an in-depth analysis based on of a larger set of Alpine populations. Small differences in genetic diversity in Tyrol shown by the decrease of the number of private alleles from South (12) to North Tyrol (7) could be due to genetic bottlenecks affecting diversity during migration into the Alps. We did not find any evidence of recent bottlenecks in the studied populations; however, this does not mean that there were no former bottlenecks, as our sample sizes might have been too small to detect them (Luikart et al. 1998).

Generally, SGS did not deviate significantly from random expectation. Only in very few cases, statistically significant—but very weak—clumping was observed. In the planted population R3, no deviation from random SGS was expected. However, we found a statistically significant Loiselle’s coefficient in the first distance class in R3 (adults). We interpret this significant deviation as well as those in other populations as biologically not meaningful as Loiselle’s coefficients never exceeded 0.025. Our sampling was not ideal to compare SGS among populations as distance classes could not directly be compared, as adults and juveniles were partly found in patches rather than evenly distributed across space. Disregarding these methodical limitations, it is noteworthy that larch seeds are already released in spring and may be blown over long distances on smooth snow surfaces. Also low, but significant Loiselle’s kinship coefficients (< 0.05) in distance classes up to 40 m have been already reported in L. decidua (Pluess 2011). Seed dispersal in tamarack (L. laricina) (Brown et al. 1988) is probably more limited than in L. decidua, which may explain a more pronounced non-random spatial structure in naturally regenerated tamarack stands (Knowles et al. 1992).

4.2 Genetic diversity among populations

A high level of genetic diversity and moderate differentiation among populations and groups were found in the AMOVA analysis. This pattern is in accordance with other outcrossing, wind pollinating conifers (Porth and El-Kassaby 2014). L. decidua is similarly or only slightly higher differentiated than P. abies (Androsiuk et al. 2013) or P. sylvestris (Belletti et al. 2012). Our G′_ST of approximately 30% is similar to that found in Larix lyallii (G′_ST = 0.26; 19 populations studied), but much higher than in L. occidentalis (G′_ST = 0.13; 9 populations studied) (see Supplemental material in Heller and Siegismund 2009, based on data from Khasa et al. 2006). Considering the allele size, R_ST could not explain more of the observed variance than F_ST, suggesting the absence of a phylogeographic pattern. Our results also indicate (corroborated by UPGMA, PCoA and Structure analyses) that the genetic differentiation is higher in the Southern and Eastern Carpathians than in the Tyrolean Alps. While in Tyrol a continuous distribution is found, the distribution of larch in the Southern and Eastern Carpathians is restricted to small and isolated populations probably limiting gene flow and thus enhancing population differentiation.

Our analysis confirms the strong range wide structure reported by Wagner et al. (2015a) by analysing material collected in a provenance trial covering a gradient over the native distribution range on the basis of mitochondrial and nuclear markers. Our study reveals a more detailed genetic differentiation in the Tyrolean Alps and in the Carpathians. We have provided evidence that in putative autochthonous Tyrolean populations two slightly contrasting genetic clusters exist. These two genetic clusters (northern vs. southern populations) reflect the post-glacial immigration of larch (Wagner et al. 2015b). Today’s Tyrolean populations survived probably close to Verona and migrated northwards through the Etsch valley and then split into different valleys. A major split into a western (following the river Etsch) and eastern route following the river Eisack probably occurred where today the city of Bozen is located. This hypothetical migration pattern is supported by our data.

4.3 Translocation evidence

It is obvious from our results that population R3 is of north-western Tyrolean origin. The Eastern and Southern Carpathians were in particular affected by seed transfers originating from Alpine (Tyrolean) sources (Jansen and Geburek 2016). This transfer started already in the middle of the nineteenth century (Gava 1963; Rubţov 1965; Rubţov and Mocanu 1958). Therefore, we expected that the natural regeneration of Romanian putatively autochthonous stands may have been affected by the Tyrolean gene pool. This was the case to a minor extent for R1 and R4 as shown by Structure (2.9% allochthonous individuals) and GeneClass indicated 2.2% migrants from Tyrolean origin into the Southern and Eastern Carpathians. However, a general impact is impossible to be assessed from our study as any genetic introgression depends on the geographical and ecological distance as well as on the size of introduced populations. However, for R1, at least seven trees of Tyrolean origin (SC1 and SC2) were in close or relative close neighbourhood to R1 and at least three allochthonous trees (SC3) were in close proximity to R2; hence, we expected a higher gene flow from these sources. Unfortunately, no detailed information on the extent of allochthonous larch plants used in this area was available. Between 5 and 10% of intraspecific hybrids were detected in the offspring in P. abies stands in Norway, although only 4% of allochthonous trees was initially present (Dietrichson 1991). Introgression was also studied both for Spanish autochthonous Pinus pinaster and P. sylvestris populations surrounded by 15-fold and 3-fold, respectively, area of allochthonous plantations in close proximity. For both species, a male gametic gene flow of 6–8% was estimated (Unger et al. 2014). In our study, neither Structure nor GeneClass indicated an influence of autochthonous Romanian populations on the juvenile population in R3. According to the forest records, in the neighbourhood of R3, other larch stands originating from unknown (very probably Alpine) sources are located but no autochthonous populations. But if R3 was surrounded by autochthonous larch populations, a much higher introgression rate would be expected at least in the R3 juveniles unless flowering between Romanian and Tyrolian populations is badly synchronised.

5 Conclusions

In this study, we show that populations of European larch from the Alpine and Carpathian regions differ both in genetic diversity and in population genetic structure among and within regions. The anthropogenic transfer of larch from the Alpine region (Tyrol) to the Carpathians was clearly demonstrated. We show that translocated populations as well as introgressed individuals can be identified with relatively little effort. Our data indicate a minimal influence of this introduced material on the native gene pool. Hence, genetic swamping has been probably insignificant and, therefore, the intrinsic value of the gene pool of autochthonous, scattered larch populations as rear-edge populations in the Southern and Eastern Carpathians remains high. However, a general assessment of the intraspecific introgression in this region must remain unanswered, as our samples size was too small, providing only a snapshot of the situation in the region. To fully assess, the impact of translocations on the natural gene pool in a first step maternally inherited markers should be used to help identifying allochthonous populations and in a second step a gene flow study between native and introduced stands should be done. Ideally, this would be supplemented by an assessment of the flowering phenology and other mating system parameters of translocated Tyrolean compared to Romanian larch trees.