The organellar genomes of Pellidae (Marchantiophyta): the evidence of cryptic speciation, conflicting phylogenies and extraordinary reduction of mitogenomes in simple thalloid liverwort lineage

Organellar genomes of liverworts are considered as one of the most stable among plants, with rare events of gene loss and structural rearrangements. However, not all lineages of liverworts are equally explored in the field of organellar genomics, and subclass Pellidae is one of the less known. Hybrid assembly, using both short- and long-read technologies enabled the assembly of repeat-rich mitogenomes of Pellia and Apopellia revealing extraordinary reduction of length in the latter which impacts only intergenic spacers. The mitogenomes of Apopellia were revealed to be the smallest among all known liverworts—109 k bp, despite retaining all introns. The study also showed the loss of one tRNA gene in Apopellia mitogenome, although it had no impact on the codon usage pattern of mitochondrial protein coding genes. Moreover, it was revealed that Apopellia and Pellia differ in codon usage by plastome CDSs, despite identical tRNA gene content. Molecular identification of species is especially important where traditional taxonomic methods fail, especially within Pellidae where cryptic speciation is well recognized. The simple morphology of these species and a tendency towards environmental plasticity make them complicated in identification. Application of super-barcodes, based on complete mitochondrial or plastid genomes sequences enable identification of all cryptic lineages within Apopellia and Pellia genera, however in some particular cases, mitogenomes were more efficient in species delimitation than plastomes.

