Genome-wide assessment of putative endemism and phylogeography of Cladonia sandstedei (Ascomycota: Cladoniaceae) in the Caribbean

Abstract

Cladonia sandstedei forms cushion-shaped lichens that colonize open environments and is distributed throughout the Caribbean and the south-eastern United States. It co-occurs in parts of its range with C. subtenuis, a morphologically similar taxon that is distinguished from the former by the presence of usnic acid. Preliminary phylogenetic analysis with the RPB2 and TEF-1α loci revealed that these taxa were closely related, but relationships were inconsistent among markers. Here, we combined phylogenetic and population genomic analyses based on RADseq data to clarify the evolutionary relationships and phylogeography of these taxa. Both approaches indicate that the taxa cannot be separated based on secondary metabolites, as previously proposed, but instead form a complex composed of several lineages, largely unrelated to chemistry but with a strong geographical structure in their genetic variation. Continental populations formerly separated under the names C. sandstedei and C. subtenuis were closely related to each other. A similar pattern was observed in the Jamaican counterparts of these taxa, suggesting homoplasy of secondary chemistry. Discriminant Analysis of Principal Components (DAPC) hinted at potential conspecificity between populations in Cuba and Puerto Rico on one hand, and between Jamaica and the continental US on the other; however, phylogenetic analysis and other population-level analyses (PCA and fineRADstructure) suggested that both insular and continental populations were more likely to be reproductively isolated from each other. Based on this, we propose to recognize only one species for the entire complex, under the older name C. sandstedei, with the four spatially structured clades as subspecies: C. sandstedei subsp. sandstedei (restricted to Jamaica) and C. sandstedei subsp. subtenuis comb. nov. (restricted to continental North America) exhibit several chemosyndromes variably containing usnic acid and/or atranorin. The two additional subspecies described here as new, C. sandstedei subsp. cubana and C. sandstedei subsp. landroniana, exhibit the atranorin chemosyndrome and are restricted to Cuba and Puerto Rico, respectively. Our work reaffirms the power of combining RADseq-based phylogenetics and population genetics to disentangle taxonomic and evolutionary histories in poorly understood, closely related and phenotypically similar lichen-forming fungal species.

Introduction

A defining character of our planet's biosphere is the uneven distribution of species diversity (Boenigk et al. 2015). Global biodiversity is shaped by diverse climates and regions rich in endemic species (Cowling 2001; Freeman & Pennell 2021), but limited taxonomic knowledge, especially in biodiverse areas such as the tropics, can obscure these patterns. The insular Caribbean serves as an ideal case of this issue. Well known for its high degree of endemism in animals and vascular plants (Losos 2009; Nieto-Blázquez et al. 2017), the Caribbean is also emerging as a region with unappreciated diversity, particularly in less-studied groups such as lichenized fungi (Mercado-Díaz et al. 2014, 2020; Moncada et al. 2018).

Cladonia sandstedei (Abbayes) Ahti is a cushion-forming ‘reindeer’ lichen that has a primarily Caribbean distribution but also occurs in the south-eastern United States, predominantly in Florida (Ahti 1984, 2000) (Fig. 1A). It is a soil-dwelling taxon that mostly colonizes open environments but can also be found in partly shaded situations. Within Puerto Rico, a peculiar aspect of C. sandstedei is its disjunct distribution, with populations known from high-elevation forests in the Maricao State Forest and sea-level scrub-type ecosystems of the Tortugero Lagoon in the Vega Baja and Manatí municipalities (Ahti 2000). The peculiar ecology of the second population corresponds with its ecology in the south-eastern United States but contrasts with its tendency to occur in high-elevation environments on other islands of the Greater Antilles (i.e. Cuba, Jamaica and Hispaniola). In addition to elevational differences and concomitant climate variation, the two areas in Puerto Rico also exhibit contrasting soil types. Specifically, Maricao features a serpentinite-derived soil (Ricart-Pujals & Padrón-Vélez 2010), whereas the Vega Baja site is characterized by white siliceous sand (Lugo et al. 2001). In the south-east USA and the Bahamas, the species is mostly distributed in low-elevation habitats with sandy soils or other soil types.

Figure 1. A, habit of Cladonia sandstedei from the population in Vega Baja, Puerto Rico. B & C, phylogenetic trees (ML) based on gene markers RPB2 (B) and TEF-1α (C) for a selection of C. sandstedei individuals from main populations from Puerto Rico (Vega Baja and Maricao) and continental samples of C. subtenuis. Thick branches indicate strongly supported clades (BS > 70). In colour online.

Cognizant of the importance of environmental differences as drivers of diversification in lichenized fungi (Kraichak et al. 2015; Huang et al. 2019), it was initially hypothesized that morphologically and chemically indistinguishable populations of Cladonia sandstedei within Puerto Rico represented separate, cryptic lineages. To test this, exploratory phylogenetic analysis using the RPB2 (RNA Polymerase II subunit 2) and TEF-1α (Elongation Factor 1-alpha) barcoding loci was conducted (see Supplementary Material File S1, available online), but the results were inconclusive. Gene histories showed that both populations were either separate genetic entities (RPB2; Fig. 1B) or genetically indistinguishable (TEF-1α; Fig. 1C). More importantly, it was found that C. sandstedei was consistently nested within C. subtenuis (Abbayes) Mattick, another taxon with a similar distribution in the Caribbean and south-eastern US, but further extending its range to Central and Eastern North America (Fig. 2).

Figure 2. Sampling localities (red circles) used in this study for populations of Cladonia sandstedei and C. subtenuis in the Caribbean and continental United States. Triangles show the distribution of C. subtenuis (black) and the distribution of C. sandstedei (orange) according to Global Biodiversity Inventory Facility (GBIF) occurrence records. Cladonia subtenius has also been documented in several localities in Central America, but these records are not included in the GBIF database. In colour online.

Although morphological characters such as cortex texture and branching type have been suggested to differ between C. sandstedei and C. subtenuis (Ahti 2000), in practice these features are challenging to use due to their wide variation. Both taxa contain fumarprotocetraric complex substances as major products, limiting their discriminatory power. Alternatively, the presence of atranorin in C. sandstedei versus usnic acid in C. subtenuis, which results in a white-grey colour versus a light yellow-green colour (Ahti 2000), has been used as the primary method to distinguish between them. Since concentrations of these substances can be low and colour can be influenced by other factors, the separation of these taxa often relies on using potassium hydroxide (KOH) spot tests: a specimen exhibiting a K+ yellow reaction is typically identified as C. sandstedei, whereas a K− reaction identifies individuals as C. subtenuis. However, secondary chemistry in this complex is not as straightforward, and spot tests often fail to produce definitive results (J. A. Mercado-Díaz, personal observation). Furthermore, chemical signatures associated with these species are not consistent in taxonomic treatments, since individuals lacking atranorin have been identified as C. sandstedei and samples with both usnic acid and atranorin (or lacking either substance) have been identified as C. subtenuis (Ahti 2000). The form lacking both atranorin and usnic acid and containing fumarprotocetraric acid has been described as C. subtenuis f. cinerea and is found mostly in the continental United States (Ahti 1984, 2000). Apart from their secondary metabolites, these taxa are morphologically very similar, with a presumed tendency of the branching pattern in C. sandstedei to be slightly denser and more isotomic than in C. subtenuis (Ahti 2000; Rosentreter et al. 2015).

Several genetic markers have proved useful for resolving species-level relationships in some groups of Cladonia (Kanz et al. 2015). However, others suffer from insufficient resolution to confidently place boundaries between species (Stenroos et al. 2002; Pino-Bodas et al. 2011, 2013). Considering these limitations and issues inherent in using phenotypic characters that might lack sufficient discriminatory power, resolving species boundaries between C. sandstedei and C. subtenuis, and understanding the phylogeographical structure of populations throughout the Caribbean remains challenging. One alternative is to reassess phylogenetic relationships using a wider sampling of genomic regions. This approach was implemented in a multilocus study by Stenroos et al. (2019), which found support for 11 major clades within Cladonia. Nonetheless, species-level relationships were poorly resolved in several clades, suggesting that more extensive genome sampling might be needed to further clarify relationships at shallower taxonomic levels.

Next-generation sequencing (NGS) approaches, particularly those categorized as reduced-representation sequencing or ‘genotype-by-sequencing’ techniques, are increasingly becoming the option of choice for generating comprehensive genome-wide datasets, particularly in non-model organisms. Among these, Restriction-Site Associated DNA sequencing (RADseq) is prominent as it allows the discovery and genotyping of high-throughput single nucleotide polymorphisms (SNP) at a reasonable cost and without requiring prior genomic information of the taxa under study (Narum et al. 2013; Andrews et al. 2016). RADseq approaches have proved useful for disentangling the genetic basis of poorly understood ecological and evolutionary phenomena in diverse taxa (see table 1 in Narum et al. (2013)). Notable examples include the identification of markers linked to key ecological traits in three-spine sticklebacks (Baird et al. 2008), studies assessing patterns of introgression and hybridization between native and invasive trout species in the western USA (Hohenlohe et al. 2011, 2013) and groundbreaking work resolving species-level relationships in recently diverged Lake Victoria cichlid species (Wagner et al. 2013). By relaxing the prerequisite of prior identification of loci, RADseq-based methods have allowed researchers to perform diverse analyses that were previously unattainable for non-model organisms. This has led to an increase in studies of many groups, including investigations in lichenized fungi, such as the genera Cladonia, Pseudocyphellaria and Usnea (Grewe et al. 2017, 2018; Alonso-García et al. 2021; Widhelm et al. 2021, 2023; Barcenas-Peña et al. 2023; Otero et al. 2023).

