Genetic patterns in Pinus nigra from the central Balkans inferred from plastid and mitochondrial data

Pinus nigra J.F. Arnold, European black pine, is a typical component of Mediterranean and subMediterranean coniferous forests with highly fragmentary distribution. Western Mediterranean populations of this species have been studied genetically to date, while eastern populations from the central Balkans, which are larger and more abundant, are still genetically understudied. We analyzed seven populations of P. nigra representing all infraspecific taxa recognized within the central Balkans (subspecies nigra with varieties nigra and gocensis Đorđević; and subspecies pallasiana (Lamb.) Holmboe with varieties pallasiana and banatica (Endl.) Georgescu et Ionescu), with three chloroplast microsatellites (cpDNA SSRs) and one mitochondrial (mtDNA) locus. Although our molecular data failed to support circumscription of studied infraspecific taxa, we found that genetic patterns at both genomes are in accordance with those found previously in westward populations of this species, that is – exceptionally high levels of genetic diversity (HT = 0.949) and low genetic differentiation (GST = 0.024) at the cpDNA level, and moderate levels of genetic diversity (HT = 0.357) and genetic differentiation (GST = 0.358) at the mtDNA level. Based on genealogical relations of mtDNA types currently present in Balkans’ and Iberian/ African populations, we inferred that the ancestral gene pool of P. nigra already harbored polymorphism at position 328 prior to the divergence to two lineages currently present in westward and eastward parts of the species range distribution. Subsequent occurrence of three mutations, which distinguish these two lineages, suggests their long-term isolation.


Introduction
Pinus nigra J.F. Arnold, European black pine, is one of the economically and ecologically most important pine species in Europe. It is a typical component of Mediterranean and sub-Mediterranean coniferous forests, with highly fragmentary distribution extending from the North Africa, the northern Mediterranean Basin and eastwards to the Black Sea (Gaussen et al. 1993). Due to its Tertiary origin, relict nature and a prominent phenotypic plasticity throughout the range distribution, its taxonomic status and infraspecific classification are still somewhat unclear (e.g. Vidaković 1991;Gaussen et al. 1993;Christensen 1993; see also references in Naydenov et al. 2006 andGanopoulos et al. 2013).
To date, fragmentary distributed natural populations of P. nigra from the western part of the species range have been studied genetically (e.g. Afzal-Rafii and Dodd 2007;Jaramillo-Correa et al. 2010;Soto et al. 2010). However, those from the Balkans and Turkey, which are larger and more abundant, have occasionally been studied from various aspects (e.g. Liber et al. 2003;Lučić et al. 2010;Gülsoy et al. 2014) but rarely from a population genetics aspect (e.g. Naydenov et al. 2006Naydenov et al. , 2015Šarac et al. 2015). The genetic knowledge on these populations is required not only for shedding new light on species infraspecific circumscription and overall genetic and/or phylogeographic structure, but also for implementation of suitable conservation measures and future breeding efforts (Frankham et al. 2002). Studies of this kind, however, are not straight-forward because, on the one hand, P. nigra is one of the most commonly planted conifers in Europe and the origin of the planted material cannot be tracked in all cases, and, on the other hand, this predominantly outcrossing species is characterized by a rather high dispersal ability (pollination and seed dispersal by wind, Afzal-Rafii and Dodd 2007;Jaramillo-Correa et al. 2010;Soto et al. 2010) suggesting possible introgression of alien gene pools into native populations.
Within the central Balkans, two subspecies of P. nigra (ssp. nigra and ssp. pallasiana (Lamb.) Holmboe), and their varieties (var. nigra and gocensis Đorđević of the former subspecies, and var. pallasiana and banatica (Endl.) Georgescu et Ionescu of the latter subspecies), have been recognized (Jovanović 1992). It has been shown recently that these infraspecific taxa differ at the biochemical level (Bojović et al. 2012;Šarac et al. 2013). In the frame of the present study, we employ: 1) paternally inherited and potentially more variable chloroplast microsatellites (cpDNA SSRs) shown to be informative for assessing genetic patterns in P. nigra and other pines (Naydenov et al. 2006;Afzal-Rafii and Dodd 2007;Jaramillo-Correa et al. 2010;Soto et al. 2010), and 2) a more slowly evolving maternally inherited mitochondrial (mtDNA) locus nad7 intron 1, proven to be suitable for depicting phylogeographic structure in westward populations of P. nigra (Jaramillo-Correa et al. 2010), and test them in individuals of P. nigra from seven populations representing all recognized infraspecific taxa of this species from the central Balkans in order to: 1) assess whether cpDNA and/ or mtDNA markers are suitable for distinguishing infraspecific taxa of P. nigra from the central Balkans; 2) estimate levels of genetic diversity and genetic structuring at cpDNA and mtDNA levels in these genetically understudied populations of P. nigra; and 3) use available variation at the studied mtDNA locus to make inferences on historical events that have led to the current genetic patterns in westward and eastward parts of the species range.