with nine chromosomes were deliminated on either N or S based on presence different alleles in peroxidase locus. The samples with 18 chromosomes and peroxidase band pattern specific for Pellia epiphylla N were classified as P. borealis. Specimens with poor band staining or without confirmed chromosome numbers were excluded from analysis. For delimitation of Apopellia cryptic species intron-exon splice markers (ISJ) were used 24 according to previously published protocol 40 . Experimental research and field studies on plants (either cultivated or wild), including the collection of plant material, complied with relevant institutional, national, and international guidelines and legislation. DNA extraction. DNA from the fresh samples was isolated by the modified CTAB procedure. The liquid nitrogen-ground thalii were thoroughly mixed with 3 ml preheated CTAB isolation buffer (2% CTAB, 100 mM Tris-HCl, pH 8.0, 20 mM EDTA, 1.4 M NaCl and 2% β-mercaptoethanol) and incubated at 55 °C for 1 h. After three chloroform extractions, the DNA was precipitated and dissolved in sterile, deionized H 2 O. Total genomic DNA from herbarium specimens was extracted using ZR Plant/Seed DNA MiniPrep™ kit (Zymo Research Corp., Irvine, CA, USA) according to the manufacturer's protocol. The purity of DNA samples was assessed spectrophotometrically using Cary 60 spectrophotometer (Agilent). DNA quantity was estimated using the Qubit fluorometer and Qubit™ dsDNA BR Assay Kit (Invitrogen, Carsbad, NM, USA). DNA quality was checked by electrophoresis in 0.5% agarose gel stained with Euryx Simple Safe (Eurx, Gdańsk, Poland). The extracted DNA prior to long-read sequencing was carefully examined using Tapestation Genomic DNA kit (Agilent) and additionally cleaned using Genomic DNA Clean and Concentrator kit (Zymo, Irvine, USA).
Library preparation, sequencing and assembly. The genomic libraries were constructed with DNA TrueSeq (Illumina) and were sequenced using HiSeqX (Illumina) to generate 150 bp paired-end reads at Macrogen Inc. (Seoul, Korea) with 350 bp insert in size between paired-ends. Afterwards, sequencing reads were cleaned by removing the adaptor sequences and low quality reads with Trimmomatic v0.36 41 . The long-read libraries were constructed using Ligation Sequencing Kit SQK-LSK109 (Oxford Nanopore Technologies) and NEBNext ® Companion Module for Oxford Nanopore Technologies ® Ligation Sequencing (New England Biolabs) according to manufacturer's protocol and sequenced using MinION MK1B portable device (ONT) and R.9.4.1 Flow Cell (ONT). The Flow Cell was prepared for sequencing with Flow Cell Priming Kit EXP-FLP002 (ONT). Sequence reads were basecalled using high-accuracy (HAC) mode of Guppy 6.0.7 basecaller.
The previously sequenced Pellidae genomes, including plastome of Apopellia endiviifolia (NC_019628) and mitogenomes of Fossombronia cristula (MK645818) and Makinoa crispata (MK230958) couldn't be used as direct references due to low sequence similarity, especially considering Pellia sensu stricto genus. Therefore to assemble Pellia reference genomes we applied previously published approaches 42 . To assemble the Fossombronia crispula plastome SSR8246023 reads archive was used. The specimen details including GenBank accession numbers of newly assembled genomes are given in Supplementary Table S1.
The reference mitochondrial genomes for Pellia were assembled using Flye 2.8 43 with default settings to overcome long regions of short tandem repeats. The remaining Pellia organellar genomes were assembled using NOVOplasty 4.3.1 44 . The GeSeq software was used for annotation of organellar genomes 45 and genome maps were drawn using OGDRaw web server 46 . The junctions between chloroplast inverted repeats and single copy regions were visualised using IRscope software 47 .
Complete chloroplast and mitogenome genome sequences of 75 specimens, including 17 from Pellidae were used for phylogenetic analysis. Organellar genomes were initially aligned using MAFFT and ambiguously aligned regions were trimmed by Gblocks 0.91 48 . However, due to significant evolutionary divergence, reliable alignment wasn't possible for most plastid intergenic and intronic regions, therefore for the main phylogenetic analyses only protein-coding genes (PCG) sequences were selected. The mitogenome dataset was also reduced to PCG, due to the loss of most intergenic spacers in genus Apopellia and short tandem repeat richness of Pellia intergenic spacers.
Phylogenetic analysis was carried out using the maximum likelihood (ML). The optimal model for the plastome dataset was identified as GTR + F + I + G4 by ModelFinder based on the Bayesian information criterion (BIC). ML analysis was performed using IQ-tree 49 . Mitochondrion and plastome trees were compared using cophylo function of phytools 1.0-3 R package 50 .
In order to explore plastome and mitogenome gene-tree conflict across liverworts phylogeny the single gene trees were calculated for each plastid and mitochondrial PCG.
The consistency among single genes was evaluated by the maximum-likelihood (ML) analysis for every locus that contained parsimony informative characters using the RAxML v7.2.3 plugin 51 for Geneious Prime 2022 with the GTR + I + G model and bootstraps estimated using 1000 replicates. Potential conflicts among genes were visualised by constructing a supernetwork using SplitsTree4 52 according to the approach given by Liu et al. 37 , except the number of replicates was increased to 1000 and the software used for creating consensus trees was Geneious Prime 2022. The supernetworks were created separately for chloroplast and mitochondrial datasets.
Genetic variation, barcoding and molecular species delimitation. Comparative analysis of chloroplast genomes was carried out in the Spider 1.1-5 program 53 based on inter-and intraspecific distances that were calculated using the Kimura 2-parameter model (K2P) of nucleotide substitution. Due to reasons mentioned above Apopellia and Pellia datasets were analysed separately. Barcoding analyses of entire Pellia/Apopellia organellar genomes and their 500 bp-long fragments generated by sliding window (step by 100 bp) were made in Spider. The discrete molecular diagnostic characters (MDCs) for each species were calculated according to the Jörger and Schrödl 54 approach using FASTACHAR software 55 . Nucleotide diversity analyses were carried out using the PopGenome package 56  www.nature.com/scientificreports/ among Apopellia and Pellia species was calculated using simplot function of ggmsa package 57 , with sliding windows parameters as in previous analyses and accession numbers OL654073 (cpDNA), OQ236455 (mtDNA) and OQ280821 (cpDNA), OQ236469 (mtDNA) as references for Apopellia and Pellia respectively. Next to species delimitation analysis in Spider, ASAP: assemble species by automatic partitioning 58 was applied. This method generates ten partitions, computing for each partition the asap-score, which is the average of two components: probability of panmixia (p-value) and barcoding gap width. The lower the score, the better the partition for the dataset. As input data alignments of plastomes and mitogenomes were used separately in each analysis.
Codon usage analysis. The relative synonymous codon usage (RSCU) indicator was calculated based on the number of times a particular codon is observed, relative to the number of times that the codon would be observed assuming a uniform synonymous codon usage 59 . The RSCU values were calculated using the seqinr 4.2 R package 60 and visualised using ggplot2 61 .
To calculate difference in codon usage bias between Apopellia and Pellia protein-coding genes we applied the Measure Independent of Length and Composition (MILC) with hard filtering option to remove sequences shorter than 80 AA, as implemented in coRdon 1.15.0 package 62,63 . Differences in codon usage (CU) between genera were visualised using intraBplot function of coRdon package 63 .

