Introduction
Cladonia is the largest genus in the lichen-forming fungal family Cladoniaceae, consisting of c. 459 accepted species (T. Ahti, August 2015, unpublished data). Based on morphology and secondary chemistry, Ahti (Reference Ahti2000) divided the neotropical species of Cladonia into seven taxonomic sections and three sections were recognized in the segregate genus Cladina (Ahti Reference Ahti2000). The division was applicable to most species in the world. The group called Cladina (known as reindeer lichens) is most abundant in the coniferous belt of the Northern Hemisphere and in the Nothofagus regions in the Southern Hemisphere, but is also known in sandy areas of the south-eastern United States and elsewhere, as well as at high altitudes in many mountain ranges. While the lack of competitive ability of lichens with plants is well known, the reindeer lichens have adapted better than almost all other lichens to the terrestrial niches uninhabited by vascular plants and bryophytes. Some species, such as Cladonia arbuscula, C. rangiferina, C. stygia and C. stellaris, are important components of northern ecosystems where they provide vast areas of ground cover (Auclair & Rencz Reference Auclair and Rencz1982; Shaver & Chapin Reference Shaver and Chapin1991) and form a major component of the winter food for caribou and reindeer (Svihus & Holand Reference Svihus and Holand2000; den Herder et al. Reference den Herder, Kytöviita and Niemelä2003). Knowledge of their species status would inform ecosystem management and maintenance of biodiversity.
In recent times the genus Cladina has not been recognized by most authors because in phylogenetic studies it is nested within Cladonia and is not monophyletic (Stenroos et al. Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002, Reference Stenroos, Pino-Bodas, Weckman and Ahti2015; Guo & Kashiwadani Reference Guo and Kashiwadani2004). The group of reindeer lichens constitute Cladonia, sections Crustaceae (=Cladina section Cladina), Tenues, and Impexae, which typically have highly branched ecorticate podetia and a soon disappearing crustose primary thallus. Section Tenues typically has anisotomic branching with dichotomous branches, section Crustaceae has anisotomic branching with two, three or four divisions to the branches, and section Impexae has isotomic or subisotomic branching with two, three, four, or five divisions to the branches. Phylogenetic relationships among all sections of Cladonia have been examined previously (Stenroos et al. Reference Stenroos, Ahti and Hyvönen1997, Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002), but usually with only one to a few specimens of each species included in the analysis. While Stenroos et al. (Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002) recommended studying a larger number of specimens for each species, knowledge is lacking of the number of monophyletic groups, the species composition of each monophyletic group, and their phylogenetic relationships. Earlier studies were focused on species identification (Ruoss Reference Ruoss1987a , Reference Ruoss b ; Ruoss & Ahti Reference Ruoss and Ahti1989) and recent studies show relationships among distant geographical collections of Cladonia arbuscula in the broad sense (Myllys et al. Reference Myllys, Stenroos, Thell and Ahti2003; Piercey-Normore et al. Reference Piercey-Normore, Ahti and Goward2010), but knowledge concerning phylogenetic relationships is insufficient. Monophyly and diagnosability are considered to be important criteria for species delimitation (e.g. Bacon et al. Reference Bacon, McKenna, Simmons and Wagner2012). While monophyly may suggest low levels of intraspecies variation, the absence of monophyly may suggest ongoing speciation, incomplete lineage sorting, or interbreeding among species.
The goals of this study were: 1) to examine the monophyly of the sections Crustaceae, Tenues and Impexae; 2) to reconstruct the phylogenetic relationships among selected species of reindeer lichens in the genus Cladonia; and 3) to examine species delimitations.
Materials and Methods
Lichen specimens were collected from Canada or borrowed from herbaria, and additional sequences were obtained from the NCBI GenBank (Table 1). Sixty-two representative specimens are deposited in the University of Manitoba Herbarium (WIN) or the Botanical Museum, Finnish Museum of Natural History, Helsinki (H) (Table 1).
Table 1 Collection location, collection numbers, and accession numbers for the Cladonia specimens used in this study. Specimens with collection numbers were used to generate sequences and those with references were obtained from GenBank
For the phylogenetic analysis, either spore cultures or pieces of dry thalli (10–20 mg) from the apical region of each lichen sample were selected and visually inspected for contaminating debris. The DNA was isolated using a modified CTAB (hexadecytrimethylammonium bromide) protocol (Grube et al. Reference Grube, DePriest, Gargas and Hafellner1995). The polymerase chain reaction (PCR) of fungal DNA on the internal transcribed spacer 1 and 2 (ITS1 and ITS2) and the 5.8S of the nuclear ribosomal DNA (rDNA) was performed using the primers 1780F-5’ (Piercey-Normore & DePriest Reference Piercey-Normore and DePriest2001) and ITS2KL-3’ (Lohtander et al. Reference Lohtander, Myllys, Sundin, Källersjö and Tehler1998), and on the mitochondrial small subunit (mtSSU) DNA using the primers mrSSU2 and mrSSU3R (Zoller et al. Reference Zoller, Scheidegger and Sperisen1999). Where there were problems in amplifying across both ITS regions for some samples, the primers ITS1F, ITS2, ITS3 and ITS4 (White et al. Reference White, Bruns, Lee and Taylor1990) were used to amplify the ITS region in two fragments. PCR reaction mixtures (20 µl) contained 20 ng of template DNA, 1× PCR buffer (50 mM KCl, 20 mM Tris), 0·5 µM of each forward and reverse primer, 3·0 mM of MgCl2 (2·0 mM MgCl2 for mtSSU), 200 mM of each dNTP (Invitrogen Life Technologies, California, USA), and 0·1 U of Taq polymerase (Invitrogen Life Technologies, California, USA). Amplification was carried out in a Biometra® TGradient thermocycler (American Laboratory Trading Inc., Connecticut, USA). The PCR conditions were as follows: initial denaturation at 95 °C for 5 min; 30 cycles of denaturation at 95 °C for 1 min, annealing at 54 °C for 1 min, and an extension at 72 °C for 1 min 30 s for all primers. For samples for which we had difficulties with PCR amplification, touchdown cycles were as follows: initial denaturation at 95 °C for 5 min; 1 cycle of denaturation at 95°C for 1 min, annealing at 60 °C for 1 min, and an extension at 72 °C for 1 min 30 s. Annealing temperature of the following 3 cycles was dropped by 2 °C at each cycle (58, 56, 54) followed by 26 cycles with an annealing temperature of 52 °C.
Four to six identical 50 μl reaction volumes of PCR product were pooled for DNA sequencing and gel was purified using the Wizard® SV Gel and PCR Clean-Up System (Promega, Wisconsin, USA) following the manufacturer’s instructions. Cycle sequencing reaction volumes were 20 μl, containing 60–70 ng of purified DNA, BigDye V3.1 (Applied Biosystem, California, USA) and the same PCR primers that were used for sequencing. Reactions were cleaned using the ethanol/EDTA precipitation method (Applied Biosystem Handbook) according to the manufacturer’s instructions. The dried product was dissolved in 20 μl formamide and loaded into a 96-well plate for sequencing on a 3130 Genetic Analyzer (Applied Biosystems, California, USA). The sequences were edited using Sequencher® version 4.6 (Gene Codes Corporation, Michigan, USA). In addition, 57 accessioned DNA sequences were retrieved from NCBI GenBank and were included in the phylogenetic analysis. All sequences were automatically aligned using the ClustalX (Jeanmougin et al. Reference Jeanmougin, Thompson, Gouy, Higgins and Gibson1998) program and manually edited.
To infer the relationships of the sections Crustaceae, Tenues and Impexae within Cladonia we constructed a matrix with 69 sequences (belonging to 51 species) of ITS rDNA. In this alignment, one sequence per major clade of the phylogeny presented by Stenroos et al. (Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002) was included, along with one or two sequences per species of Crustaceae, Tenues and Impexae. The type species of each section were included. Cladonia wainioi was assigned as the outgroup taxon because of its basal position in the phylogenetic tree of the genus Cladonia (Stenroos et al. Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002). The ambiguous positions were identified and removed by Gblocks version 0.91b (Castresana Reference Castresana2000) using the less stringent options. Three separate ITS rDNA alignments were then constructed, one per monophyletic Cladina group identified in the previous analysis including multiple specimens per species. All nucleotide sequences generated in this study have been deposited in the NCBI GenBank (Table 1).
The alignments were subjected to maximum parsimony (MP), maximum likelihood (ML) and Bayesian analyses. Maximum parsimony analyses were run using PAUP* 4.0b10 (Swofford Reference Swofford2003), ML using RAxML 7.0.3 (Stamatakis Reference Stamatakis2006) and Bayesian analyses using MrBayes v3.2 (Ronquist & Huelsenbeck Reference Ronquist and Huelsenbeck2003; Ronquist et al. Reference Ronquist, Teslenko, van der Mark, Ayres, Darling, Höhna, Larget, Liu, Suchard and Huelsenbeck2012). MP analyses were performed using tree bisection and reconnection branch swapping, heuristic searches with 1000 random addition replicates and bootstrap searches of 1000 resamplings (Felsenstein Reference Felsenstein1985) using the heuristic option. ML analyses were performed using the GTRGAMMA model and fast boostrap searches of 500 replicates. The best models of nucleotide substitution were estimated using jModelTest 2.1.1 (Darriba et al. Reference Darriba, Taboada, Doallo and Posada2012). These models were applied for the Bayesian analyses (Table 2). The Bayesian analysis used 10 000 000 generations with two parallel runs performed. In each run, four chains were used and initiated with a random tree. The convergence between runs was diagnosed by standard deviation of split frequencies <0·01. At the end of the runs, in all cases the split frequencies fell below 0·01. TRACER 1.5 (Rambaut & Drummond Reference Rambaut and Drummond2007) was used to plot the log-likehood scores of sample points against generation time. Afterwards, the first 2 500 000 generations were removed and the 50% majority-rule consensus tree was calculated. Bootstrap values greater than 70% and posterior probabilities greater than 0·95 are reported in the phylogenies.
Table 2 Data from the alignments and analyses of the different clades of Cladonia
n=number of sequences including the outgroup; PI=parsimony-informative; CI=consistency index; RI=retention index; MP=maximum parsimony; –Lnl=likelihood values; ML=maximum likelihood.
Haplotype analysis was performed using 34 mtSSU sequences belonging to 14 species and one mtSSU sequence obtained from GenBank. A haplotype network analysis is recommended when polymorphism is low (Clement et al. Reference Clement, Posada and Crandall2000). The TCS program version 1.21 (Clement et al. Reference Clement, Posada and Crandall2000) was used to construct a parsimony network. The parsimony probability criterion (Templeton et al. Reference Templeton, Crandall and Sing1992) with gaps coded as fifth character state and a 95% parsimony threshold for network relationships was used in the analysis. Chromatograms were examined to rule out miscalled bases. Haplotype analyses using the same procedure were also conducted on three pairs of species (C. rangiferina-C. stygia; C. arbuscula-C. mitis; and C. portentosa-C. terrae-novae) using the ITS rDNA alignments. An analysis of molecular variance (AMOVA) was conducted in GenAlEx ver. 6.5 (Peakall & Smouse Reference Peakall and Smouse2012) with 999 permutations to determine the extent of shared polymorphism among species duplets. The same ITS rDNA alignments used for the haplotype networks were also used for AMOVA analysis. AMOVA populations were defined as species duplets and the species in each of the duplets were defined as subpopulations.
Results
Forty-four ITS rDNA sequences generated in this study and additional sequences from GenBank were used in the phylogenetic analyses. The alignment constructed to study the phylogenetic relationships of Crustaceae, Tenues and Impexae contained 537 unambiguously aligned sites, of which 190 were parsimony-informative. The ITS rDNA maximum parsimony analysis generated 1000 equally parsimonious trees with 847 changes. The consistency index (CI) and retention index (RI) were 0·4864 and 0·7540, respectively. The ML and Bayesian analyses (using GTRGAMMA model) produced a tree with a topology that was consistent with that of the MP tree and with likelihood values of −Lnl=5167·301 and 5248·02, respectively. For the three other constructed phylogenies, Table 2 summarizes the data used in the analyses and models selected by jmodeltest.
The phylogeny shows a number of clades, among which four are notable: A, B, C, and D. Clades A and B represent section Impexae, are highly supported, and are basal to the taxa in the tree (Fig. 1). Clade C represents some species of section Crustaceae (C. arbuscula, C. densissima, C. mitis and C. submitis) with support values of 69% MP bootstrap (BS), 74% ML BS and 0·99 posterior probability (PP). Most of clade D represents species of sections Crustaceae and Tenues together, and is moderately supported (81% MP BS, 79% ML BS and 0·95 PP). Clades C and D form a polytomy with the Unciales.
Fig. 1 Phylogenetic tree generated from ML analysis based on ITS rDNA sequences. It shows the placement of the sections Crustaceae, Tenues, and Impexae in the genus Cladonia. The taxonomic sections are represented as A, B, C, and D. The values on the branches are ≥70% bootstrap for MP and ML analyses and ≥0·95 posterior probability for Bayesian analysis. The thick branches represent branches supported in the three analyses.
The analysis of clades A and B with additional ITS rDNA sequences yielded some highly supported subclades (Fig. 2A). Clade A contains C. stellaris, which is sister to the only specimen of C. delavayi. However, C. delavayi (Himalayan) and C. stellaris are very different from one another in morphology and chemistry, as well as their ITS rDNA sequences. Clade B contains three well-supported subclades. Subclade B1 represents C. portentosa s. str. (including the recently synonymized species C. azorica and C. macaronesica according to Pino-Bodas et al. Reference Pino-Bodas, Pérez-Várgas, Stenroos, Ahti and Burgaz2016) and C. terrae-novae. Neither C. terrae-novae, C. portentosa s. str., nor C. portentosa subsp. pacifica could be supported as monophyletic groups. Cladonia evansii alone forms subclade B2 (97% MP BS, 95% ML BS and 1·0 PP). The subclade B3 is represented by C. pycnoclada. Cladonia confusa is not monophyletic.
Fig. 2 Phylogenetic trees from ML analyses based on ITS rDNA sequences showing a wide sampling of sequences per species. A, analysis from clades A and B; B, analysis from clade C. Numbers on the branches indicate bootstrap support ≥70% in MP and ML analyses and posterior probabilities ≥0·95 from the Bayesian analysis. The subclades that are supported in the analyses are numbered in each tree.
The analysis of clade C, representing C. arbuscula, C. submitis, C. densissima and C. mitis, shows two unsupported subclades, one with all the sequences of C. mitis and the other with most of the C. arbuscula sequences (Fig. 2B). It is remarkable that other chemically and morphologically distinct segregates, namely C. submitis and C. densissima, are also nested within C. arbuscula. The analysis thus suggests that even the rest of C. arbuscula may not be taxonomically uniform. One subclade containing two specimens of C. mitis (subclade C1) is supported. Two other subclades with specimens of C. arbuscula (subclades C2 and C3) are also supported.
The analysis of clade D with additional sequences shows several supported subclades (Fig. 3). Section Tenues is represented by two monophyletic groups within subclades D1 and D2; subclade D1 contains six sequences from C. ciliata (82% MP BS, 92% ML BS and 1·0 PP). However, there is no support for the separation of C. ciliata var. tenuis and var. ciliata. Subclade D2 contains eight sequences from C. subtenuis (only 54% MP BS, 84% ML BS and 1·0 PP). Cladonia rangiferina and C. stygia are not monophyletic. Subclade D3 consists of specimens representing C. argentea, C. dendroides, C. rangiferina, and C. rotundata (87% MP BS, 93% ML BS and 1·0 PP).
Fig. 3 Phylogenetic tree from ML analyses based on ITS rDNA sequences for clade D showing a wide sampling of sequences per species. Numbers on the branches indicate bootstrap support ≥70% in MP and ML analyses and posterior probabilities ≥0·95 from the Bayesian analysis. The subclades that are supported in the analyses are numbered.
Thirty-four mtSSU sequences were generated in this study and one sequence was taken from GenBank. The alignment consisted of 477 positions with two single-nucleotide gaps at positions 31 and 223. The analysis was conducted using alignments with and without gaps and produced the same results. Additionally, missing bases were present near the beginning and end of 4 sequences (1 base in 2 sequences, 7 bases in 1 sequence, and 83 in 1 sequence). The mtSSU haplotype network produced 24 haplotypes from 33 sequences. The haplotype network is congruent with the ITS rDNA results, showing that the section Crustaceae is divided into two groups. The species included in the section Tenues are closely related to C. rangiferina and C. stygia (Fig. 4). The haplotypes of C. arbuscula and C. mitis were mixed but they appeared to be monophyletic in the ITS phylogeny, despite a lack of statistical support. One haplotype was shared between C. rangiferina and C. stygia. All other haplotypes are represented by a single species each.
Fig. 4 Haplotype network of mtSSU sequences showing the relationship between Cladonia taxa (see legend). Thallus branching pattern is indicated beside each species with the number representing the branching, and the number in brackets representing the less common branching. The small solid black dots indicate distance between haplotypes and the size of the circles represents number of haplotypes (1 to 4 haplotypes). The numbers with double dashes on the lineages represent number of changes between the dots. MB=Manitoba; NL=Newfoundland and Labrador; SC=South Carolina; GE=Georgia; NJ=New Jersey; AL=Alaska; LI=Lithuania; DN=Denmark; BO=Bolivia; SP=Spain; PT=Portugal.
Three duplet haplotype networks are shown that represent paraphyletic species groups from Fig. 2 and 3 (Fig. 5). The haplotype network of C. portentosa-C. terrae-novae was based on 560 aligned positions and produced 11 haplotypes from 14 sequences and no haplotypes were shared between species (Fig. 5A). The haplotype network of C. arbuscula-C. mitis was based on 564 aligned positions and produced 19 haplotypes from 25 sequences, with two of the haplotypes shared between species and incomplete clustering of the haplotypes within species (Fig. 5B). One loop (dotted line) in the network indicated homoplasy. Twenty-seven sequences of C. rangiferina-C. stygia were based on 472 aligned positions and produced 23 haplotypes with no haplotypes shared between the two species, but the C. stygia haplotypes were intermixed with those of C. rangiferina (Fig. 5C). Five loops (dotted lines) in the network indicated homoplasy.
Fig. 5 Haplotype network of ITS rDNA sequences showing relationship between A) C. portentosa and C. terrae-novae, B) C. arbuscula and C. mitis, and C) C. rangiferina and C. stygia (see legends). The small solid black dots indicate distance between haplotypes and the size of the circles represents number of haplotypes (1 to 3 haplotypes). MB=Manitoba; NL=Newfoundland and Labrador; BC=British Columbia; NB=New Brunswick; NS=Nova Scotia; ON=Ontario; YK=Yukon; SC=South Carolina; GE=Georgia; NJ=New Jersey; WA=Washington; AL=Alaska; LI=Lithuania; DN=Denmark; BO=Bolivia; SP=Spain; PT=Portugal; GM=Germany; AR=Argentina; FN=Finland; GL=Greenland; CH=China; GU=Guyana; ID=India; UnK=unknown location.
The AMOVA analysis showed low to moderate levels of population differentiation between species duplets (Table 3) with low to moderate Phi values. The Phi statistic is a measure of allelic differentiation among subpopulations. A population which represents no allelic differentiation has a Phi value of 0 and a high level of differentiation when the Phi value is 1·0. In this study a population was defined as a species duplet and the subpopulations as the species. A null hypothesis of homogeneity of variance was tested for each of the three species duplets. Significant (P<0·05) genetic differences are shown for C. rangiferina-C. stygia and for C. mitis-C. arbuscula, suggesting a low level of homogeneity among the subpopulations (species), which implies that the species are differentiated from one another. The genetic differentiation for C. portentosa-C. terrae-novae was not significant, suggesting no differentiation between species of the species duplet. The partitioning of the total variance shows a higher level of variance within species than between species for all three duplets. Separate species are supported for each of C. arbuscula, C. mitis, C. rangiferina and C. stygia, but not for C. portentosa and C. terrae-novae.
Table 3 Results of AMOVA analyses between and within species for each species duplet; C. rangiferina-C. stygia, C. arbuscula-C. mitis, and C. portentosa-C. terrae-novae
d.f.=degrees of freedom; SS=sums of squares; MS=mean of squares; variance partitioning includes both observed and percent of total; Phi statistic=fixation index; P-value with significance inferred at 0·05.
Discussion
Phylogenetic relationships of the reindeer lichens
This is the most comprehensive study to date on the phylogeny of the reindeer lichens. It includes 22 species and several specimens for many of them, with an emphasis on C. rangiferina-C. stygia, C. mitis-C. arbuscula and also C. portentosa-C. terrae-novae. A number of previous phylogenetic studies proved that the former genus Cladina was not monophyletic (DePriest et al. Reference DePriest, Piercey-Normore, Sikaroodi, Kärkkäinen and Oksanen1999, Reference DePriest, Piercey-Normore, Sikaroodi, Kärkkäinen, Oksanen, Yahr and Ahti2000; Guo & Kashiwadani Reference Guo and Kashiwadani2004; Stenroos et al. Reference Stenroos, Pino-Bodas, Weckman and Ahti2015); nevertheless, the species sampling was scarce and the composition of the different groups was not complete. Guo & Kashiwadani (Reference Guo and Kashiwadani2004) and Stenroos et al. (Reference Stenroos, Pino-Bodas, Weckman and Ahti2015) concluded that the reindeer species form three independent groups, but relationships among them could not be determined. The results of the present study are consistent with these findings. While some morphological characteristics were not useful in elucidating the evolutionary relationships within certain subdivisions of Cladonia (Stenroos et al. Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002), thallus branching pattern showed a trend in species discrimination in this study. The anisotomic branching was common to both sections Crustaceae and Tenues and the isotomic branching was present in section Impexae.
The section Impexae is monophyletic and basal to the species in this phylogeny (Fig. 1). Tenues, along with most of the species of Crustaceae, forms a monophyletic group that is related to a clade consisting of species from the sections Cocciferae and Perviae, and the groups Amaurocraeae and Divaricatae (according to Stenroos et al. Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002). In the analyses by Stenroos et al. (Reference Stenroos, Pino-Bodas, Weckman and Ahti2015) this group is related to Cladonia uncialis, though this relationship was not supported by most of the analyses here. Guo & Kashiwadani (Reference Guo and Kashiwadani2004) obtained a similar result. The third group is formed by C. arbuscula, C. densissima, C. mitis, and C. submitis. The phylogenetic relationships within this group of the genus Cladonia, however, are not statistically supported in our analyses. The phylogenetic analyses by Stenroos et al. (Reference Stenroos, Pino-Bodas, Weckman and Ahti2015) showed this group to be related to the group Divaricatae, but this relationship was not supported. The study of additional loci will be necessary in order to more accurately understand the phylogenetic relationships among the different groups of reindeer lichens within the genus Cladonia.
Clades A and B: section Impexae
The section Impexae sensu Ahti (Reference Ahti2000) is monophyletic, including C. confusa, C. delavayi, C. evansii, C. mediterranea, C. portentosa, C. stellaris, and C. terrae-novae. The characters used to define this section were the presence of dichotomous or trichotomous branching, hyaline slime in conidiomata (exception: C. stellaris has red slime), and the presence of perlatolic acid (Ahti Reference Ahti1984, Reference Ahti2000). Within section Impexae, C. portentosa (17 sequences) and C. terrae-novae (4 sequences) are morphologically similar to one another. They both have thallus branches that diverge in threes and sometimes twos, and both grow in a boggy habitat. They differ in that C. terrae-novae produces atranorin and perlatolic acid, while C. portentosa produces perlatolic acid alone (both contain usnic acid in addition although it may occasionally be absent; Ahti Reference Ahti1961; Orange Reference Orange1993), and Ahti (Reference Ahti1961) mentions minor differences in branching and surface structure. The two species are closely related but they are geographically separated from one another in North America where C. portentosa is distributed along the west coast (recognized as C. portentosa subsp. pacifica; Ahti Reference Ahti1961, Reference Ahti1984) and C. terrae-novae along the east coast. Allopatry may encourage divergence between these species but if the period of time has not been sufficient for complete divergence, they would show incomplete genetic divergence and an absence of monophyly. Monophyly of C. portentosa was reported by Smith et al. (Reference Smith, Alphandary, Arvidson, Bono, Chipman, Corkery, DiMegli, Hansen, Isch and McAlpine2012) but they did not include C. terrae-novae in their analysis for comparison. While they could not be distinguished by the ITS rDNA phylogeny, the mtSSU haplotype network clearly separated the species, showing C. confusa and C. evansii to fall between them. The AMOVA analysis suggests that in addition to gene flow, another possible explanation is shared ancestral polymorphism where characters may be part of the reaction norm in both species but those characters do not exist in the current niches. The sample size for this species duplet was low, which may have biased the analysis (Fitzpatrick Reference Fitzpatrick2009). These conflicting results suggest that further study of these two species is still needed.
Cladonia stellaris is strongly supported as monophyletic and C. delavayi is basal to C. stellaris, which is consistent with Stenroos et al. (Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002). Cladonia stellaris occupies a basal position in the Impexae (Fig. 1) and is a derived species in the evolution of the Cladoniaceae (Stenroos et al. Reference Stenroos, Ahti and Hyvönen1997). Cladonia stellaris is more closely linked with other members of the Impexae in the mtSSU haplotype network than in the ITS rDNA phylogeny, which is consistent with the morphology. Choisy (Reference Choisy1928) postulated that C. stellaris originated from the ancestors of Cocciferae and Perviae. While C. delavayi was originally thought to be a member of Unciales (des Abbayes Reference Abbayes1958), it was proposed to move it to supergroup Crustaceae and group Cladinae because of this affiliation (Stenroos et al. Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002). Nevertheless, C. delavayi could be a new section of Cladonia on its own since both morphologically and chemically it is very different from the species included in section Impexae. Cladonia delavayi has thickly corticate podetia and, in addition to usnic acid, contains 4-O-methylcryptochlorophaeic and cryptochlorophaeic acids, characters not shared by any other species in Impexae or with C. stellaris. Molecular data would support this separation (Figs 1 & 2A).
Clade D: sections Crustaceae and Tenues
The characters used to separate Crustaceae from Tenues were the branching type and the colour of the slime in conidiomata (Ahti Reference Ahti1984, Reference Ahti2000). The species of section Tenues have dichotomous branching and conidiomata with red slime, while the species included in Crustaceae have predominantly tetrachotomous branching and hyaline slime. However, the phylogenetic results suggest that these characters may have evolved several times independently, they may represent ancestral characters with multiple losses or they may be silenced in some ecological niches.
The position of C. rotundata, C. argentea, and C. dendroides, nested within a clade of C. rangiferina, is consistent with Stenroos et al. (Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002). Cladonia rangiferina and C. stygia were shown to be closely related in this study, which is also supported by Stenroos et al. (Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002). The similar branching pattern may be explained by the non-monophyly of these two species, and the homoplasy in the mtSSU network. The two species differ by moist (C. stygia) and dry (C. rangiferina) habitats (Ruoss & Ahti Reference Ruoss and Ahti1989) but they may overlap in habitats with moderate levels of moisture. Additionally, the stereome in C. stygia is black and the slime in the conidiomata is red, whereas the stereome of C. rangiferina is grey or brown and the slime is colourless (Ahti & Hyvönen Reference Ahti and Hyvönen1985; Ruoss & Ahti Reference Ruoss and Ahti1989). Haplotype networks and AMOVA analysis imply frequent homoplasy and genetic differentiation between C. rangiferina and C. stygia. Therefore, they may not be reproductively isolated from each other and are in the early stages of speciation. The close physical proximity in geographical distribution may encourage interbreeding and obscure genetic divergence between these species.
The two species C. ciliata (six sequences) and C. subtenuis (eight sequences) are both monophyletic in the ITS rDNA phylogeny and are similar in their thallus branching pattern, but they show homoplasy with the Crustaceae in the mtSSU network. Yahr et al. (Reference Yahr, Vilgalys and DePriest2006) reported a low level of population structure and inferred that recombination was occurring within C. subtenuis. Monophyly of C. ciliata was also reported by Smith et al. (Reference Smith, Alphandary, Arvidson, Bono, Chipman, Corkery, DiMegli, Hansen, Isch and McAlpine2012); however, the two colour variants (chemotypes) C. ciliata var. ciliata (=f. ciliata) and C. ciliata var. tenuis (=f. flavicans) were not resolved, which supports the finding in this study. Moreover, the findings of this study support the synonymy of the sections Crustaceae (Cladina) and Tenues as recommended by Stenroos et al. (Reference Stenroos, Hyvonen, Myllys, Thell and Ahti2002).
Clade C: Cladonia arbuscula and related species
Cladonia arbuscula and C. mitis (together with C. submitis and C. densissima) formed a highly supported clade in Fig. 1 but not in Fig. 2B which included a larger number of sequences. Possible explanations include the potential for paralogous ITS rDNA regions (Buckler et al. Reference Buckler, Ippolito and Holtsford1997), a failure of concerted evolution in the nuclear ribosomal repeats (Ambrose & Crease Reference Ambrose and Crease2011), or divergence within C. arbuscula that may be detected using other genes or an increased number of specimens representing all the subspecies. The difficulty in separating the larger number of specimens of C. arbuscula and C. mitis was also consistent with Piercey-Normore et al. (Reference Piercey-Normore, Ahti and Goward2010) and with the mtSSU network (Fig. 4). However, Smith et al. (Reference Smith, Alphandary, Arvidson, Bono, Chipman, Corkery, DiMegli, Hansen, Isch and McAlpine2012) showed both C. mitis and C. arbuscula to be monophyletic species although they had fewer specimens. An analysis with a larger number of specimens would have a higher probability of showing paraphyly than one with fewer specimens. Multiple gene phylogenies showed that C. mitis is supported as a monophyletic species when beta-tubulin, GAPDH, a group 1 intron in nuclear 18S rDNA, and the ITS rDNA were used in the phylogeny, but it was paraphyletic when the intron was omitted from the group of genes or when either gene was used alone in the analysis (Myllys et al. Reference Myllys, Stenroos, Thell and Ahti2003). The reticulate nature of the haplotype analysis in this study, with one case of homoplasy, might suggest a lack of reproductive isolation between these two species. However, the AMOVA analysis indicates only a moderate level of genetic differentiation and a significant P-value, suggesting the species are genetically different from one another. The diagnostic characters overlap between the species, where C. arbuscula has denser branching of the apices with more browned and curved branch tips than C. mitis. Cladonia mitis produces usnic acid alone (or with rangiformic acid), whereas C. arbuscula produces both usnic acid and fumarprotocetraric acid, but this feature can be variable (Ruoss Reference Ruoss1987b ; Ruoss & Ahti Reference Ruoss and Ahti1989). The close evolutionary relationship between the two species and their physical proximity in similar habitats provides opportunities for gene flow and therefore interbreeding. Consequently, they may have had only a short history of reproductive isolation (Myllys et al. Reference Myllys, Stenroos, Thell and Ahti2003), resulting in low resolution of C. arbuscula and C. mitis where sequence divergence (speciation) has lagged behind morphological evolution. Other evolutionary processes may also have influenced the patterns observed in this study, such as incomplete lineage sorting through speciation (Knowles & Carstens Reference Knowles and Carstens2007). A coalescent-based approach using multiple loci will improve the resolution but may not remove the effects of incomplete lineage sorting, depending on the extent of speciation (Knowles & Carstens Reference Knowles and Carstens2007).
Potential bias from ITS rDNA and sampling
The lack of monophyly observed in C. arbuscula and C. mitis, and similarly in C. rangiferina and C. stygia, might suggest there is insufficient phylogenetic signal in the ITS rDNA region to resolve the morphological differences between these species. Among the several gene regions proposed for species discrimination in fungi (Taylor et al. Reference Taylor, Geiser, Burt and Koufopanou1999, Reference Taylor, Jacobson, Kroken, Kasuga, Geiser, Hibbett and Fisher2000; Myllys et al. Reference Myllys, Stenroos, Thell and Ahti2003), Pino-Bodas et al. (Reference Pino-Bodas, Martin, Burgaz and Lumbsch2013) concluded that the best combination for barcoding in Cladonia is RPB2 and ITS rDNA. The ITS rDNA region was previously supported by Schoch et al. (Reference Schoch, Seifert, Huhndorf, Robert, Spouge, Levesque and Chen2012) as a potential barcoding marker for fungi. While the ITS rDNA is widely used, species discrimination using ITS rDNA has been previously shown to be a challenge with some members of Cladonia (Fontaine et al. Reference Fontaine, Ahti and Piercey-Normore2010; Kotelko & Piercey-Normore Reference Kotelko and Piercey-Normore2010; Kelly et al. Reference Kelly, Hollingsworth, Coppins, Ellis, Harrold, Tosh and Yahr2011; Pino-Bodas et al. Reference Pino-Bodas, Burgaz, Martin and Lumbsch2011). The greater intraspecific variation in the ITS rDNA observed with species of Cladonia such as C. arbuscula, C. mitis, C. rangiferina and C. stygia, might also suggest that evolutionary processes such as incomplete lineage sorting in a recent divergence obscures species delimitation (Knowles & Carstens Reference Knowles and Carstens2007). Ribosomal DNA has been subject to interpretations of having divergent paralogs (Buckler et al. Reference Buckler, Ippolito and Holtsford1997) or failure of concerted evolution (Ambrose & Crease Reference Ambrose and Crease2011) in a diversity of organisms. These explanations cannot be ignored in the interpretation of ribosomal DNA patterns.
In conclusion, the current study supported monophyly for five of 22 species in Cladonia sections Crustaceae and Tenues, some of which have also been shown to be monophyletic by other studies, but the 17 other species were not supported and Impexae was divided into two clades using a phylogenetic analysis of the ITS rDNA and a haplotype network of the mtSSU gene. The mtSSU network also illustrated a morphological trend of the thallus branching pattern in these lichens, where members of section Impexae have isotomic branching and Tenues and Crustaceae have anisotomic branching. Incomplete lineage sorting, recombination, gene flow, and recent divergence were considered as explanations for the reticulate nature of the haplotype networks of some species. These results emphasize the importance of examining the non-monophyletic species using multiple loci for a coalescence-based approach. In addition, further investigation of gene flow and recombination between and within the species duplets reported in this study might reveal more about the evolutionary status of these species.
The authors thank T. Booth (University of Manitoba) for providing light microscopic imaging facilities, R. Kotelko for providing field collections, and S. Toni for laboratory assistance. Funding was provided by the Natural Science and Engineering Research Council (NSERC) for a Canada Graduate Scholarship to SA and Discovery Grant to MPN. This project has received funding from the European Union’s Seventh Framework Programme for research, technological development and demonstration under grant agreement No PIEF-GA-2013-625653 (RPB, SS).