Table 1. Population-level statistics of genetic differentiation for a priori groups. Statistics include the number of representative sequences for each group (Num. seq.), number of segregating sites (S), number of haplotypes (h), haplotype diversity (Hd), average number of nucleotide differences (K), nucleotide diversity (π), and nucleotide diversity with Jukes-Cantor correction (πJC). PR = Puerto Rico.

Here, we implemented a RADseq approach to evaluate phylogenetic relationships and explore the genetic structure in the Cladonia sandstedei/subtenuis complex from populations in the Caribbean and in the south-eastern United States. We assessed the degree of divergence between populations within the Caribbean and describe evolutionary relationships between these lineages and those from the continent. Ordination techniques were used to explore population structure and evaluate genetic admixture. We anticipated considerable genetic differentiation associated with geography, including the disjunct populations within Puerto Rico emerging as distinct genetic clusters, separated from other Caribbean and continental populations. Due to lower barriers to dispersal, the degree of admixture and genetic variability was expected to be higher among continental lineages. Phylogenetic patterns were used to evaluate the usefulness of secondary metabolites for species discrimination. However, no specific patterns were necessarily anticipated since secondary chemistry in Cladonia can be variable (Stenroos et al. 2002), sometimes concordant with phylogeny-based species delimitations (Pino-Bodas et al. 2012a) but occasionally overestimating diversity (Pino-Bodas et al. 2012b).

Materials and Methods

Sampling and taxonomic work

We sampled localities where individuals of Cladonia sandstedei, C. subtenuis f. subtenuis and C. subtenuis f. cinerea Ahti species have been previously reported. These species produce irregular to subglobose, cushion-shaped colonies, which we treat as separate individuals (Fig. 1A). A total of 187 samples were collected in the continental United States and the Caribbean (Fig. 2). Most specimens of C. sandstedei and the two ‘formae’ of C. subtenuis were collected in the USA, specifically in the states of Alabama (37), Florida (34), Georgia (65) and Tennessee (4). In Puerto Rico, we collected a total of 31 samples identified as C. sandstedei: 15 from populations in Maricao and 16 from Vega Baja. Seven samples from Cuba were identified as C. sandstedei. Four of these samples were collected as part of a separate effort and were in the same herbarium sheet. These were processed separately since each of these corresponded to an individual ‘cushion’. Among four samples from Jamaica, three were previously identified as C. subtenuis and one as C. sandstedei. Our sampling also included five C. rangiferina (L.) Weber specimens from Canada (4) and the Dominican Republic (1), which were used as outgroups for phylogenetic analysis. A minimum distance of 5 m between individuals was maintained, except for populations in Maricao which occur in relatively restricted patches (c. 100 m²) and for the four Cuban individuals from the same sheet. Herbarium vouchers are deposited in BSC, HAJB, UPR and B.

Preliminary taxonomic identification of samples relied on thallus colour and gross morphology assessment. Spot tests (K) were used to predetermine the presence of atranorin; to confirm the secondary chemistry, high-performance thin-layer chromatography (HPTLC) was also carried out. Besides detecting atranorin, this method can also be used to detect usnic acid and other substances present in these taxa. Solvent C was used for all HPTLC analyses, which followed Lumbsch (2002).

DNA extraction

Small fragments of each sample were soaked overnight in acetone to remove secondary substances. The acetone was removed, and all samples were air-dried before they were manually ground with a mortar and pestle. DNA was then extracted using the ZR Fungal/Bacterial DNA MiniPrep Kit (Zymo Research, Irvine, CA, USA). Other than incubating specimens overnight, the extraction process followed the manufacturer's instructions. A Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) was used to corroborate the DNA concentration of samples, while the DNA quality was assessed via gel electrophoresis. These DNA extractions were used for both single-locus PCR (Supplementary Material File S1, available online) and RAD sequencing. Of the 187 specimens collected for this work, only three samples from Cuba were not processed for phylogenomic analysis due to poor DNA quality.

RADseq library preparation and sequencing

DNA was submitted to the University of Wisconsin-Madison Biotechnology Center for library preparation and sequencing. DNA concentration was verified using the Quant-iT™ PicoGreen® dsDNA Kit (Life Technologies, Grand Island, NY, USA). Libraries were prepared as in Elshire et al. (2011) with minimal modification; in short, 50 ng of DNA was digested using the 5-bp cutter ApeKI (New England Biolabs, Ipswich, MA, USA), after which barcoded adapters amenable to Illumina sequencing were added by ligation with T4 ligase (New England Biolabs, Ipswich, MA). Adapter ligation proceeded in batches of 96 samples which were then pooled and amplified to provide library quantities suitable for sequencing. Adapter dimers were removed by SPRI bead purification. The quality and quantity of the finished libraries were assessed using the Agilent Bioanalyzer High Sensitivity Chip (Agilent Technologies, Inc., Santa Clara, CA, USA) and the Qubit® dsDNA HS Assay Kit (Life Technologies, Grand Island, NY, USA), respectively. Single-end sequencing for an initial set of 15 samples was carried out on a HiSeq 2500 sequencer (Illumina Inc.). Paired-end sequencing for two additional sets, one including 164 samples and the other including the five outgroup samples, were subsequently sequenced with a NovaSeq 6000 System (Illumina Inc.). For these efforts, sequencing targeted 1.3 M reads per sample. Images were analyzed using the standard Illumina Pipeline, v. 1.8.2.

RADseq dataset assembly

The process_radtags pipeline from Stacks v. 2.3 (Rochette et al. 2019) was used to demultiplex raw paired-end reads obtained from the sequencers. Quality control was also carried out with process_radtags. This entailed cleaning (removing uncalled bases), filtering out reads with low-quality scores (Phred score = 33), and removing adapters by trimming reads to a 55bp length. Quality control with process_radtags was also carried out for the set of 15 single-end libraries which were previously demultiplexed. Parameters and command-line settings used are provided in Supplementary Material File S2 (available online). FastQC reports (Banraham Bioinformatics, Cambridge, UK) were utilized for quality control.

We combined the 15 forward reads generated from our single-end sequencing with forward reads of the paired-end sequence datasets to increase the number of individuals used for downstream analysis. A RADseq data assembly was carried out with ipyRAD (Eaton & Overcast 2020) using a reference-based approach, which allowed the filtering of mycobiont loci. Mapping and subsequent filtering of reads was performed using the recently assembled Cladonia rangiferina genome available at NCBI (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_023065455.1/). To facilitate interpretation, raw Illumina RAD sequences are referred to as ‘reads’, and ‘loci’ refers to clustered reads per individual sample. The final matrices are alignments of homologous loci from multiple samples with nucleotide substitutions referred to as ‘SNP’. Sequences were processed in ipyRAD by setting the datatype to ‘gbs’, ploidy to haploid (‘1’), clustering of reads within and between individuals to a similarity threshold of ‘0.90’ and a minimum coverage of four. Seven samples were removed from datasets because they either had less than 1000 recovered loci or belonged to non-target taxa. Subsequent phylogenetic and population genomic analysis was carried out using filtered ipyRAD output files, such as ‘.usnps,’ ‘.alleles’ and ‘.vcf’.

Phylogenetic analysis

Unlinked SNP files (i.e. matrices limited to one SNP per RAD locus) from the filtered RADseq dataset were used for phylogenetic reconstructions with RAxML v. 8.2.11 (Stamatakis 2014). This analysis entailed searching for the best-scoring ML tree under the ASC_GTRGAMMA model with ascertainment bias correction (--asc-corr = lewis). The bootstrap convergence test using the extended majority-rule consensus tree criterion (autoMRE) was used for a posteriori bootstrapping analysis. Phylogenetic trees were first inspected in FigTree v. 1.4.3 (Rambaut 2012) and then plotted with the R package ggtree (Yu 2020).

Population genomics

For population genetic analysis, the dataset was reduced to sites with a minor allele frequency (MAF) ≥ 0.05 and a minimum coverage of 80%. This was accomplished using the program VCFtools v. 0.1.15 (Danecek et al. 2011), which was also used to remove outgroup taxa. This dataset was converted to a genlight object using the R package vcfR (Knaus & Grünwald 2017). The package adegenet v. 2.0.2 (Jombart et al. 2010; Jombart & Ahmed 2011) was then used to convert this file to a genind or genlight object. Additional information settings for haploid genomes and sample population membership were subsequently appended to the dataset.

Population structure was initially explored using Principal Component Analysis (PCA). This analysis was based on the genlight object and was run with the ‘glPca’ function from adegenet. The first three principal components were retained. PCA scores were visualized with the R package ggplot2. Colours and ellipses (level = 95%) denote major clades inferred in our RAxML tree, which correspond to the island-specific clades of Jamaica and Cuba, the main clades corresponding to populations within Puerto Rico (i.e. Maricao and Vega Baja), and the single clade recovered for continental samples (see Results). To further characterize genetic variation within these clades, basic population-level statistics were computed using the ‘Gene Flow and Genetic differentiation’ option in DnaSP (Librado & Rozas 2009).