Results
Characteristics of Pelliidae chloroplast genomes. Newly sequenced plastid genomes of Pelliidae are circular molecules consisting of regions typical for land plants. The total length of assembled plastomes ranged from 115,878 bp in Pellia neesiana ( Fig. 1) to 124,598 bp in Moerckia flotoviana. The GC content ranged from 35.8% in water form of A. endiviifolia to 41.6% in F. foveolata.
One hundred twenty-two unique genes (taking into account only one copy of inverted repeat regions) were identified in the plastomes of Apopellia and Moerckia: 81 protein-coding genes, four ribosomal RNAs, 31 transfer RNAs and six ycf genes of an indeterminate function. The plastomes of Fossombronia and Pellia genera lack cysA and cysT genes, additionally in the Fossombronia, clusters of trnS(GCU)-psbI-trnS(GCA) genes are missing. The gene order seems to be conserved in Pelliidae, with exception of translocation of rps12-rps7-ndhB to the minus strand of LSC/IR boundary.
The values within IR regions were almost five times lower than within single copy LSC and SSC parts of the plastome. In the both genera the single copy regions near LSC-IRb and IRb-SSC regions were also highly variable, especially in the case of SSC, where ndhF gene is located. The second peak of nucleotide diversity within the SSC region was within the ycf1 gene. However, the gene near the SSC-IRa border, chlL, was among less variable regions of Pellia and Apopellia plastomes.
The patterns of nucleotide diversity within Apopellia and Pellia genera usually overlap with few additional hotspots of the former found in the rpoC2 gene and igs between atpI-atpH, psbM-trnS, psbA-trnH and trnV-atpE (Fig. 2). Lower variation in Apopellia than Pellia was revealed at 5' end of ycf2, across whole accD, rpl22 and chIL genes.
The analysis of sequence similarity within the same 500 bp window revealed some regions of intraspecific variation, which was generally higher in Pellia than Apopellia (Fig. 2). In the case of the latter slight variation at cryptic species level was found within the ndhF gene and the trnG(UCC)-ycf12 spacer of the typical form. The mentioned spacer was also variable within P. epiphylla N and S. These two species shared infraspecific variation within few other plastome regions including ndhB-ycf66, atpH-atpF, ycf2-cysA, ndhC-trnV and ndhH-ycf1b. Additionally, plastomes of P. epiphylla S revealed variation in the ycf2 gene, ycf3-rps4 region and between rrn12 and rrn23 genes of IR.