Plant material and DNA extraction
Samples (young twigs with needles) were collected in early fall 2009 from 14 -15 trees per population from seven native populations of P. nigra from the central Balkans which represent all infraspecific taxa of this species (ssp. nigra represented by var. nigra and var. gocensis; and ssp. pallasiana represented by var. pallasiana and var. banatica) recognized by the Flora of Serbia (Jovanović 1992). We sampled 104 trees in total (Table 1). Trees sampled from each population were distant at least 30 m from each other, and sampled populations were not in the vicinity of known planted trees/populations of this species. Plant material, deposited in labeled polyethylene bags, was stored at -20 °C prior to DNA extraction.
Approximately 30 mg of the plant tissue of each individual, homogenized with TissueLyser II (Qiagen, Valencia, CA, USA), was used for DNA extraction following Aleksić et al. (2012).
Three mitochondrial (mtDNA) loci, the second intron of the NADH dehydrogenase subunit 1 gene -nad1 gene intron 2 (primers published by Demesure et al. 1995), nad5 intron 4 (Dumolin-Lapegue et al. 1997), and nad7 intron 1 (Jaramillo- Correa et al. 2004;Tian et al. 2010), were initially tested in the same panel of seven P. nigra individuals following PCR conditions provided in the corresponding publications. However, PCR products were obtained only at locus nad7 intron 1 using the F primer of Tian et al. (2010) and the R primer of Jaramillo-Correa et al. (2004), and the following PCR conditions: a total volume of 25 μl containing 25 ng template DNA, 1 × Taq buffer with (NH 4 ) 2 SO 4 (Fermentas UAB, Vilnius, Lithuania), 2.5 mM MgCl 2 (Fermentas UAB, Vilnius, Lithuania), 0.2 mM dNTPs, 0.1 μM of each forward and reverse primers, 0.80% BSA (Bovine Serum Albumin, Fermentas UAB, Vilnius, Lithuania), and 0.025 U Taq DNA polymerase (Fermentas UAB, Vilnius, Lithuania), and a touchdown amplification protocol: 4 min initial denaturation at 94 °C, 11 cycles of 1 min denaturation at 94 °C, 45 s annealing at 68 °C, 2 min extension at 72 °C with progressive decrease of the annealing temperature by 0.5 °C/cycle, followed by 26 cycles of 1 min denaturation at 94 °C, 45 s annealing at 56 °C, 2 min extension at 72 °C, and final extension of 10 min at 72 °C. PCR products were Sanger sequenced by Macrogen Europe using F or both (F and R) primers (see Results).

Data analysis
At the cpDNA level, we assessed the following parameters of genetic diversity at the population level and in overall sample: the number of alleles (A), the number of private alleles (PA), the number of chloroplast haplotypes (ch), the number of private haplotypes (cph), unbiased within population gene diversity (H S , equivalent to H E for diploid data, Weir, 1996), and the total gene diversity (H T ) using CONTRIB 1.02 (Petit et al. 1998). Using the same software, we assessed also the number of mitochondrial haplotypes (mh), H S and H T at the mtDNA level in populations and at the species level. The standardized coefficients of genetic differentiation among populations, G ST , at the cpDNA and mtDNA levels were also calculated in CONTRIB. Pairwise population differentiation estimates, F ST , based on cpDNA haplotype frequencies, and the analysis of molecular variance (AMOVA, Excoffier et al. 1992), based on the sum of squared number of repeat difference between cpDNA haplotypes (Slatkin 1995), were calculated using ARLEQUIN 3.11 (Excoffier et al. 2005). The significance of tests was examined by 15 000 permutations.
TESS 1.2 (François et al. 2006;Chen et al. 2007) was used to further study genetic differentiation of populations based on cpDNA data, by estimating the number of genetically similar clusters. TESS incorporates the spatial locations of all individuals into the analyses to assess genetic cluster membership. Spatial coordinates of individual trees within defined areas, corresponding to population sizes, were generated using population sample coordinates. Using no admixture model and interaction parameter of 0.6, 100 runs at each K ranging from 2 to 4 with 50 000 sweeps and a burn-in of 10 000 for each run were performed. At each K, 20 runs with the lowest deviance information criterion (DIC) were recorded and used for calculating average DICs. As recommended by Durand et al. (2009), the optimal number of clusters may be determined by plotting average DIC against K, and depicting K at the beginning of the plateau.