We also used Discriminant Analysis of Principal Components (DAPC) as implemented in adegenet to explore the grouping of samples that might be linked to population structure. This non-parametric method attempts to summarize genetic variation between groups while overlooking within-group variation (Jombart et al. 2010). In essence, the method performs a Principal Components Analysis (PCA) transformation of the data, which produces a set of uncorrelated variables (principal components) suitable for Discriminant Analysis (DA). Resembling Bayesian clustering methods, individuals are probabilistically assigned to groups (clusters).

Both a priori grouping and de novo clustering of samples was performed with DAPC. As with PCA, clades for each Caribbean island and the major clade for the continental US were used as a priori groups. Since we were interested in evaluating the potential influence of contrasting habitat preferences on Puerto Rico's genetic structuring, the Maricao and Vega Baja populations were kept separate. For de novo clustering, the function ‘find.clusters’ was used to assign samples to groups. Group assignment was based on the ‘diffNgroup’ criterion using the lowest-score Bayesian information criterion (BIC) value for selecting the ‘best’ number of populations (K). For both a priori and de novo clustering analyses, an a-score optimization approach was used to select an ideal number of principal components (PC) to retain for DAPC analysis. This entailed running an initial DAPC using 30 PCs and then using the function ‘optim.a.score’ to determine the final number of PCs to retain. Differentiation between populations was achieved by visualizing discriminant functions and principal components of the DAPC using the function ‘scatter’. Group membership bar plots from DAPC were generated in ggplot2.