Characteristics of Pelliidae mitochondrial genome. The newly assembled mitochondrial genomes
have almost identical gene content and order, as previously published (MK230958, MK749462 and MK230936) with exception of the trnR-UCG loss in the Apopellia genus (Fig. 3a). The mitogenomes of Pellidae ranged from 108,928 bp in typical form of Apopellia endiviifolia to 179,862 bp in Moerckia flotoviana (Table S1), but the length of the molecule doesn't correlate with intron content (Fig. 3b).
The species of Fossombroniales order lack both introns of the atp1 gene, and in the genus Fossombronia lost two introns of cox1, being the second and third largest mitogenomes of Pellidae (179,676 bp and 175,222 bp respectively). In contrast, over 70 kbp smaller mitogenomes of Apopellia contain complete introns set. The reduction of Apopellia mitogenome in comparison to the closest related Pellia is due to the loss of nucleotides in intergenic regions (Fig. 4). The intergenic spacers like atp4-trnMf, trnV-trnD, sdh4-nad4L, nad4L-tatC were over 80% shorter, while in the most of igs reduction fall within 40-60% range. Among the larger igs the smallest reduction, by 28.4% were in the atp1-cox1. Seven very short intergenic spacers were longer in Apopellia than in Pellia and in the one longer between cob-nad9 genes with the difference lower than 5%. Despite extraordinary reduction of igs, the identity between genes remains high, mostly above 90% (Fig. 4).  www.nature.com/scientificreports/ genera of Pallaviciniales. The order Pelliales further splits into two clades corresponding to genera Apopellia and Pellia. Within the former, all cryptic lineages including typical, water and spring forms were resolved as monophyletic clades with maximum bootstrap support. The water habitat lineages (water and spring) form a common clade sister to the typical form. The second clade of Pelliales were formed by species of the genus Pellia, revealing P. neesiana as early divergent within this group. Cryptic species of P. epiphylla formed monophyletic, maximally supported groups with a common clade. Allopolyploid species P. borealis was included in P. epiphylla N clade. The concentrated mitochondrial gene dataset consists of 45,852 nucleotides including 18,257 variable and 13,051 parsimony informative. The mitochondrial dataset generated a tree with similar, but not identical topology as the plastid dataset. Analysis of mitochondrial genes didn't support monophyly of Pelliidae and resolved Pelliales as earliest divergent among Jungermanniopsida, however support for common Fossombronialles/Pallaviciniales-Jungermanniideae clade wasn't maximal (BS 91%), while the plastid dataset supported monophyly  www.nature.com/scientificreports/ of Pellidae with maximum value. Another significant difference was revealed within suborder Lophocoleineae where species of Herbertus were placed as the earliest divergent. (Fig. 5). Minor differences comprise relationships within complex thalloid liverworts: (1) the mitochondrial dataset placed Asterella and Dumortiera within common clade (but poorly supported-BS 69%); (2) relationships at intraspecific level within Pellia neesiana and P. epiphylla N. The supernetworks ( Supplementary Fig. 1) based on single-gene trees of the mitochondrial and plastid dataset didn't reveal significant phylogenetic conflicts among genes, especially during diversification of the major liverwort lineages (Fig. 6). The competing splits within Pelliidae appeared only in the Pellia clade, where some minor incongruences in the mitochondrial dataset were detected among Pellia epiphylla N and P. neesiana lineages.

