Basidiomata of Sebacina epigaea (Berk. & Broome) Bourdot & Galzin and S. incrustans (Pers.) Tul. were collected in Southern Germany from August to October 2010 and 2011 from plots in diverse forest stands that were between 20 km and 250 km apart from each other (Figure 1). Collection areas are categorised as follows: (1) the Bavarian Alps (BA) included five sites separated from each other by a few hundred meters to 20 km. Site BA1 is a pure montane Picea abies plantation with a northern exposure. BA2, BA3 and BA4 sites are dominated by montane mixed forests composed mainly of P. abies. Site BA5 is a Fagus-dominated mixed forest on a lakeside, separated by 20 km from the remaining BA sites. (2) The Western Swabian Alb (WA) included three sites, each separated by 1 km, comprising a submontane slope with mixed forests dominated by Fagus sylvatica with a south-western exposure (WA1) and an eastern exposure (WA3). Site WA2 is a plateau dominated by old specimens of Abies alba. (3) The Eastern Swabian Alb (EA) included two sites (EA1 and EA2) in a planar landscape separated by 1 km. Both sites contain mixed forests dominated by old trees of Fagus sylvatica. (4) The Neckar Valley (NV) included six sites separated by a few hundred meters to 30–10 km. Sites NV1, NV2 and NV3 comprise mixed forest separated from each other by 200 m (NV1 is located on a small riparian lake). Sites NV4 and NV5 correspond to a plantation of Picea abies, and the sampling plots were separated by 2 km. Site NV6 represents a coniferous forest dominated by Pinus sylvestris. (5) Hohenlohe (HO) included two sites separated by 50 m within a mixed forest dominated by Carpinus betulus in a moist basin.
For each collection site, details such as the location, altitude, substrate and ectomycorrhizal plant species nearby were recorded (Additional file 6).
DNA extraction, PCR amplification and sequencing
Total genomic DNA was extracted from dried basidioma fragments using the InnuPREP Plant DNA Kit (Analytik Jena, Jena, Germany) following the standard protocol. Fungal portions contained in Eppendorf tubes were deep-frozen in liquid nitrogen and then ground several times with a sterile plastic pestle.
The internal transcribed spacer (ITS1 and ITS2), 5.8S and D1/D2 regions of the nuclear ribosomal DNA were amplified with the primer combination ITS1F and NL4 (for oligonucleotide primer sequences, see Additional file 8). PCRs were performed in a total volume of 25 μl, containing 5.00 μl GC buffer 5x (including 7.5 mM MgCl2), 15.25 μl water, 1.00 μl dNTP mix (5 mM), 0.50 μl of each primer (25 pmol/μl), 0.25 μl Phusion™ High-Fidelity DNA Polymerase (Finnzymes Oy, Vantaa, Finland) (2 U/μl) and 2.50 μl undiluted DNA. PCR amplification was conducted as follows: 35 cycles of 10 s at 98°C, 20 s at 55°C and 30 s at 72°C, with a final extension of 10 min at 72°C. In the case of negative or weak amplification, PCRs were repeated as follows: 5.00 μl buffer 10x, 0.75 mM MgCl2 (50 mM), 14.50 μl water, 1.00 μl dNTP mix (5 mM), 0.50 μl of each primer (25 pmol/μl), 0.25 μl MangoTaq™ DNA Polymerase (Bioline, Luckenwalde, Germany) (2 U/μl), and undiluted 2.50 μl DNA under the following cycling profiles: 10 cycles of 30 s at 94°C, 45 s at 60°C (−1°C per cycle), and 75 s at 72°C, followed by 26 cycles of 30 s at 94°C, 45 s at 50°C and 75 s at 72°C, with a final extension of 7 min at 72°C.
The second largest subunit of RNA polymerase II (RPB2) regions 5–7 was amplified with the primer combination fRPB2-5F and bRPB2-7.1R using Phusion DNA Polymerase and the PCR conditions indicated above. Subsequently, samples that yielded no or only weak products were amplified with MangoTaq Polymerase and the primer pairs bRPB2-5F and bRPB2-7.1R, or sRPB2-5.1F and sRPB2-7R following the thermocycling program by Matheny .
The mitochondrial adenosine triphosphatase subunit 6 (ATP6) was amplified using the primer combinations ATP6-3/ATP6-4 or sATP6-3/ATP6-4 and MangoTaq Polymerase with the PCR concentration reactions indicated above and the cycling profiles described by Kretzer and Bruns .
The presence of PCR products was monitored by 0.7% agarose gel at 140 V for 15 min in 1× Tris-acetate-EDTA buffer and made visible by ethidium bromide staining and UV light at 254 nm wavelength. The amplified DNA fragments were purified using a ExoSAP-IT® reagent (USB Corporation, Cleveland, OH, USA) diluted 1:20 according to the manufacturer’s instructions. Purified PCR products were cycle sequenced in both directions with a 1:6 BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) on an ABI Prism 3130xl Genetic Analyzer (Applied Biosystems). The primers used for DNA sequencing in both directions are listed in Additional file 8. Fungal vouchers/basidiomata used in this study are deposited in the Herbarium Tubingense (TUB; Additional file 6).
Sequence editing, identity and alignments
Forward and reverse sequence fragments were assembled and edited using Sequencher 4.10.1 (Gene Codes Corporation, Ann Arbor, MI, USA). The ITS + 5.8S + D1/D2 and RPB2 sequences were screened against those available in the GenBank database (http://www.ncbi.nlm.nih.gov) using the Basic Local Alignment Search Tool (BLAST; ). ATP6 sequences were analysed with Basidiomycota sequences from our lab (unpublished data by SG) because there are no available sequences in GenBank for Sebacinales.
In our study, phylogenetic and demographic population analyses were inferred from morphological species concepts. We assembled seven datasets: datasets 1 and 2, containing 63 Sebacina epigaea and 63 S. incrustans ITS + 5.8S + D1/D2 sequences which were automatically aligned using MAFFT 6.815b under the E-INS-i algorithm  and POA 2 ; datasets 3 and 4, including 85 (S. epigaea) and 78 (S. incrustans) RPB2 (exon 6) sequences aligned with MAFFT; datasets 5 and 6, comprising 11 (S. epigaea) and 42 (S. incrustans) ATP6 sequences aligned using MAFFT; and dataset 7, containing our own ITS haplotype sequences of both species, and those from GenBank and UNITE . Datasets containing sequences obtained from Sebacina basidiomata and ectomycorrhizae from Europe were aligned using MAFFT und POA. For datasets 1, 2 and 7, selection of the most consistent alignments (MAFFT) was performed using trimAl 1.4 . Nucleotide alignments from datasets 3 to 6 were improved manually from amino acid codon sequences in Se-Al 2.0a11 Carbon .
The DNA sequences (original datasets including heterozygous sequences) used in this study have been submitted to GenBank (accession numbers JQ665465–JQ665564 for ITS, 5.8S and D1/D2; JQ665565–JQ665657 for RPB2 exon 6; and JQ665658–JQ665710 for ATP6). Alignments are deposited at TreeBASE under submission ID S13159.
Haplotype reconstruction, detecting recombination and selection pressure
Because all sequences are derived from dikaryotic (n + n) isolates, we used PHASE 2.1 [33, 34] as implemented in DnaSP 5.10.01  using the Markov Chain Monte Carlo option followed by 1000 iterations under a hybrid model to infer haplotype phases. From haplotype sequences within each species, we removed indels and excluded infinite-site violations using Map as implemented in SNAP Workbench 2.0 [36, 37]. Population recombination parameters were estimated following the method of Hudson and Kaplan  based on the minimum number of recombination events in a sample (R
M) using DnaSP and RecMin  as implemented in SNAP Workbench. The significance of the RM estimation was calculated by performing 10,000 coalescent simulations  in DnaSP. Estimations of the selection pressure on coding sequences were based on the ω = dN/dS ratio by comparing the rates of non-synonymous and synonymous mutations. We calculated ω ratios using the Synonymous Non-synonymous Analysis Program (SNAP) .
Phylogenetic relationships, congruence among gene phylogenies, neutrality tests and demographic structure analyses
Maximum likelihood analysis with combined rapid bootstrapping under the GTRCAT model was computed from 1000 runs with RAxML 7.0.4 [42, 43]. The phylogenetic trees were midpoint rooted using FigTree 1.3.1.
Conflicting phylogenetic signals of the different datasets were checked using a partition homogeneity test/incongruence length difference (ILD) test  as implemented in PAUP* . The number of ILD replicates was set to 1000, setting one tree per replicate and branch swapping with tree bisection–reconnection (TBR). In order to test if the topologies of the different analyses and datasets were significantly different we used the maximum parsimony based Shimodaira-Hasegawa (SH) test as implemented in PAUP, using 1000 replicates with TBR swapping.
DnaSP was used to calculate the standard indices of population diversity for each lineage and sampling area. To detect departures from a constant population size under the neutral model, Tajima’s D , and Fu’s and Li’s D* and F*  statistics were calculated using Arlequin 22.214.171.124 . The significance of these values was obtained in neutrality tests with 1000 permutations.
Intraspecific gene genealogies were inferred using the median joining  and statistical parsimony networks. Haplotype genealogies were constructed using Network 4.6.1 (http://www.fluxus-engineering.com) with the parameter ϵ = 10 and the ‘MP’ option to identify and eliminate unnecessary median vectors and links . Network graphics were generated using Network Publisher 1.3 (http://www.fluxus-engineering.com). Analyses of genetic differentiation among species using a 95% parsimony limit reconstruction criterion  suggest that biological species often form unconnected parsimony networks. Based on these measurements, we reconstructed a parsimony network of haplotypes with a 95% connection probability limit using TCS 1.21 . These analyses were run separately for sequences with and without recombination blocks. To examine the amount of genetic structure within and among population groupings, hierarchical analyses of molecular variance (AMOVA; [54–56]) were conducted in Arlequin. For this purpose, we tested two hypothetical scenarios: (i) differentiation explained by geographic distribution: North (HO) vs. Central (WA, EA, NV) vs. South (BA) (scenario 1); and (ii) differentiation between populations separated by a geographic barrier. For this, we compared populations from the north of the Swabian Alb (WA, EA, NV, HO) vs. those south of the Swabian Alb (BA) (scenario 2). To search for genetic differentiation between population samples and major lineages, exact tests of population differentiation  and FST values were computed in Arlequin. The significance of these estimators was assessed using 1000 permutations. Furthermore, we examined whether a correlation existed between genetic (FST) and geographic distance matrices with a Mantel test using zt 1.1 .
The microscopic structures of the basidiomata were analysed with a light microscope using standard procedures. Microscopic structures were analysed from dried specimens using 3% KOH and stained with aqueous 1% phloxine. Measurements of basidiospores (n = 51) are given as length and width, uncommon extreme values appear in brackets. Basidiospore length/width quotient (Q), mean and standard deviation were calculated.