Lastly, we used fineRADstructure (Malinsky et al. 2018) to better understand recent shared ancestry between individuals. This program uses haplotype linkage to summarize nearest-neighbour haplotype relationships between individuals and infer populations. A co-ancestry matrix is used to visualize relationships. As a first step, the ‘.alleles’ file from ipyRAD was converted into a fineRADstructure format using the ‘finerad_input.py’ script from fineRADstructure-tools (https://github.com/edgardomortiz/fineRADstructure-tools). During this process, the dataset was reduced to contain only unlinked loci (default parameter) and a minimum sample number of four (--minsample 4). Following the authors’ recommendations, the ‘sampleLD.R’ script from fineRADstructure was then used to reorder loci. A co-ancestry matrix for a haploid dataset (-p 1) was generated with RADpainter, and individuals were assigned to populations using fineSTRUCTURE Markov chain Monte Carlo (MCMC) clustering algorithm with the following arguments: -x 100 000, -z 100 000, and -y 1000. This clustering algorithm also generated a simple coalescent tree to explore the relationships between inferred populations with the following arguments: -m T and -x 10 000. Visualization and generation of the co-ancestry figure were carried out in the program ‘Finestructure GUI’ (Lawson et al. 2012) after uploading the co-ancestry matrix, the inferred MCMC clusters (populations), and the coalescent tree.

Molecular-based taxonomic diagnoses

ITS1 and ITS2 regions for representative specimens of the two new subspecies introduced here were examined to detect the presence of nucleotide variants that may serve as unequivocal identifiers for these subspecies (see Supplementary Material File S1). Short, 11 nucleotide-long barcodes with the variant nucleotide at the centre (position 6) are provided in taxonomic descriptions below to allow unambiguous identification of these subspecies. PCR and cycle sequencing parameters of ITS follow Mercado-Díaz et al. (2023).

Results

RADseq data processing

Of the 187 samples collected, 177 samples were included in the final genomic datasets, five corresponding to Cladonia rangiferina outgroup samples. An average of 2 076 171 (SD = 387 788; range 747 109–2 997 624) raw reads per sample were recovered (Supplementary Material File S3, available online). An average of 18% of reads (SD = 6%) successfully mapped to the reference genome. A statistically significant association between within-sample clusters and the number of mapped reads was found (R² = 0.53, P < 0.001). Conversely, the number of within-sample clusters was strongly correlated with the final number of loci used for analysis (R² = 0.70, P < 0.001) (Supplementary Material File S4, available online). Retained samples yielded an average of 12 744 loci (SD = 3487) which were used for downstream analysis. Samples excluded from genomic datasets were still used for chemical characterization.

Phylogenetic analysis

Our phylogenetic analysis of 172 C. sandstedei/C. subtenuis specimens + five C. rangiferina outgroup samples split C. sandstedei/C. subtenuis individuals from the Caribbean and the south-eastern US into two monophyletic clades with bootstrap support of 100 and 66, respectively (Fig. 3). The Caribbean clade showed strong stratification by island, with strongly supported subclades representing each island. Jamaican samples diverged the earliest within this clade and were sister to all remaining Caribbean samples. A monophyletic Cuban clade was found to be sister to all samples from Puerto Rico. Within this Puerto Rican clade, Vega Baja individuals formed a strongly supported monophyletic clade, while Maricao individuals were paraphyletic. All Caribbean samples contained atranorin, including the three C. subtenuis individuals from Jamaica, which also contained usnic acid.

Figure 3. RAD-seq based phylogenetic analysis (RAxML) of Cladonia sandstedei and C. subtenuis samples from the Caribbean and south-eastern United States. Presence/absence of main substances are shown to the left of tip labels. Coloured maps of geographical areas sampled are used to indicate geographical origin of specimens: United States (grey), Vega Baja, Puerto Rico (yellow), Maricao, Puerto Rico (green), Cuba (magenta), Jamaica (navy blue). Black stars on geographical areas indicate approximate sampling localities. Except for branches with 100% statistical support, which are thickened, bootstrap values for other strongly supported clades (BS 70–99) are placed above branches. The bootstrap value for the continental clade is also placed above its corresponding branch but its marginal support is indicated with an asterisk. Note that while individuals representative of both infraspecific ‘formae’ of C. subtenuis are present in our sampling, we group them all under the name ‘Cladonia subtenuis’.

As indicated, the monophyly of the continental clade was supported by a bootstrap value of 66, and it included samples identified both as Cladonia sandstedei and C. subtenuis (Fig. 3). Both species in this clade were intermixed, and most relationships within the clade received low bootstrap support. Only two subclades, one with five C. subtenuis samples from Georgia and another comprising five C. subtenuis individuals from Alabama, were strongly supported. Outside these and several pairs of genetically similar individuals, relationships within this clade were unresolved. Contrary to the samples in the Caribbean clade, which either contained atranorin or atranorin and usnic acid, a broad diversity of chemosyndromes was evident in the continental clade. This diversity was exemplified by samples producing either atranorin or usnic acid and other samples producing both substances, or lacking them altogether. Samples without usnic acid or atranorin were identified as C. subtenuis f. cinerea. This finding reflected a secondary chemistry pattern in these species with no clear separation between atranorin- versus usnic acid-producing samples.

Population genetic structure

A total of 6697 SNPs were included in the reduced dataset used for population genomic analyses. Most variation in our PCA is associated with PC axis 1 (17.3%) (Fig. 4), which showed a geographical gradient in genetic variation. Continental individuals exhibit high genetic variability, with samples spanning a broad range of values over both PC axes, particularly along PC axis 2 (3.2%). Jamaican and Cuban samples show low within-island variation and cluster between individuals from the continent and Puerto Rico. Populations from Puerto Rico show moderate genetic variation and a strong overlap in PCA space. Additionally, the continental clade yielded the highest values for most population-level statistics estimated with DnaSP (i.e. number of segregating sites (S), number of haplotypes (h), average number of nucleotide differences (K), nucleotide diversity (π), and nucleotide diversity with Jukes-Cantor correction (πJC)) (Table 1). Only haplotype diversity (Hd) was the same for both continental US and Jamaican populations (Hd = 1).

Figure 4. Principal component analysis based on a reduced 6697 SNPs dataset for Caribbean and continental US samples of Cladonia sandstedei and C. subtenuis. Most variation is associated with PC axis 1 (17.3%) which showed a geographical gradient in genetic variation (Vega Baja and Maricao are sampled areas in Puerto Rico). Continental individuals exhibit high genetic variability, with samples spanning a broad range of values over both axes, particularly PC axis 2 (3.2%).

De novo clustering at K = 3 yielded the lowest BIC score for selecting the ‘best’ number of populations (Supplementary Material File S5: Fig. S5.1, available online). Six PCs and two (all) discriminant functions were retained for this analysis. In the resulting DAPC, genomic variation associated with the continent was split into two weakly separated clusters: cluster 1, which merged several continental samples and individuals from Jamaica, and cluster 2, which grouped the rest of the continental samples (Fig. 5). Several continental samples represented admixed individuals. Conversely, samples from Cuba and the two populations from Puerto Rico (Maricao and Vega Baja) were grouped in a third cluster (cluster 3). These clusters show moderate scatter along discriminant functions.

Figure 5. Results from de novo clustering of Cladonia sandstedei and C. subtenuis samples with DAPC at K = 3. Upper part shows the scatterplot for discriminant functions whereas the lower part shows a bar plot with assigned membership probabilities. Each dot and bar represent an individual. PR = Puerto Rico.

We explored de novo clustering at K = 2 and K = 4 since these also yielded comparatively low BIC scores. One PC and one discriminant function were retained for the K = 2 analysis. The resulting DAPC merged Jamaican and continental samples into a single cluster and samples from Cuba, Maricao and Vega Baja into a second cluster (Supplementary Material File S5: Fig. S5.2). At K = 4, six PCs and three discriminant functions were retained. The resulting DAPC showed a clustering of continental samples similar to K = 3, with samples either assigned to one of two clusters (clusters 1 and 3) or representing admixed individuals (Supplementary Material File S5: Fig. S5.3). Jamaican individuals were unambiguously assigned to the subset of continental samples grouped under cluster 3. Samples from Cuba and Puerto Rico formed two weakly separated clusters, one containing six samples from Maricao and all individuals from Cuba (cluster 4) and the other composed of all samples from Vega Baja (cluster 3) and six samples from Maricao (cluster 2). One sample from Maricao exhibited mixed Maricao/Vega Baja ancestry.

Five PCs together with four (all) discriminant functions were retained for DAPC analysis based on a priori groupings. This analysis shows a clear separation between continental samples and Caribbean samples (Supplementary Material File S5: Fig. S5.4). Within the Caribbean, samples from Jamaica and Cuba are unambiguously assigned to separate groups. Although several Maricao samples were assigned with high probability to Vega Baja, group membership for most samples was concordant with geographical location. However, individuals from these two populations failed to separate along discriminant function space. Continental samples formed a single group with no observed admixture.

Finally, fineRADstructure analysis inferred a total of 38 populations (Fig. 6). Levels of shared co-ancestry mirrored patterns of clade support in our RAxML tree. Thus, shared ancestry was higher in samples within each island and in continental samples that formed strongly supported clades. It is worth noting that co-ancestry between Cuban samples (‘Pop 28’) and Puerto Rican samples (‘Pop 3’, ‘Pop 5’, ‘Pop 17’, ‘Pop 21’, ‘Pop 23’, ‘Pop 26’, ‘Pop 35’) was higher than co-ancestry between Jamaican samples (‘Pop 38’) and Puerto Rican samples. Moreover, levels of shared ancestry between Jamaican and Cuban samples resembled the degree of co-ancestry observed between Cuban and Puerto Rican samples.

Figure 6. fineRADstructure co-ancestry matrix showing recent shared ancestry between Cladonia sandstedei and C. subtenuis individuals from the Caribbean and south-eastern US. Populations inferred are numbered from 1–38; these broadly correspond to strongly supported clades in the phylogenetic analysis. Indication of the geographical origin of samples follows the same scheme used in Fig. 3 and Fig. 4. The scale to the right depicts levels of shared co-ancestry, ranging from high (blue) to low (yellow). In colour online.

Discussion

Phylogenetic patterns

We found a remarkably strong geographical pattern in our phylogenetic analysis. This is especially evident in the Caribbean, where strongly supported island-specific clades were inferred for C. sandstedei from Cuba and Puerto Rico, and C. sandstedei/C. subtenuis from Jamaica. The strong divergence associated with these clades resembles patterns observed in other metazoan lineages and plant groups from this region, many of which are considered putative endemics (Judd 2001; Michelangeli et al. 2008; Rodríguez et al. 2010; Alonso et al. 2012; Reynolds et al. 2013; Matos-Maraví et al. 2014; Nieto-Blázquez et al. 2020). Phylogenetic studies of lichen-forming fungal species with predominantly Caribbean distributions are scant, but evidence to date suggests that lineages likely restricted to the Caribbean exhibit similar degrees of phylogenetic uniqueness (Lücking et al. 2020; Mercado-Díaz et al. 2020, 2023). It should be noted, however, that inferring the degree of divergence of a clade recovered from phylogenetic analysis based on SNP data should be interpreted with caution. A simulation study demonstrated that, without proper correction, SNP-based phylogenetic analysis might introduce systematic errors in phylogenetic inferences, including biases in branch lengths (Bertels et al. 2014). Overestimated branch lengths might result even when ascertainment bias corrections, similar to the one performed in this work (i.e. asc_corr=lewis), are carried out (Leaché et al. 2015). This type of issue might also affect topological patterns, which can negatively impact downstream analyses (Lewis 2001; Leaché et al. 2015).

The contrasting patterns of phylogenetic support in continental versus island clades and the lack of clear phylogenetic signals in secondary chemistry deserve further consideration. Despite representing different chemotypes, C. sandstedei and C. subtenuis individuals from the continent were not recovered as separate clades. Similarly, Jamaican samples from both currently accepted species were more closely related to each other than they were to any other sample, suggesting that separation based on secondary metabolites does not reflect phylogenetic relationships between these species. This homoplasy in chemical signature might lead to patterns of overestimated diversity, as evidenced in Pino-Bodas et al. (2012b). In contrast, C. sandstedei from Cuba and Puerto Rico exhibit homogeneous secondary chemistry, yet they were placed in distinct island-specific clades. These patterns suggest that chemical signature in this group may also result in underestimated diversity (see Taxonomy section below). Chemical races with and without usnic acid are also known in species closely related to C. sandstedei/C. subtenuis (e.g. C. ciliata), but they have been used at infraspecific ranks (i.e. C. ciliata var. ciliata Stirt. vs C. ciliata var. tenuis (Flörke) Ahti). These observations further demonstrate that, at least at the species level, secondary metabolites are of little taxonomic value in this complex (Ahti 2000).

The ITS-based study by Yahr et al. (2006) of C. subtenuis in the US revealed patterns of statistically supported clade substructure that are not straightforwardly reflected in our continental clade. We suspect that limited geographical and taxonomic sampling overlap between that study and ours obscure these apparent differences. In the first place, the taxonomic sampling of Yahr et al. (2006) excludes C. sandstedei individuals, which limits the possibility of adequately assessing topological congruence between phylogenies. Furthermore, areas close to only two of the 11 sampling areas that Yahr et al. (2006) visited during their work (i.e. FL3 and GA, in the states of Florida and Georgia) were visited for this study. It is, therefore, possible that haplotypes that showed those patterns of clade substructure were missed during our efforts. Nonetheless, the only continental samples that showed some degree of clade substructure in our study belonged to C. subtenuis individuals from our Georgia and Alabama sites (the latter being geographically close to our Florida sites). This supports the hypothesis of potential genetic correspondence between these samples and the samples showing strong clade substructure in Yahr et al. (2006).

Population genomics

Analysis of population genomic structure provides additional details about genetic variation. As suggested by the PCA and further validated in the phylogenetic analysis, genetic differences along an island-to-continent gradient represent the strongest axis of variation in our data. Continental samples showed the largest scatter of PC scores and yielded the highest haplotype and nucleotide diversity values, which suggests that genetic diversity and variation are highest in continental versus insular lineages. Such observations broadly follow expectations from population genetics theory, which predicts higher base pair differences in gene copies from larger (e.g. continent) versus smaller (e.g. island) populations (Kimura 1983). However, conclusions are hard to draw due to the limited number of island samples included. Despite this shortcoming, varied genetic diversity at local and regional scales has been previously documented for continental populations of C. subtenuis. For instance, Beard & DePriest (1996) detected restriction site polymorphisms associated with the small subunit (SSU) ribosomal DNA marker in individuals from the same locality and between individuals from different populations in the United States. Between-population differences in the accumulation of nucleotide changes in one group I intron were also observed (Beard & DePriest 1996). Yahr et al. (2006) found high diversity in fungal ITS genotypes (32) in a set of 79 samples collected throughout the Eastern United States. Further characterization of genetic diversity in island lineages is still required, which highlights the need to increase sampling efforts, particularly of individuals with the C. sandstedei phenotype, for which genetic data was unavailable before the present work.

Clustering analysis with DAPC further helped visualize potential boundaries of population subdivision in our data. Both a priori and de novo clustering agreed with the separation of most insular populations from continental individuals. A priori clustering mirrored patterns observed in our phylogeny and PCA analysis. Except for several samples from Maricao, this method assigned 100% membership probabilities to individuals in their respective a priori groups.

Groupings recovered when exploring optimal K values for de novo clustering provided additional insight into potentially hidden patterns of genetic structuring. Conclusively, the most noteworthy feature in this respect was a clear separation of Cuban + Puerto Rican populations from individuals in Jamaica and the continental US. This feature is also consistent in the DAPC plots generated using slightly less optimal K values (Supplementary Material File S5, available online), which suggests a true genetic divide between these populations. Furthermore, the genetic variation found at the continental level is split into two clusters in plots at K = 3 and K = 4. Several evolutionary processes might be invoked to explain such patterns, including sympatric speciation processes at their early stages or genetic structuring associated with isolation by distance between populations. To address the latter, we performed a post -hoc analysis evaluating whether genetic dissimilarity correlated with geographical distance between individuals (Supplementary Material File S6, available online). Although the correlation was weak, we found the association between these variables to be statistically significant, suggesting that distance between populations might partly explain the observed splitting. Genetic partitioning at the continental level is somewhat surprising because landscape connectivity should promote substantial gene flow in groups such as Cladonia which have putatively high dispersal capabilities (Myllys et al. 2003; Alonso-García et al. 2021). This view is supported by weak genetic structure (Fst) found in an ITS-based population genetic analysis on widely dispersed populations of C. subtenuis in the continental US (Yahr et al. 2006) and the lack of significant genetic differentiation between populations of C. stellaris (Opiz) Pouzar & Vězda observed along a continental-level latitudinal gradient in Canada (Alonso-García et al. 2021). Additional population genetics work will be needed to dissect further the factors promoting this spatial pattern in genetic structuring.

Lastly, candidate populations inferred by fineRADstructure corresponded to strongly supported clades in our phylogenetic analysis, reaffirming that genetic variation in our data is geographically structured. Since the method is optimized for detecting recent coalescence, the analysis suggests that observed phylogeographical patterns might be linked to recent shared ancestry between individuals from populations within each island and the continent. Levels of co-ancestry between island versus continental populations were relatively low, further supporting the separation of these populations. In agreement with phylogenetic analysis, co-ancestry between Jamaican samples and continental US populations was also low, suggesting that the onset of genetic differentiation between these populations is somewhat distant. As with the Jamaican population, the degree of shared co-ancestry among Cuban samples was similar to levels observed in a study that concluded that Usnea aurantiacoatra and U. antarctica populations were different species (Grewe et al. 2018). Yet, as suggested in the DAPC results, it was clear that shared ancestry between Cuban and Puerto Rican populations was higher than that between Puerto Rican and Jamaican populations. On the other hand, lower co-ancestry associated with the rest of the samples resembled patterns observed by Alonso-García et al. (2021) for clade-defined populations of Cladonia stellaris in north-eastern North America.

Considering their habitat preferences, it was somewhat surprising to find that island-level populations of Cladonia sandstedei/C. subtenuis exhibited some degree of genetic isolation from the continent. While occasionally found under moderately closed canopies, these species are often found at low- to mid-elevation open habitats which are typically associated with higher levels of gene flow and, thus, wider distributions (e.g. Sticta scabrosa; Moncada et al. 2021). Underlying this pattern are presumed high dispersal capacities in these species, which result from small propagules (e.g. spores, minute thallus fragments) assumed to be carried away easily by wind currents and/or animal vectors. In our case, we suspect that limitations associated with dispersal modes and ecological factors collude to restrict gene flow among populations. For C. subtenuis, patterns of association between mycobionts and photobionts have been shown to align better with apparent landscape-level spore dispersal (Yahr et al. 2006). If spore movement underlies gene flow, we suspect that genetic differentiation between these populations could indirectly result from unsuccessful establishment, which may happen either via ineffective spore dispersal (e.g. spores being damaged by physical factors) and/or from failure to re-lichenize with a suitable photobiont partner. In fact, this latter point is in line with findings from Yahr et al. (2006), who posed that re-lichenization in C. subtenius and its associated algal partners is determined by the composition of compatible algal partners at a site and local environmental conditions, suggesting a high degree of ecological specialization. Conversely, genetic differentiation between these populations may also emerge if gene flow is most strongly modulated by vegetative dispersal. As recently demonstrated for C. stellaris (Alonso-García et al. 2022), thallus fragmentation also constitutes the main vegetative dispersal mechanism for this group (Yahr et al. 2006). Thallus fragments are, in most cases, larger (and heavier) than other more common vegetative propagules (e.g. soredia, isidia), suggesting that stronger biotic and/or abiotic forces would be required for successful dispersal to other areas. In this context, genetic isolation of insular populations from the mainland due to vegetative dispersal constraints would seem plausible especially because, unlike other studies (e.g. Alonso-García et al. 2022), our populations are separated by stretches of sharply contrasting habitats (i.e. ocean vs land). These views are supported by work on C. subcervicornis (Vain.) Kernst. occurring in open environments along several islands off the west coast of Norway. This study found that populations on one island were strongly genetically differentiated from the rest (Printzen & Ekman 2003). Considering the proximity between these islands (< 50 km), their data suggest that Cladonia species from open environments might also experience effective barriers to dispersal, even over short distances. Furthermore, being restricted to sites with relatively rare soil types might also help explain moderate genetic divergence associated with the Maricao and Vega Baja populations within Puerto Rico.

Caution should be exerted when interpreting our findings, as methods are prone to biases that derive from different issues. For instance, our work is constrained by relatively limited sampling in areas where these species have been recorded. In the Greater Antilles, this is particularly notable for the Dominican Republic, from where we have no representative specimens, and for Cuba and Jamaica, where the sampling was limited. As such, the genetic variability in these populations remains somewhat obscure. We also lack samples from populations in central Florida, where the type of C. subtenuis was collected, and from the southern tip of this state, where both species have been documented (see Fig. 2). Also, potential methodological artifacts might underlie the grouping of samples with DAPC. Previous work has shown that clustering methods and similar tools for the visualization of population structure could be sensitive to uneven sampling (Shringarpure & Xing 2014; Puechmaille 2016; Wang 2017). Studies on sampling bias sensitivity in DAPC are still required, but a recent study found that de novo clustering was inaccurate in situations with high migration rates (Miller et al. 2020). We did not explicitly quantify migration but the effects of this type of bias on our analyses cannot be ruled out. Furthermore, potential biases alluded to above that derive from using SNP data to reconstruct phylogenies also limit the extent to which these could be used to infer phylogeographical patterns.

Taxonomic implications

Upon reviewing the findings presented, it becomes evident that the phenotype-based species delimitation approach that has been traditionally applied to this complex fails to capture underlying patterns of genetic variation. This deficiency is particularly apparent in our phylogenetic analysis, highlighting the necessity for a revised taxonomic framework that better aligns with the evolutionary history of these populations.

The significant impact of geography on the genetic structure of these populations suggests that a more suitable approach would be to adopt an alternative taxonomy that recognizes this geographical structuring. Ideally, this taxonomic framework would recognize the continental clade and individual island-level clades as distinct taxonomic entities. However, such a solution must be cognizant of the results from our DAPC de novo clustering analysis, which indicate that certain populations may share genetic similarities (e.g. PR + Cuba and Jamaica + continental US; Fig. 5). Additionally, marginal support for the continental clade must be considered.

Added to the limited genetic variability observed in island-level populations, we consider that an ideal avenue to circumvent these taxonomic complexities would be to accept these island and continental populations as geographical races within one species. This approach aligns with Mayr's definition of a subspecies (Mayr 1942) and acknowledges the significant role of geography in shaping the genetic identity of these populations.

As a result, we introduce a revised classification in the Taxonomy section below which includes a new combination and the description of two new subspecies. Furthermore, the single nucleotide position within ITS1 that was found to vary between the two new subspecies is included as part of their diagnoses (see Taxonomy below). It is worth noting that the subspecies recognized herein account for the nomenclatural precedence of Cladonia sandstedei over C. subtenuis (Abbayes 1938, 1939). Further sampling efforts will be needed to determine if populations in Hispaniola, which were not included in our sampling, warrant taxonomic recognition.

The taxonomy proposed here is tentative. An alternative would be to continue accepting the continental population at the species level. For those preferring such a classification, the name Cladonia subtenuis (Abbayes) Mattick is available in a different circumscription as discussed above.

Taxonomy

Cladonia sandstedei Abbayes

J. Bot. 76, 349 (1938).—Cladina sandstedei (Abbayes) Ahti, Beih. Nova Hedwigia 79, 40 (1984); type: Jamaica, St Andrew, vicinity of Cinchona, 1490 m, 1896, W. Harris 10007 (G—lectotype; UCWI—isolectotype).

Cladonia rangiferina var. intricata Müll. Arg., Flora 69, 253 (1886); type: Jamaica, Gorton town, J. Hart 71 comm. Joshua 1885 (G—holotype).

Cladonia sandstedei subsp. sandstedei

Notes

According to our results, C. sandstedei subsp. sandstedei has a variable chemistry, with individuals either exhibiting the typical atranorin chemosyndrome or concurrently containing atranorin and usnic acid. It is geographically restricted to Jamaica, and therefore endemic to this island.

Samples from C. sandstedei subsp. sandstedei used in this study were collected in montane broadleaf forests in the Blue Mountains National Park in Portland Parish, specifically along the trail leading to the Blue Mountains Peak from Portland Gap. These specimens were found among bryophytes, growing directly on rocks or soil in open to partly shaded situations. Based on locality data in Ahti (2000), C. sandstedei subsp. sandstedei has a much broader distribution within this island. Some of these areas include different sites in the parish of St Andrew, particularly high-elevation forests (c. 1490 m) near the vicinity of Cinchona (type specimen) and other mid- to high-elevation (915–1500 m) areas in this parish (Ahti 2000). It is also expected to occur in other high-elevation areas (1220–2250 m) in St Thomas Parish (Ahti 2000).

Based on the taxonomy proposed here, definitive morphological features to distinguish C. sandstedei subsp. sandstedei from the other subspecies described below are still lacking. Individuals showing other combinations of the chemosyndromes associated with this subspecies have not been found but cannot be ruled out. It is worth noting that phylogenetic analysis placed the single C. sandstedei subsp. sandstedei specimen showing the typical atranorin chemosyndrome (Mercado-Díaz 3441 (UPR)) outside a strongly supported clade including specimens with both substances (Fig. 3). Since none of our population-level analyses provided a strong indication of divergence of this specimen from the rest, we treat it as belonging to this subspecies.

With the exception of our DAPC de novo clustering analysis, most of the analyses performed supported the separation of C. sandstedei subsp. sandstedei and C. sandstedei subsp. subtenuis. Additional sampling and extended genomic analysis will be needed to further explore the genetic similarities that this method suggested might exist between these populations.

Selected specimens examined

Jamaica: Cedar Valley: Portland Parish, Blue Mountains National Park, trail to Blue Mountain Peak from Portland Gap, 18°2ʹ41ʺN, 76°35ʹ24ʺW, 1934 m, 2018, J. A. Mercado-Díaz 3437 (UPR), 3441 (UPR); 18°2ʹ38ʺN, 76°35ʹ7ʺW, 2121 m, 2018, J. A. Mercado-Díaz 3446 (UPR); 18°2ʹ45ʺN, 76°35ʹ1ʺW, 2165 m, 2018, J. A. Mercado-Díaz 3460 (UPR).

Cladonia sandstedei subsp. subtenuis (Abbayes) Merc.-Díaz, Grewe & Lumbsch comb. nov.

MycoBank No.: MB 856207

Cladonia tenuis subsp. subtenuis Abbayes, Bull. Soc. Sci. Bretagne 16(fasc. hors, ser. 2), 106, 108 (1939).—Cladonia subtenuis (Abbayes) Mattick, Repert. Spec. Nov. Regni Veg. 49, 165 (1940); type: USA, Florida, Seminole Co., near Sanford, Rapp, Kryptog. exs. mus. Vindob. 3066 (H—lectotype; H, UPS, US—isolectotypes).

Cladonia subtenuis f. cinerea Ahti, Ann. Soc. Zool.-Bot. Fenn. ‘Vanamo’ 32, 69 (1961).—Cladina subtenuis (Abbayes) Hale & W. L. Culb., Bryologist 73, 510 (1970).—Cladina subtenuis f. cinerea (Ahti) Ahti, Beih. Nova Hedwigia 79, 41 (1984); type: USA, New York, Suffolk Co., Long Island, Promised Land, 1951, Latham 27630 (NYS—holotype).

Notes

Cladonia sandstedei subsp. subtenuis has a variable secondary chemistry, exhibiting all combinations of the atranorin and usnic acid chemosyndromes that have been previously linked to C. sandstedei s. lat. and the two infraspecific ‘formae’ associated with C. subtenuis s. lat. (Ahti 2000). It is so far known from tropical to cool temperate continental zones in eastern North America. No major morphological features distinguish C. sandstedei subsp. subtenuis from the other subspecies described here.

Specimens of Cladonia sandstedei subsp. subtenuis used in this study were collected in a diversity of ecosystems in the states of Tennessee, Georgia, Alabama and Florida. Elevation in these areas ranged from sea level to 512 m and included temperate deciduous forests near the Big South Fork National River and Recreation Area in Tennessee, Florida Rosemary Bald ecosystem and longleaf pine-wiregrass woodlands in Georgia, mixed pine/hardwood temperate forests in Alabama and mixed scrub forest with Palmetto palms and Myrtle Oak in Florida. At higher elevations, C. sandstedei subsp. subtenuis tends to grow directly on or among rocky outcrops, although it is commonly found on sandbanks or other sandy soils at lower elevations (Ahti 2000). As with the other subspecies, C. sandstedei subsp. subtenuis prefers open to partly shaded situations.

As demonstrated in this study, it is expected that additional C. subtenuis s. lat. individuals collected in the Caribbean islands will cluster under island-specific clades in subsequent phylogenetic studies. However, according to Ahti (2000), other continental populations of C. subtenuis s. lat. exist in areas outside of eastern North America (i.e. Guatemala, Nicaragua and Mexico). Additional studies will be needed to corroborate the phylogenetic position and taxonomic status of these populations. Ahti (2000) further highlights a single collection of C. sandstedei made in 1907 at Eight Mile Rock in Grand Bahama (Britton & Millspaugh 2606 (NY)). Assuming homoplasy in chemical signature and considering the short distance between this site and the continent (c. 77 miles) and the relatively young age of Grand Bahama's substrata (Pleistocene to Holocene), we suspect these populations will cluster with C. sandstedei subsp. subtenuis in subsequent phylogenetic analyses. Such a finding would effectively extend the distribution of C. sandstedei subsp. subtenuis to the insular Caribbean. Unfortunately, the areas surrounding Eight Mile Rock were subject to considerable damage by Hurricane Matthew in October 2016, thus the conservation status of these populations remains unknown. Furthermore, additional sampling and extended genomic analysis will be needed to clarify the potential conspecificity between C. sandstedei subsp. sandstedei and C. sandstedei subsp. subtenuis that was found in our DAPC de novo clustering analysis.

Selected specimens examined

USA: Alabama: Covington, roadside population along Eagle Rd (State Road 14), Conecuh National Forest, 31°9ʹ32ʺN, 86°36ʹ16ʺW, 68 m, 2019, J. A. Mercado-Díaz 3824, 3825, 3826, 3827, 3828, 3829, 3830, 3831, 3832 (all UPR); ibid., roadside population along road going to the Solon Dixon Forestry Education Center (Center Rd), 31°9ʹ4ʺN, 86°41ʹ54ʺW, 84 m, 2019, J. A. Mercado-Díaz 3821, 3822, 3823 (all UPR); Escambria, sandbank to the east of dorms parking lot, Solon Dixon Forestry Education Center, 31°9ʹ54ʺN, 86°42ʹ5ʺW, 88 m, 2019, J. A. Mercado-Díaz 3795, 3796, 3797, 3798, 3799, 3800, 3801, 3802, 3803, 3804, 3805, 3806 (all UPR); ibid., side of dirt road (Blue Pond Rd) leading to Blue Pond aka ‘the sinkhole’, 31°9ʹ42ʺN, 86°41ʹ52ʺW, 68 m, 2019, J. A. Mercado-Díaz 3807, 3808, 3809, 3810, 3811, 3812, 3813, 3814, 3815, 3816, 3817, 3818, 3819 (all UPR). Florida: Okaloosa, roadside populations on W. College Blvd, beside entrance to Twin Cities hospital Emergency Room, 30°32ʹ4ʺN, 86°29ʹ52ʺW, 12 m, 2019, J. A. Mercado-Díaz 3863, 3864, 3865, 3866, 3867, 3868, 3869, 3870, 3871 (all UPR); ibid., scrub forest at the Eglin Matterhorn Beach Access Point, 30°23ʹ24ʺN, 86°32ʹ41ʺW, 5 m, 2019, J. A. Mercado-Díaz 3858, 3859, 3860, 3862 (all UPR); ibid., surroundings of parking lot at the Gulf Islands National Seashore Okaloosa area, 30°23ʹ46ʺN, 86°34ʹ50ʺW, 1 m, 2019, J. A. Mercado-Díaz 3833, 3834, 3835, 3839, 3840, 3841 (all UPR); Walton, along nature trail in Henderson State Park, 30°23ʹ9ʺN, 86°26ʹ57ʺW, 5 m, 2019, J. A. Mercado-Díaz 3855, 3856, 3857 (all UPR); ibid., north of Henderson State Park parking lot, 30°23ʹ6ʺN, 86°26ʹ42ʺW, 5 m, 2019, J. A. Mercado-Díaz 3842, 3843, 3844, 3845, 3847, 3848 (all UPR); ibid., north of Henderson State Park parking lot, close to abandoned trail, 30°23ʹ7ʺN, 86°26ʹ33ʺW, 3 m, 2019, J. A. Mercado-Díaz 3850, 3851, 3852, 3853, 3854 (all UPR). Georgia: Coffee, Flat Rock Altamaha Formation Outcrop, longleaf pine-wiregrass sandhill (woodland), Broxton Rocks Preserve, 32°44ʹ4ʺN, 82°51ʹ19ʺW, 74 m, 2019, J. A. Mercado-Díaz, M. Hodges & S. Q. Beeching 3748, 3749, 3750, 3751, 3752, 3753 (all UPR); ibid., 32°44ʹ16ʺN, 82°51ʹ8ʺW, 63 m, 2019, J. A. Mercado-Díaz, M. Hodges & S. Q. Beeching 3754, 3755, 3756, 3757, 3758, 3759, 3760, 3761, 3762, 3763, 3764 (all UPR); ibid., roadside population along Rd 169 (Old River Rd), on Altamaha Formation Outcrop embedded in pine plantation, c. 0.75 miles north of entrance to Broxton Rocks Reserve, 32°44ʹ0ʺN, 82°52ʹ55ʺW, 86 m, 2019, J. A. Mercado-Díaz, M. Hodges & S. Q. Beeching 3787, 3788, 3789, 3790, 3791, 3792, 3793, 3794 (all UPR); Jasper, Charlie Elliot Wildlife Center, 33°27ʹ43ʺN, 83°44ʹ15ʺW, 226 m, 2005, J. A. Mercado-Díaz 2138 (UPR); Jeff Davis, roadside population along State Rd 107, on Altamaha Formation Outcrop embedded in pine plantation, Flat Tub WMA, 32°45ʹ42ʺN, 82°49ʹ47ʺW, 72 m, 2019, J. A. Mercado-Díaz 3767, 3768, 3769, 3770, 3771, 3772, 3773, 3774, 3775, 3776, 3777, 3778, 3779, 3780, 3781, 3782, 3783, 3784, 3785, 3786 (all UPR); Swainsboro, dirt road at intersection with trail entrance at hunter info. kiosk and Halls Bridge Rd, going north to Florida Rosemary Bald, 32°32ʹ4ʺN, 82°27ʹ32ʺW, 62 m, 2019, J. A. Mercado-Díaz 3736 (UPR); ibid., end of dirt road at McCleod Br. Rd and Robin Rd intersection, on loop at the end, east of Ohoopee River, 32°36ʹ12ʺN, 82°25ʹ45ʺW, 79 m, 2019, J. A. Mercado-Díaz 3700 (UPR), 3709 (UPR); Florida Rosemary Bald, besides trail, 32°32ʹ16ʺN, 82°27ʹ36ʺW, 64 m, 2019, J. A. Mercado-Díaz 3723, 3724, 3725, 3727, 3728, 3730, 3734 (all UPR); ibid., open area, roadside on US 221 N near milepost ‘Mile 1’, 32°28ʹ27ʺN, 82°29ʹ24ʺW, 75 m, 2019, J. A. Mercado-Díaz 3693, 3695, 3696, 3698 (all UPR); ibid., roadside population on Keens Crossing Rd (dirt road), west of Scott community, 32°33ʹ49ʺN, 82°43ʹ9ʺW, 86 m, 2019, J. A. Mercado-Díaz 3743, 3744, 3745 (all UPR); ibid., roadside population on Road 15 c. 2 miles north of city of Adrian, 32°33ʹ14ʺN, 82°36ʹ21ʺW, 95 m, 2019, J. A. Mercado-Díaz 3739, 3740, 3741, 3742 (all UPR). Tennessee: Scott, Big South Fork National River and Recreation Area, overlook near south arch, 36°32ʹ28ʺN, 84°44ʹ8ʺW, 512 m, 2018, J. A. Mercado-Díaz 3616 (UPR); ibid., Sawmill trail, after trail intersection along trail near an informal overlook, 36°31ʹ48ʺN, 84°45ʹ36ʺW, 512 m, 2018, J. A. Mercado-Díaz 3617 (UPR); ibid., Twin Arches trailhead, loop intersection with trail from parking, 36°32ʹ24ʺN, 84°44ʹ24ʺW, 512 m, 2018, J. A. Mercado-Díaz 3618 (UPR); Spencer, Big South Fork National River and Recreation Area, Sawmill trail, after trailhead intersection, right after first wooden bridge in trail, in camping area, 36°31ʹ12ʺN, 84°46ʹ12ʺW, 512 m, 2018, J. A. Mercado-Díaz 3613 (UPR); ibid., Gorge outlook trail, along Rocky Point Overlook trail in Fall Creek Falls, 35°39ʹ59ʺN, 85°21ʹ15ʺW, 514 m, 2018, J. A. Mercado-Díaz 3609 (UPR), 3610 (UPR).

Cladonia sandstedei subsp. cubana Merc.-Díaz, Rivera-Queralta & Motito-Marín subsp. nov.

MycoBank No.: MB 856208

Subspecies cubana can be distinguished from the rest of the subspecies within the C. sandstedei complex by the presence of thymine in the position within parentheses in the following ITS1 region section: RCCCY(T)AGCGTT. It further differs in its distribution since it is found only in Cuba.

Type: Cuba, Santiago de Cuba, Sierra Maestra, Paisaje Natural Protegido La Gran Piedra, 21 km E of Santiago, Cafetal La Isabelica, 20°00ʹ19ʺN, 75°37ʹ06ʺW, 1110 m, secondary vegetation surrounded by montane rainforest, on the ground, 2019, Lücking & Moncada 46737 (BSC—holotype; B, HAJB, UPR—isotypes).

Thallus fruticose, well developed, consisting of loosely to highly branched podetia that form irregular to subglobose, cushion-shaped colonies. Podetia 4–8 cm tall, 0.3–1 mm thick, whitish to light greenish grey when fresh, whitish grey to light tan in herbarium conditions, darkening towards the substratum but not becoming melanotic, often browned at tips; branching type anisotomous to subisotomous dichotomy, rarely trichotomous or tetrachotomous; main stems occasionally distinguishable towards the base but most often not well differentiated, absent from juvenile podetia; axils rarely perforated; mature internodes 2–5 mm long, 0.5–1 mm wide, ultimate apical branchlets erect, divaricate, acute to deflected. Podetial surface ecorticate, arachnoid, dull, of smooth appearance towards tips but becoming slightly rugulose to verrucose towards the base where it tends to disintegrate; soredia and squamules lacking. Podetial wall 80–150 μm; ectal layer 30–100 μm; stereome c. 45–85 μm.

Apothecia brown, darkening with age, stipitate, infrequent, but sometimes abundant in portions of the thallus, 0.1–0.3 wide, spores simple, fusiform, 7.5–11 × 2.5–3.25 μm.

Conidiomata not seen.

Etymology

The epithet reflects the restricted distribution of this subspecies to the island of Cuba.

Secondary chemistry

Atranorin and fumarprotocetraric acid complex as major products. P+ red, K+ yellow, C−, UV−. Refer to Ahti (2000) for additional information on minor products.

Ecology and distribution

Cladonia sandstedei subsp. cubana is restricted to the island of Cuba. Specimens examined for this study derived from previously known populations at the Paisaje Natural Protegido La Gran Piedra and from the Parque Nacional Turquino, both in the Santiago de Cuba Province. Within the Paisaje Natural Protegido La Gran Piedra, C. sandstedei subsp. cubana was collected in open to partly shaded areas in the montane rainforests associated with Pico Mogote (c. 1000 m elev.). This area is characterized by mean annual temperatures ranging from 16–18 °C and a mean annual precipitation ranging from 1300–1900 mm (Birdlife International 2024). This subspecies was also collected in semi-open areas along the trail leading to Pico Turquino at the Parque Nacional Turquino, specifically near the Pico Turquino summit (1922 m elev.) and along the section leading to Loma Redonda from the Aguada de Joaquin camp. Moist broadleaf forests near Pico Turquino are characterized by mean annual temperatures ranging from 5–16 °C and relatively high precipitation (1500–1700 mm year⁻¹) (Birdlife International 2024).

While not visited for this study, Cladonia sandstedei subsp. cubana is also expected to occur in other C. sandstedei s. lat. sites previously highlighted by Ahti (2000) for Cuba. These sites are environmentally and ecologically diverse and include areas in the province of Guantánamo, near Monte Cristo and Cupeyal, the province of Pinar del Río, near Viñales, and La Isla de La Juventud, near Los Indios.

Remarks

No major morphological features distinguishing C. sandstedei subsp. cubana from the other subspecies described here have been identified. Compared to C. sandstedei subsp. landroniana (see below), C. sandstedei subsp. cubana showed a slight tendency for podetial tips to become more frequently browned, but the limited number of Cuban individuals evaluated suggest this might represent a sampling artifact. Individuals showing other combinations of the chemosyndromes associated with this subspecies have not been found but cannot be ruled out. As suggested by our fineRADstructure analysis, it is plausible that a relatively recent genetic divergence process might explain the clustering of C. sandstedei subsp. cubana and C. sandstedei subsp. landroniana (below) in our DAPC de novo clustering analysis. Further scrutiny using an extended species sampling and divergence time estimates will be needed to further clarify these issues.

Selected specimens examined

Cuba: Santiago de Cuba: Paisaje Natural Protegido La Gran Piedra, near entrance to Pico Mogote trail, 2016, J. A. Mercado-Díaz, Y. Rivera-Queralta & Y. Tomé 2557 (BSC, UPR); Turquino National Park, segment between Loma Redonda and the Aguada de Joaquín camp towards Pico Turquino summit, 2016, J. A. Mercado-Díaz & Y. Rivera-Queralta 2712 (BSC, UPR); near Pico Turquino summit, 2016, J. A. Mercado-Díaz & Y. Rivera-Queralta 2793 (BSC, UPR).

Cladonia sandstedei subsp. landroniana Merc.-Díaz, Grewe & Lumbsch subsp. nov.

MycoBank No.: MB 856209

Subspecies landroniana can be distinguished from subspecies cubana by the presence of cytosine in the position within parentheses in the following ITS1 region section: RCCCY(C)AGCGTT. It is only distinguishable from the rest of the subspecies within the C. sandstedei complex by its distribution, since it is found only in Puerto Rico.

Type: USA, Puerto Rico, Maricao, Bosque Estatal de Maricao, short trail that leads to small pine plantation along PR-120, 18°9ʹ48ʺN, 66°59ʹ39ʺW, 653 m, secondary vegetation surrounded by montane wet forest, on the ground above serpentine-derived soils, 2017, J. A. Mercado-Díaz 3686 (UPR—holotype).

Thallus fruticose, well developed, consisting of loosely to highly branched podetia that form irregular to subglobose, cushion-shaped colonies. Podetia 4–10 cm tall, 0.3–1 mm thick, whitish to light greenish grey when fresh, whitish grey to light tan in herbarium conditions, darkening towards the substratum but not becoming melanotic, weakly to moderately browned at tips, but often of the same colour as branches; branching type anisotomous to subisotomous dichotomy, rarely trichotomous or tetrachotomous; main stems occasionally distinguishable towards the base but most often not well differentiated, absent from juvenile podetia; axils rarely perforated; mature internodes 2–5 mm long, 0.5–1 mm wide, ultimate apical branchlets acute, divaricate, erect but occasionally becoming deflected with age. Podetial surface ecorticate, arachnoid, dull, of smooth appearance towards tips but becoming slightly rugulose to verrucose towards the base where it tends to disintegrate; soredia and squamules lacking. Podetial wall c. 80–120 μm; ectal layer 25–90 μm; stereome c. 50–90 μm.

Apothecia and conidiomata not seen.

Etymology

This subspecies is dedicated to the late Dr Ismael Landrón-Concepción for his contributions to the study of lichens in Puerto Rico. He was the first Puerto Rican with formal training in lichenology.

Secondary chemistry

Atranorin and fumarprotocetraric acid complex as major products. P+ red, K+ yellow, C−, UV−. Refer to Ahti (2000) for additional information on minor products.

Ecology and distribution

Cladonia sandstedei subsp. landroniana is restricted to Puerto Rico and has a disjunct distribution within this island. Areas where C. sandstedei subsp. landroniana populations are found differ markedly in habitat characteristics. One of these areas is the central-west mountains, within the Maricao State Forest. This is a high-elevation forest, with relatively high humidity and serpentinite-derived soils that are generally nutrient limited (Ricart-Pujals & Padrón-Vélez 2010). The second habitat is near the north-central coast in the Laguna Tortuguero area in the Vega Baja and Manatí municipalities. Laguna Tortuguero features relatively mild to high coastal temperatures (24–31 °C), variable humidity and precipitation ranging between 1000–2000 mm year⁻¹. Areas of Laguna Tortuguero where C. sandstedei subsp. landroniana is found are characterized by the presence of white silicious sands that are also relatively poor in nutrient content (Lugo et al. 2001). Unfortunately, efforts carried out to locate C. sandstedei subsp. landroniana in areas with habitat characteristics similar to Laguna Tortuguero (e.g. Ciénaga Prieta, Vega Alta) have failed to yield new populations of this subspecies.

Remarks

As with the other subspecies, there are no clear morphological attributes that unambiguously separate C. sandstedei subsp. landroniana from the other subspecies in the complex. One aspect worth mentioning in this regard is that we failed to detect any reproductive structures in the specimens evaluated, a finding potentially linked to environmental factors specific to the sites where these individuals were collected.

Phylogenetic analysis along with PCA and fineRADstructure results provided evidence that as a geographical unit, C. sandstedei subsp. landroniana is genetically distinct from other C. sandstedei subspecies within this region. In fact, there is a further indication of potential genetic subdivision within this island as phylogenetic analysis placed the Vega Baja population in a strongly supported monophyletic clade. However, since the Maricao population was associated with a paraphyletic grade and considering that the rest of the population-level analyses failed to separate Maricao from Vega Baja, we opted to refrain from providing further subspecies recognition. In this vein, Ahti (2000) documented a single occurrence of ‘Cladina subtenuis f. subtenuis’ in Guadeloupe which suggests that C. sandstedei individuals with chemical make-ups other than the typical atranorin chemosyndrome exist beyond eastern Puerto Rico. Based on the patterns observed in this work, we suspect that irrespective of chemical signature, populations in the Lesser Antilles will be phylogenetically closer to C. sandstedei subsp. landroniana. Resolving the taxonomic status of these populations is contingent on extended sampling in these islands followed by additional phylogenetic analysis.

Based on the published number of specimens and the vouchers recorded in this study, Puerto Rico stands out among other islands in the Caribbean in terms of surveying efforts associated with these subspecies. It is therefore unlikely that additional collecting efforts for C. sandstedei subsp. landroniana on this island will yield individuals with chemosyndromes other than the typical atranorin chemosyndrome that has been linked to this subspecies.

Selected specimens examined

USA: Puerto Rico: Manatí, Reserva Natural Laguna Tortuguero, 18°27ʹ35ʺN, 66°26ʹ37ʺW, 2 m, 2017, J. A. Mercado-Díaz, 3322, 3323, 3324, 3325, 3326 (all UPR); Maricao, Bosque Estatal de Maricao, PR-120, Km 17, Hm. 5, 18°9ʹ0ʺN, 66°59ʹ0ʺW, 800 m, 2017, J. A. Mercado-Díaz 280 (UPR); ibid., short abandoned trail to the left of old quarry on PR-120, 18°9ʹ27ʺN, 66°59ʹ54ʺW, 736 m, 2020, J. A. Mercado-Díaz 4117, 4118, 4119 (all UPR); ibid., short trail that leads to small pine plantation along PR-120, 18°9ʹ48ʺN, 66°59ʹ39ʺW, 653 m, 2017, J. A. Mercado-Díaz 3313, 3315, 3316 (all UPR); ibid., 2018, J. A. Mercado-Díaz 3685, 3686, 3687, 3688 (all UPR); ibid., 18°9ʹ48ʺN, 66°59ʹ40ʺW, 655 m, 2020, J. A. Mercado-Díaz 4120 (UPR); ibid., 18°9ʹ49ʺN, 66°59ʹ39ʺW, 655 m, 2020, J. A. Mercado-Díaz 4121, 4122, 4123, 4124 (all UPR); Vega Baja, Reserva Natural Laguna Tortuguero, 18°27ʹ0ʺN, 66°26ʹ0ʺW, 4 m, 2003, J. A. Mercado-Díaz 1837 (UPR), 1839 (UPR); 18°27ʹ26ʺN, 66°26ʹ12ʺW, 0 m, 2017, J. A. Mercado-Díaz 3334 (UPR), 3335 (UPR); ibid., 18°27ʹ25ʺN, 66°26ʹ11ʺW, 0 m, 2020, J. A. Mercado-Díaz 4099 (UPR); ibid., 18°27ʹ24ʺN, 66°26ʹ12ʺW, 0 m, 2020, J. A. Mercado-Díaz 4100 (UPR); 18°27ʹ25ʺN, 66°26ʹ15ʺW, 0 m, 2020, J. A. Mercado-Díaz 4102 (UPR); ibid., 18°27ʹ23ʺN, 66°26ʹ13ʺW, 0 m, 2020, J. A. Mercado-Díaz 4104 (UPR); ibid., 18°27ʹ25ʺN, 66°26ʹ12ʺW, 0 m, 2020, J. A. Mercado-Díaz 4105 (UPR); ibid., 18°27ʹ26ʺN, 66°26ʹ11ʺW, 0 m, 2020, J. A. Mercado-Díaz 4106 (UPR); ibid., 18°27ʹ27ʺN, 66°26ʹ10ʺW, 0 m, 2020, J. A. Mercado-Díaz 4107 (UPR); ibid., 18°27ʹ27ʺN, 66°26ʹ9ʺW, 0 m, 2020, J. A. Mercado-Díaz 4108 (UPR); ibid., 18°27ʹ28ʺN, 66°26ʹ7ʺW, 0 m, 2020, J. A. Mercado-Díaz 4112 (UPR); ibid., 18°27ʹ29ʺN, 66°26ʹ7ʺW, 0 m, 2020, J. A. Mercado-Díaz 4114 (UPR); ibid., 18°27ʹ29ʺN, 66°26ʹ8ʺW, 0 m, 2020, J. A. Mercado-Díaz 4116 (UPR).

Conclusion

Patterns of genetic variation in C. sandstedei and C. subtenuis were highly variable in the region studied. We found no clear distinction between C. sandstedei versus C. subtenuis populations within the continent, which suggests that the current phenotype-based delimitation of these taxa is arbitrary. In contrast, stronger patterns of geographical structuring were found in the Caribbean. Island-specific clades that were clearly separated from a continental clade were recovered during phylogenetic analyses, irrespective of current taxonomic placement. Similarly, most population-level analyses performed suggested that all island-level and continental populations associated with these species represented genetically distinct units. Considering challenges associated with conserving the current taxonomic concept of these species, four subspecies were designated to account for these findings. A potential taxonomic subdivision of populations within Puerto Rico is plausible, but the evidence was inconclusive. Extended sampling efforts, particularly in Jamaica, Hispaniola, Cuba, central/southern Florida, and other continental areas where C. subtenuis s. lat. and C. sandstedei s. lat. have been reported will be needed for further refinement of the genetic and taxonomic boundaries between these subspecies. Our work illustrates the importance of using phylogenetic and population genetic frameworks to assess the phylogeography of poorly understood, sympatrically distributed species of phenotypically similar lichens. Adherence to these frameworks leads to more efficient ways to delimit species, which ultimately allows for better characterization of the diversity present in our ecosystems.