Mitochondrial genome variation.
Codon usage differentiation of Apopellia and Pellia. Comparative codon usage analysis revealed a difference between Apopellia and Pellia plastome genes, despite sharing by these genera identical sets of plastid trnS (Fig. 7). However, in the case of mitochondrial protein-coding genes, the CU patterns of both genera were similar and the loss of one trnR(ugc) gene by Apopellia mitogenomes didn't impact CU. The differences in codon usage among Apopellia and Pellia aren't equally distributed among amino acids ( Supplementary Fig. 2, Fig. 7). In the Apopellia plastomes most common stop codon 'taa' is preferred over 'tag' , while in Pellia all stop codons have similar frequency. Patterns of CU for amino acids like tyrosine, glutamine, glutamate, asparagine, aspartate and www.nature.com/scientificreports/ cysteine were similar in both genera. Where amino acids were coded by more than two codons, the differences in CU usually comprised two of them, like in the case of alanine, where 'gct' codon is more frequent than 'gcc' in Apopellia, while the opposite situation was revealed in Pellia.
The general pattern of CU in the mitogenome and plastome was similar (Fig. 7) and as mentioned above, without clear differences between studied genera. Within Apopellia differences were found in 'water' form, in which serotonine 'tcg' and 'tct' codon frequencies are different than in remaining species. This kind of variation   Table 2). The lowest number of MDCs was revealed for P. borealis, while the highest for P. neesiana. The cryptic species of the P. epiphylla complex have 2979 and 3307 MDCs for N and S respectively. Among Apopellia endiviifolia cryptic species the biggest number of MDCs was identified for the typical form (20,260), while the water and spring form revealed 309 and 264 MDCs respectively.
Comparison of infra-and interspecific distances based on the plastome confirmed the presence of a significant barcoding gap among all species, including cryptic ones (Supplementary Fig. 3). In the case of mitochondrial genomes, the barcoding gaps are smaller but still significant. The share of nucleotide diagnostic characters is not equally distributed among analysed species (Fig. 2). ASAP analysis performed for data based on the chloroplast genome generated four identical partitions with low values of asap-score ranging from two to three (Fig. 8a). These four proposed divisions indicated six species among testing liverwort groups. These results fully reflect the species split among the cryptic complexes in the genera Pellia and Apopellia. Within the Pellia genus, cladogram indicated that P. epiphylla N and P. epiphylla S are definitely separated species (black circle, Fig. 8a). As different species, cladogram also clearly showed two forms of Apopellia endiviifolia: water and spring (black circle, Fig. 8a).
The mitogenome-based dataset gave slightly weaker results. Only one partition (asap-score = 7.5) proposed by ASAP was congruent with real division of Pellia and Apopellia, indicating a total of six species within these genera (Fig. 8b). On the other hand, one partition with quite low asap-score with a value of 3.5 grouped together only water and spring form of Apopellia. Although two ASAP-partitions put forms of Pellia epiphylla together,  www.nature.com/scientificreports/ the cladogram next to the bar chart suggests that Pellia epiphylla is not a homogeneous group and the division into N and S types is justified (red circle, Fig. 8b).