Results
Three cpDNA SSRs were successfully amplified in 104 P. nigra individuals, yielding 5-9 alleles per locus, 22 alleles in total, and 38 haplotypes ( Table 1). The number of PA and cph were high (7 and 22, respectively) as well as average H S (0.926) and H T (0.949). Genetic differentiation of populations was low (G ST = 0.024). MtDNA locus nad7 intron 1 was successfully amplified and sequenced with F primer in 99 P. nigra individuals while sequencing with R primer was successful in 73 individuals. The distribution of mtDNA haplotypes among populations was not random, and two detected haplotypes (T and G, see later) were found at more or less equal frequencies in populations of P. nigra ssp. nigra var. nigra (I, II and III), while populations IV to VII, representing the other studied taxa, were fixed for haplotype T. Thus, the average H S and H T were low to moderate (0.229 and 0.357, respectively), while G ST was moderate (0.358). Populations I, II and III represent P. nigra ssp. nigra var. nigra, populations IV and V represent P. nigra ssp. nigra var. gocensis, population VI represents P. nigra ssp. pallasiana var. pallasiana, and population VII represents P. nigra ssp. pallasiana var. banatica (see also Table 1); the size of circles in haplotype network rooted with P. sylvestris (B) and unrooted haplotype network (C) is proportional to the relative frequency of haplotypes T and G in eastward populations (see Table 1), and haplotypes h1 to h4 in westward populations, but the relative frequency of haplotypes in these regions is not comparable; mtDNA haplotype labels are given in the legend of Table 2; D -clustering of populations at K = 2 and K = 3 using TESS.
All population pairwise F ST values were insignificant, and ranged from 0.066 (populations IV and V representing P. nigra ssp. nigra var. gocensis) to 0.000 or were negative (data not shown). Although AMOVA revealed that 1.72% and 98.28% of overall molecular variation were attributed to among and within population variation, respectively, these values were insignificant. When two subspecies of P. nigra were compared (populations I to V representing ssp. nigra vs. populations VI and VII of ssp. pallasiana), somewhat increased among population variation (2.06%), decreased within population variation (97.20%) and 0.73% of variation among populations within groups were observed. All values, however, were insignificant. The lack of genetic structuring of populations at the cpDNA level was revealed by TESS analysis as well, because at each tested K, all individuals were strongly assigned to a single gene pool, with a few individuals from populations VI and VII assigned to alternative gene pools (Fig. 1d). Therefore, calculations required for assessing the optimal number of clusters were not performed because individuals from all populations belong to a single gene pool.
The length of aligned matrices comprising nad7 intron 1 sequences of P. nigra generated in the present study, and those reported by Jaramillo-Correa et al. (2010), as well as sequences of P. sylvestris (Naydenov et al. 2007), were 985 bp, 985 bp and 991 bp, respectively. The positions of variable sites (point and length mutations) in the latter matrix are shown in Table  2. P. nigra mtDNA haplotypes from the Balkans harbored a single polymorphic site at position 328 (G/T transversion), yielding two mtDNA haplotypes: T (harboring T at position 328), and G (G at position 328), which were deposited in the GenBank (accessions KT343352 and KT343353). The Balkans' and Iberian/African mtDNA haplotypes of P. nigra shared variability at position 328, and differed in three mutations, at positions 435, 472 and 520 (Table 2). Nucleotide characters G and A at positions 435 and 520, respectively, are synapomorphic for individuals from Iberia/Africa, while a single-bp insertion (C) at position 472 is synapomorphic for individuals from the central Balkans. All P. sylvestris sequences harbored T at position 328, and differed from P. nigra sequences in 8 mutations, at positions 334,337,396,445,565,775,889, and 931 (see Table 2).
Haplotype networks constructed with and without outgroup, P. sylvestris, are shown in Fig. 1b and 1c, respectively.

