Yakushima (Yaku Island) is a roughly circular island with a circumference of approximately 130 km and is located approximately 60 km south of the main island of Kyushu. Since the initial taxonomic review of 45 species of vascular plants endemic to Yakushima (Yahara et al. 1987), six additional new species have been described in the same region (Yahara and Tsukaya 2008; Katsuyama 2009; Chen et al. 2014; Hori et al. 2015; Suetsugu and Fukunaga 2016; Suetsugu et al. 2016). Furthermore, an unknown taxon of Cardamine was also discovered in the mountainous area of Yakushima (Kudoh 2017). Between July 19 and July 24, 2020, a field trip was made to the mountainous area of Yakushima in order to collect this unknown species of Cardamine, whereby we serendipitously discovered a plant with cleistogamous flowers of the unknown taxon of Stellaria. Later, fruiting specimens of this taxon were collected on September 4, 2020, and a specimen with a chasmogamous flower was collected on May 3, 2021 in the same area by K. Fuse. In addition, we collected samples of S. alsine subsp. alsine, S. aquatica (L.) Scop., S. diversiflora Maxim., S. media (L.) Vill., S. monosperma Buch.-Ham. ex D. Don, and S. neglecta Weihe (Table 1), and examined phylogenetic relationships.

Total DNA was extracted from dried leaves using the CTAB method (Doyle and Doyle 1990). De novo SNP discovery was performed using MIG-seq (Suyama and Matsuki 2015). Based on the methodology described by Suyama and Matsuki (2015), a MIG-seq library was prepared via a two-step PCR amplification process with minor modifications, namely the annealing temperature of the first PCR was altered from 48 °C to 38 °C. Subsequently, the second PCR products were purified in the size range of 300–800 bp and sequenced on an Illumina MiSeq platform (Illumina, San Diego, CA, USA) using a MiSeq Reagent Kit v3 (150 cycle, Illumina). Sequencing of the first 17 bases of reads 1 and 2 (SSR primer regions and anchors) was bypassed using the ‘DarkCycle’ function of the MiSeq platform. Additionally, low-quality reads and extremely short reads containing adapter sequences were removed using the Trimmomatic 0.39 software (Bolger et al. 2014). The Stacks 2.41 pipeline (Catchen et al. 2013; Rochette et al. 2019) was used to obtain individual genotypes with the following parameters: minimum depth of coverage required to create a stack (m) = 3, maximum distance between stacks (M) = 2, and maximum mismatches between loci when building the catalog (n) = 2. Three different filtering criteria were applied for quality control of the SNP data. First, any SNP site where one of two alleles had less than three counts was filtered out due to the difficulty in distinguishing polymorphisms from sequencing errors that arise when the minor allele count of SNPs is too low (Roesti et al. 2012). Second, loci containing SNPs with high heterozygosity (Ho ≥ 0.6) were removed as the excess heterozygosity may have resulted from artifactual loci constructed from several paralogous genomic regions. Third, the SNPs retained by three or more samples were included in the SNP dataset.

Maximum likelihood phylogeny based on SNPs was inferred using the RAxML 8.2.10 software (Stamatakis 2014). A GTRCAT model was applied during this process and 1,000 replicates of parallel tree search bootstrapping were performed.

Three chloroplast and two nuclear genomic regions were sequenced using next-generation DNA sequencing. In this regard, rbcL, trnL intron, psbA-trnH, ITS1, and ITS2 were initially simultaneously amplified using the Multiplex PCR Assay Kit Ver. 2 (Takara Bio, Kusatsu, Japan) (first PCR). The first set of primers consisted of tail sequences and locus-specific primers (Suyama et al. 2022). Subsequently, the first PCR products were purified and used for the second PCR. The second PCR was conducted using primer pairs, including tail sequences, adapter sequences for Illumina sequencing, and the index sequence to identify each individual sample. In this step, the second PCR product from each sample was mixed and sequenced on the Illumina MiSeq platform using a MiSeq Reagent Nano Kit v2 (500 cycle, Illumina). The sequencing of the first three bases of reads 1 and 2 (anchor region for the 2^nd PCR primer) was bypassed using the ‘DarkCycle’ option of the MiSeq system. Both ends of the fragments and index sequences were read by paired-end sequences (reads 1 and 2) and index sequencing. The number of bases per read was 251 bases for read 1 and 251 bases for read 2.