Discussion
Organellar genomes of Pellidae. The liverwort organellar genomes are known from its conserved structure and gene content and results of Pelliidae analysis confirm this observation. In contrast to the chloroplast genomes of angiosperms, the structural heteroplasmy connected with SSC subunit orientation wasn't detected in Pellidae, as well as in previous studies on leafy liverwort Scapania undulata 64 . Long read sequencing identified only one plastome haplotype also in gymnosperms and pteridophytes, suggesting that alternative haplotype is specific for angiosperms 65 . The gene order of Pelliales and Fossombroniales is typical for most of liverworts, but plastome structure of Moerckia confirms that the translocation of rps12-rps7-ndhB cluster seems to be characteristic for Pallaviciniales 19 .
Plastomes of the Pellia genus lost two genes (cysA and cysT), but gene order remains stable. These genes are present in the sister genera Apopellia, Makinoa and Moerckia 19,28 , but in the remaining genera of Pelliidae (Fossombronia, Pallavicinia, Pellia) are pseudogenised, as well as in many other leafy liverworts 19 . Further gene losses include trnS(GCU)-psbI-trnS(GCA) in Makinoa and Fossombronia. However, plastome length in Pellidae doesn't correlate with gene losses, while the second largest plastome of Pallavicinia (124,103 bp) and the smallest of Pellia neesiana (115,875 bp) have identical gene content, whereas plastomes of Apopellia with complete gene sets (120,537-120,947 bp), are smaller than Pallavicinia without cysA and cysT genes. Independently of gene loss, plastome reduction was also found in leafy Cephaloziineae, which plastomes ranged between 118,571 and 114,423 bp, despite identical gene sets 28,64 . Plastid gene losses are scattered along different evolutionary lineages www.nature.com/scientificreports/ including loss of cysA and cysT genes in Ptilidium pulcherrimum and Cheilolejeunea xanthocarpa 33,35 and additional four genes (ndhF, rpl21, rps32 and ccsA) in Cololejeunea lanciloba 66 . Single genes are also missing in Aneura pinguis (ycf66), Haplomitrium blumei (psaM) and Schistochila macrodonta (ndhF) 66,67 . Heterotrophic Aneura mirabilis has a reduced plastome due to loss of six genes (ndhA-E-G-K-I and ycf66) or pseudogenization (19 genes) of genes involved in the process of photosynthesis. Two liverwort species increase the number of plastid genes by transferring them to IRs: Haplomitrium blumei transferred ndhF gene from SSC 66 , while Conocephalum salebrosum transferred rps7 and rps12 from LSC 68 . Minor changes in genomes are related to loss of introns: complex thalloid liverworts, with the exception of Blasiales, lost introns of ycf3 and rps12 genes, which phenomenon was also reported in species from distantly related genera like Herbertus, Ptychantus, Radula and Metzgeria 19,35,68 .
Single loss of introns of ndhA and ndhB genes were found in Radula japonica and Treubia lacunosa 19 .
In contrast to plastomes, evolutionary losses of mitochondrial introns seem to be common in liverworts, but mainly concern leafy liverwort lineages, while the intron content of thallus liverworts is stable 29,30,69 . Out of ten intron losses only the loss of trnL could happen once and remaining losses were independent across different evolutionary lineages of leafy liverworts 70 . Unlike previously sequenced genomes of Pelliidae (Fossombronia cristula and Makinoa crispata), the mitogenomes of Pellia and Apopellia have introns in the atp1 gene and in contrast to Fossombronia, keep complete cox1 intron set 30,70 . The gene content of liverwort mitogenomes is more stable than that of the plastomes, beside the nad7 loss, which is present only in early divergent lineages of Haplomitriales and Treubiales and the loss of ccmC and ccmFn in Treubia lacunosa 71 . The variation of tRNA genes among liverworts comprises the loss of trnT-GGU after divergence Marchantiopsida and Jungermanniopsida and trnR-UGC after divergence of Pellidae. The trnR-UGC seems to be independently lost in Apopellia, even considering incongruent results on monophyly or paraphyly of Pelliae 14,66,70,72 , which doesn't question monophyly of Apopellia/Pellia clade 14,26 . The loss of the trnR-UGC gene didn't impact codon usage pattern of mitochondrial protein coding genes-both genera revealed similar CU pattern (Figs. 7, 8). On the other hand, Apopellia and Pellia differ in CU by plastome CDS, despite identical tRNA gene content. The lack of correlation between codon usage and the www.nature.com/scientificreports/ number of tRNA was previously reported in bacterial genomes 73,74 . Procaryotes seem to respond to gains or losses of tRNA genes by modulating the concentration of tRNAs rather than modifying its frequency on codon usage 73 .
Despite the presence of some repeats, which are considered as a main factor of mitochondrial structure dynamics in vascular plants 75,76 , the structure of liverwort mitogenomes is stable 37 . The exceptions were found for the mitogenome of Gymnomitrion concinnatum belonging to Jungermanniales, where two recombination events were detected 34 and in Dumortiera hirsuta (Marchantiales) where one 96 kb inversion was found 77 . The possibility of recombination was also suggested for 12 other liverworts, however most frequent variants support typical mitogenome structure 29 .
Despite the nearly identical gene content (loss of trnR-UCG in Apopellia) and order the mitogenomes of Pellia and Apopellia differ in size. The newly sequenced Apopellia mitogenomes were smallest among known liverworts, over 49 kbp smaller than those found in Pellia and falling within 109,022-109,828 bp in range. The size of Apopellia mitogenomes is closer to the values found in mosses 37,78,79 , than to those known in liverwort mitogenomes ranging from 142,510 bp in Tritomaria quinquedentata to 187,628 bp in Monosolenium tenerum 29,69 . The size reduction of Apopellia mitogenome is not caused by gene or intron losses, but by the reduction of intergenic spacers, which affects all regions with exception of only few spacers like cob-nad9, cox3-nad9 and usually long atp1-cox1, which was reduced only by 31%. The reduction or increase of mitogenome length in leafy liverworts is often caused by indels between genes nad3 and rpl10. This region contains the pseudogenized (with exception of Haplomitrium and Treubia where nad7 is functional) nad7 gene and three tRNAs: trnV-UAC, trnD-GUC, trnS-GCU and is a main source of mitogenome size variation among Nowellia and Scapania species 28 . Comparative analyses of Pelliidae mitogenomes also confirmed this observation revealing over 11 kbp difference in nad3-rpl10 region between Makinoa and Fossombronia species, with differ by 8 kbp in total mitogenome length (Makinoa contains two more cox1 introns). This region is also twofold smaller in Apopellia (5.3 kbp) than in Pellia (11.4 kbp), but it comprises only 11% of total mitogenome reduction in the former.
Phylogenetic placement of Pelliidae based on plastid sequence is mostly congruent with previous studies 28,35,70 . One minor difference is the position of Conocephalum conicum, which in the work of Dong et al. was resolved as an unusual clade, probably due to plastome misassembly (GB accession number MK645816). The latest study on complex thalloid phylogeny confirms the position of Conocephalum as sister to Riccia 80 . Chloroplast dataset supports monophyly of Pelliidae, however the mitogenome dataset resolved Fossombroniales and Pallaviciniales in common clade with remaining Jungermanniopsida, while Pelliales was resolved as the earliest divergent from that class. Cyto-nuclear incongruences are common across phylogeny of vascular plants 81 and were also confirmed within Pellidae, where the application of transcriptomic data revealed that mito-phylogenetic results support nuclear gene topology 19 . Incongruences among nuclear, mitochondrial and plastid dataset are often explained by incomplete lineage sorting 82 , horizontal gene transfers 83 and hybridization events 84 , however in the case of Pellidae aforementioned mechanisms are unlikely. The organellar genomes of Pellidae seem to evolve much faster than any other subclass of liverworts.
The topology of the obtained trees is congruent with previous studies on Pelliales and fully supports splitting of Pellia sensu lato by placing Pellia endiviifolia in the separate genus Apopellia 26 . The splitting of the Pellia genus is not only supported by typical sequence phylogenetic divergence but also by differences in the structures and genome gene content as discussed above. Moreover, the mitochondrial dataset revealed a paraphyly of Pelliaceae and resolved Apopellia as a sister group to Fossombroniaceae. Super-barcoding and molecular diagnostic characters. Application of organellar genomes as superbarcodes enabled identification of all analysed evolutionary lineages of Pellia s.l. with the exception of allopolyploid P. borealis which inherited organellar genomes from P. epiphylla N 85 .
The number of plastid MDCs characterising Apopellia species decreased in comparison to previous study 28 , due to expanding the dataset by a spring form, which reveals some intermediate characters between water and terrestrial forms, but it still, like P. epiphylla S and N, falls within the range known from Calypogeia-a taxonomically good liverwort species 86,87 .
Based on our data it is hard to estimate the effect of multiple specimen sampling on presence of diagnostic nucleotides or MDCs. In four out of six species (including cryptic) the impact of sample number within 2-6 range on plastid MDCs was below 1% (species of Apopellia and P. neesiana). The decrease in plastid MDCs by 31% was detected in P. epiphylla S (N = 3) and by 15% in P. epiphylla N (N = 6), but it was not correlated with sample size. The mitogenomes reveal contrasting patterns, but again, not correlated with sample size. The decrease in MDCs in the case of P. epiphylla cryptic species was 7% and 9% for S and N respectively, but in the case of A. endiviifolia water form decrease in MDCs was 35%.
The mitogenomes deliver more MDCs than the plastome, however two factors should be discussed. The number of detected MDCs depends on alignment size, which in the case of mitogenomes was as much as 71% longer, however it doesn't explain the over 5.6-fold higher number of MDCs. In this particular case, the main reason for the outnumber of MDCs in mitogenome is a bigger number of large indels in mitochondrial dataset (1385 indels with mean length of 146.2 bp) in comparison to plastid dataset (1629 indels of 9.1 bp mean length). Thus the mitogenomes were more efficient in delimitation of cryptic species within the P. epiphylla complex, while plastomes performed better in the case of A. endiviifolia complex. These results can be explained by differences in genome sizes, since the mitogenome of Apopellia is smaller than plastome, which is quite unusual in plants, with exception of mosses 37,78,88 .
Comparative organellar genomics support previous hypotheses about allopolyploid origin of Pellia borealis, which pointed out Pellia epiphylla N as a maternal species 85 . Both organellar genomes found in P. borealis were almost identical to those found in P. epiphylla N. The lack of plastome MDCs and only four MDCs in mitogenomes suggest the late Pleistocene origin of P. borealis, when ranges of the N and S lineages could form a natural www.nature.com/scientificreports/ hybridization zone. The allopolyploid Calypogeia muelleriana, which carries the maternally inherited plastome of C. integristupula, revealed specific SNPs enabling molecular delimitation of this species based on the plastome and mitochondrial datasets 86,87 . The distribution of Pellia epiphylla cryptic species can't be described as northern (N) and southern (S) lineages and their ranges should be revisited. Previous studies revealed geographical vicariance despite similar ecological requirements 22,89 , but nowadays we can't confirm this pattern, since N lineage was also found in Carpathians. Coexistence of cryptic species lineages is common in other known complexes as Aneura pinguis 3 or Apopellia endiviifolia 22 .
Our results confirm the presence of the third cryptic lineage named "C" within Apopellia endiviifolia, described previously by Polok et al. 24 as "well-headed" form. Phylogenetic analyses resolved species C as a sister to A which corresponds to similar ecological preferences of these forms. However, based on current knowledge on distribution of species C, its habitat preferences are much narrower than in the case of species A. All known specimens of species C occurred directly in the springs, they cover the bottoms of forming watercourses. The earlier studies suggested that species A is known mainly from the mountains in the direct neighbourhood of running water, with very restricted occurrence in the lowlands, which corresponds to preferred microhabitat distribution 90 . The sites with clear and fast running shallow water, which this lineage depends on, are more common in the mountains. The geographic distribution of species B is concentrated in lowlands with rare occurrence in the mountains and the species prefers denuded soil and rock detritus 90 . The species C is currently known from five sites, including three lowland and two mountain locations 24 and its scattered distribution can be explained by presence of spring areas. Morphological studies on A. endiviifolia revealed significant differences between cryptic species in quantitative traits like size of basal and marginal cells, even after six months of greenhouse culture 24 . The species C had the smallest basal and marginal cells within A. endiviifolia complex. However, further ecological and morphological studies are needed to better characterise cryptic species within A. endiviifolia which could lead to formal taxonomic recognition. A species called Apopellia megaspora was identified in the early 1980s as one of the cryptic North American species of A. endiviifolia 91 . In addition to differences in spore diameters, it is distinguished from other Apopellia species by the absence of thallus proliferation in autumn. Another cryptic species of simple thalloid liverwort, Aneura pinguis, has been more thoroughly studied, but despite the presence of various lineage-specific DNA markers 3,67,92 and many biochemical compounds 93 , it lacks significant morpho-anatomical differentiation. Furthermore, the cryptic lineages of A. endiviifolia and some of A. pinguis differ in their preferred microhabitats, which can be characterised not only by general ecological descriptions, but also by specific parameters such as pH and concentrations of calcium, magnesium, potassium, and sodium 94 .
Results presented in this study provide novel insights on evolution of liverwort organellar genomes, revealing high differentiation of mitogenomes among closely related taxa. Both plastid and mitochondrial genomes proved their usefulness in molecular delimitation of cryptic species of Apopellia and Pellia. Surprisingly, effectiveness of mitogenomes in species discrimination within Pellidae was on a similar level as plastome.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.