Discussion
High levels of cpDNA SSR diversity (H T = 0.949) and low genetic differentiation of populations (G ST = 0.024) found in populations of P. nigra from the central Balkans are in the range of those detected in westward populations of this species (Afzal-Rafii and Dodd 2007;Jaramillo-Correa et al. 2010;Soto et al. 2010). These findings, along with upper Miocene fossil records of P. nigra in central and southern Serbia (Petković 1987), would suggest a long persistence of this species within the Balkans, supported also by findings of Naydenov et al. (2006). Interestingly, these authors found somewhat lower levels of cpDNA diversity and increased genetic differentiation of P. nigra populations in Bulgaria, with four groups of populations whose distinctiveness was explained by historical events, specific topography of the region and anthropogenic impact. The alternative clustering of these populations found at the nuclear genome, and discordance of cpDNA and nuclear data, was attributed to different migration patterns and effective radius of dispersal of P. nigra pollen and seed (Naydenov et al. 2015). On the other hand, we found moderately high levels of mtDNA diversity in eastward populations of P. nigra (H T = 0.357), and a non-random distribution of two mtDNA haplotypes, which has led to the moderate levels of genetic differentiation (G ST = 0.358). Since the observed levels of genetic differentiation at the mtDNA level are T and G -mtDNA haplotypes found in P. nigra individuals from the central Balkans; h1, h2, h3 and h4 -Iberian (h1 -h3) and north African (h4) haplotypes of P. nigra reported by Jaramillo-Correa et al.
(2010); PS1, PS2, and PS3 -P. sylvestris haplotypes available in GenBank (accessions DQ665913, DQ665914, and DQ665915, Naydenov et al. 2007); Δ -The positions of mutations which may be used for distinguishing westward and eastward lineages of P. nigra (nucleotide characters "G" and "A" at positions 435 and 520 are synapomorphic for P. nigra individuals from the Iberian Peninsula/northern Africa; one bp insertion (C) at position 472 is synapomorphic for individuals from the central Balkans) usually expected at much larger studied areas than ours (e.g. Jaramillo-Correa et al. 2010;Tian et al. 2010; Geburek 2014 and references therein), we assumed that the Tara Mountain, with P. nigra populations characterized by increased levels of mtDNA diversity, may represent an interglacial refugium of this species, similarly as observed in another cold-adapted species, Picea omorika (Panč.) Purk. (Aleksić and Geburek 2014).
Although it has been shown recently that cpDNA variation in closely related spruces is more species-specific than mtDNA variation, and that the latter is often shared among related conifer species (Du et al. 2009), we were not able to use cpDNA nor mtDNA polymorphisms to support infraspecific circumscription of studied P. nigra taxa which has recently been implied at the biochemical level (Bojović et al. 2012;Šarac et al. 2013). The sole conclusions with regard to this aspect of our study are that at the cpDNA level, two studied populations of P. nigra ssp. nigra var. gocensis display the highest but insignificant genetic differentiation, and that populations VI and VII, representing varieties pallasiana and banatica of P. nigra ssp. pallasiana, respectively, may harbor genetically distinct individuals. It is worth mentioning that these populations are at the range-edge of P. nigra ssp. pallasiana distribution (Vidaković 1991) and thus, may be genetically different from populations of these taxa found within their main range (e.g. Bjedov et al. 2015 and references therein). At the mtDNA level, populations of P. nigra ssp. nigra var. nigra are genetically diverse while populations of the other studied infraspecific taxa are invariable. Nonetheless, mtDNA variation enabled new insights into historical events that have led to the current distribution of mtDNA diversity in westward and eastward parts of the P. nigra range.
Based on the haplotype network rooted with P. sylvestris (Fig. 1b), we may assume that P. nigra mtDNA haplotypes harboring G at position 328 (Table 2), which are currently less abundant in eastward but prevalent in westward populations (h3 and h4), are derived. This would imply that the ancestral gene pool of P. nigra, which initially spread throughout the Mediterranean Basin, harbored T at position 328, and that the same mutation (T → G at position 328) must have had occur twice, i.e. in the ancestral P. nigra gene pool within the Balkans, and in the ancestral gene pool of this species within the western Mediterranean Basin. Alternatively, due to the long persistence of P. nigra in the Mediterranean (Naydenov et al. 2006(Naydenov et al. , 2015Jaramillo-Correa et al. 2012) and its high dispersal ability, it is possible that the above mentioned mutation occurred at one site, and expanded later on to the other site. However, the more plausible and parsimonious scenario is that the polymorphism at position 328 was already present in the ancestral gene pool of P. nigra which inhabited the Mediterranean Basin, and that three mutations, which distinguish two lineages of this species confined to the Balkans and the Iberian Peninsula/northern Africa, occurred after their divergence. It is worth mentioning the nucleotide characters G and A at positions 435 and 520, respectively, and one bp insertion (C) at position 472 may be used as molecular diagnostic characters for distinguishing westward and eastward P. nigra lineages. These findings altogether further support long persistence of P. nigra in southern Europe, and long-term isolation of two lineages of this species detected in the present study.