The sequences of the five regions were determined using the Claident pipeline (Tanabe and Toju 2013, http://www.claident.org/, Tanabe, A.S., Claident, Date of access: 05/01/2021). The raw MiSeq BCL data were first converted into FASTQ data using the BCL2FASTQ program provided by Illumina, followed by demultiplexing of the raw FASTQ data based on index and primer sequences, using the clsplitseq program in Claident. Subsequent analysis of the data was performed per region per individual. In ITS1 and ITS2, the pair-end reads were merged since reads 1 and 2 were overlapping. In rbcL, trnL intron, and psbA-trnH, reads 1 and 2 were independently analyzed because the length of the sequenced reads was too short to allow for merging. Additionally, the low-quality 3' tail was trimmed, and the low-quality sequences were filtered out using the clfilterseq program. The noisy sequences were removed using the clcleanseqv program. Finally, the remaining reads were clustered with a cut-off sequence similarity of 99%. An OTU that had the most observed reads within the individual was treated as a representative OTU sequence.

Multiple alignments were performed using the MAFFT 7.313 program (Katoh and Standley 2013), and alignment columns containing gaps were trimmed using a heuristic selection method based on the similarity statistics of trimAl 1.4.rev15 (Capella-Gutierrez et al. 2009). The Kakusan 4.0 software (Tanabe 2011) was used to find suitable nucleotide substitution models and partitioning strategies for the nucleotide datasets by independently running the chloroplast and nuclear genomic regions through this program. The AICc criterion (Sugiura 1978) was used to compare the Nonpartitioned, the Partitioned_Equal_Mean_Rate, as well as the Separate models. For chloroplast genomic regions, the Partitioned_Equal_Mean_Rate model (GTR + Γ), which assumes an equal rate of nucleotide substitutions across arbitrarily specified partitions, proved optimal. Contrastingly, for nuclear genomic regions, the Nonpartitioned model (GTR + Γ) proved optimal. Maximum likelihood phylogenies were further inferred using RAxML 8.2.10 (Stamatakis 2014), whereby 1,000 replicates of parallel tree search bootstrapping were conducted.

A total of 15,551,282 raw reads (1,196,252 ± 87,288 reads per sample) were obtained via MIG-seq, of which 14,923,278 reads (1,147,944 ± 84,758 reads per sample) remained after quality control. Following the de novo SNP detection and filtering, a dataset comprised of 881 SNPs from 703 loci was obtained. Three samples – Stmon-JPN2119, Stmon-JPN2998, and Stdiv-JPN1352 – contained more than 90% of the miscalled SNPs; thus, these samples were eliminated from further SNP analysis.

Fewer shared loci were observed between S. alsine and S. diversiflora, as well as between S. alsine and S. monosperma, resulting in the exclusion of S. diversiflora and S. monosperma from the phylogenetic reconstruction. In the tree obtained (Fig. 1), the sister relationship between S. alsine subsp. alsine and subsp. nana was supported by a 100% bootstrap value. None of the three species included in the phylogenetic reconstruction, i.e., S. aquatica, S. media, and S. neglecta, are directly related to S. alsine.

Figure 1. 

A phylogenetic tree of Stellaria alsine and its relatives reconstructed using MIG-seq data. Stals: S. alsine subsp. alsine. Staqu: S. aquatica. Stmed: S. media. Stnan: S. alsine subsp. nana. Stneg: S. neglecta. Bootstrap values are shown on the nodes. Branch length represents the average number of substitutions per SNP site.

A total of 64,378 reads (4,952 ± 387 reads per sample, ITS1) and 68,362 reads (5,258 ± 357 reads per sample, ITS2) were obtained. After the gaps were trimmed, the total length of the remaining sequences was 662 bp. The monophyly of S. alsine subsp. alsine and subsp. nana was supported by a 100% bootstrap value (Fig. 2). Moreover, none of the five other species were directly related to S. alsine. Additionally, the sister relationship between S. diversiflora and S. monosperma was supported by a 92% bootstrap value.

Figure 2. 

A phylogenetic tree of Stellaria alsine and its relatives reconstructed using ITS sequences. Stals: S. alsine subsp. alsine. Staqu: S. aquatica. Stdiv: S. diversiflora. Stmed: S. media. Stmon: S. monosperma. Stnan: S. alsine subsp. nana. Stneg: S. neglecta. Bootstrap values are shown on the nodes. Branch length represents the average number of substitutions per site.

A total of 23,206 reads (1,785 ± 129 reads per sample, rbcL), 2,142 reads (165 ± 24 reads per sample, trnL intron), and 55,274 reads (4,252 ± 569 reads per sample, psbA-trnH) were obtained. After the gaps were trimmed, the total length of the remaining sequences was 1,467 bp. The monophyly of S. alsine subsp. alsine and subsp. nana was supported by a 100% bootstrap value. As observed in the case of ITS, none of the five other species was directly related to S. alsine (Fig. 3). Furthermore, the sister relationship between S. diversiflora and S. monosperma was supported by a 100% bootstrap value.

Figure 3. 

A phylogenetic tree of Stellaria alsine and its relatives reconstructed using cpDNA sequences. Bootstrap values are shown on the nodes.