The specimens of O. stentina used for this study were sampled in the Venice Lagoon (Italy), one of the widest coastal transitional ecosystems in the Mediterranean, and a complex mosaic of habitats. In particular, they were collected in the intertidal zone at the following locations: O. stentina NEWOS06 (Torson di sotto island: 45°20’55.2” N, 12°13’46.4” E); O. stentina NEWOS18 (Darsena dell’Arsenale: 45°26’14.5” N, 12°21’18.1” E); O. stentina NEWOS81 (Faro Rocchetta: 45°20’20.7” N, 12°18’39.6” E). Muscle tissue was stored in 100% ethanol and kept at -20°C until DNA extraction.

Genomic DNA was extracted using the commercial kit Invisorb Spin TissueMini Kit (Invitek, STRATEC Biomedical, 242 Germany), quantified using Qubit dsDNA BR Assay Kit (Invitrogen–ThermoFisher Scientific), and checked for quality on agarose gel electrophoresis.

Genomic libraries were constructed using the commercial kit Illumina DNA Prep (Illumina, Inc.), quantified using Qubit dsDNA HS (High Sensitivity) Assay Kit (Invitrogen–ThermoFisher Scientific), and checked for quality on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California, USA), before sequencing. Libraries were equally pooled and sequenced on an Illumina HiSeq4000 platform with 150 bp pair-end read module at the UCDavis DNA Technologies & Expression Analysis Core (Davis, CA) in order to obtain 14 M of raw read-pairs/library.

Raw paired-reads obtained from Illumina sequencing were assessed for quality using FastQC v0.12 (Andrews, 2010), and consequently trimmed of any adaptors and low quality sequences using Trimmomatic v0.32 (Bolger et al., 2014); high quality reads of length ranging from 70 bp to 150 bp were retained. The whole mitogenome was assembled using the software Get Organelle version 1.7.7.1 with kmer values of 21, 45, 65, 85, 105 and “animal_mt” as seed database (Jin et al., 2020).

The new mitogenomes of O. stentina were annotated following the strategy described in previous works from our laboratory (Babbucci et al., 2014; Montelli et al., 2016; Basso et al., 2017). Gene nomenclature followed standard naming for animal mitogenomes (Montelli et al., 2016; Basso et al., 2017; Ghiselli et al., 2021). A comparison was done with an ESTs library of Magallana gigas (Crassostrea gigas digestive gland subtracted library, multiple accession numbers, unpublished) available in GenBank (Clark et al., 2015) to ensure a correct identification of the 5’ start of the protein coding genes that are particularly difficult to be identified in molluscs (Ghiselli et al., 2021). For the identification of the two halves of rrnL we followed Milbury et al. (2010). Finally, comparisons with other published and unpublished oyster mitogenomes were done to refine the annotations (Milbury and Gaffney, 2005; Ren et al., 2009, 2010; Wu et al., 2010; Danic-Tchaleu et al., 2011; Yu and Li, 2011; Wu et al., 2012; Xiao et al., 2015; Cavaleiro et al., 2016; Ren et al., 2016). The mitogenomes used in the analyses, and listed in Table 1 were re-annotated following the approach described above to ensure a consistent annotation across all the sequences. Different authors have assigned various names to the mitogenome strands (Basso et al., 2017). In the present paper, we refer to the strand encoding all the genes as the plus strand and the opposite strand as the minus strand. Because our study was focused on the subfamily Ostreinae we analyzed all available mitogenomes for this taxon, while for the subfamilies Crassostreinae and Saccostreinae, we restricted our analyses to selected species ( Table 1 ). However, two mitogenomes of Ostreinae (Ostrea denselamellosa ON964460 and O. edulis CM063324) became available in GenBank (Clark et al., 2015) after the main phase of our analyses had been completed. Consequently, they were not fully integrated into the study but were included in selected downstream analyses (see below).

Conserved gene blocks shared across two or more gene orders were identified through visual inspection of the gene arrangements of CraGO, MagGO, OstGO, SacGO, and TalGO. The analysis focused on gene blocks shared among taxa from different subfamilies. The ancestral GO for the Ostreidae was inferred manually, according to a principle of parsimony, by mapping GOs evolution along the reference phylogeny of the family (Salvi and Mariottini, 2017; Li et al., 2021; Salvi and Mariottini, 2021).

Multiple alignments of the protein-coding genes (PCGs) were done in two steps. Initially, for each PCG an alignment of the amino acid sequences was performed with the online version of the MAFFT program (https://mafft.cbrc.jp/alignment/software/) (Katoh et al., 2002). Successively, through the TranslatorX server (http://161.111.160.230/index_v5.html), the codons of each orthologous set of PCGs were aligned using as template the corresponding amino acid multiple alignment (Abascal et al., 2010). Multiple alignments of the orthologous tRNAs were produced in two steps. Firstly, orthologous sequences were quickly aligned with ClustalW (Thompson et al., 1994). Successively, these alignments were improved manually (Montelli et al., 2016), considering the secondary structures predicted for each tRNA with the tRNA-scan software (Chan and Lowe, 2019). To provide a figure of the conservation level of each tRNA multiple alignment, a logo vignette was created with the WebLogo software (Crooks et al., 2004). The multiple alignments of rRNAs were produced with the MAFFT program (https://mafft.cbrc.jp/alignment/software/) (Katoh et al., 2002). The secondary structure of rrnL of O. stentina was manually modeled using the published structure of Magallana gigas (Milbury et al., 2010) as template.

The occurrence of the intergenic spacers was studied in all the mitogenomes of Ostreinae and their distribution was mapped on the reference mitogenome (see below). Alignments were produced for each group of ISPs. The shortest ISPs were aligned manually, while the longest were aligned with the online version of the MAFFT software (Katoh et al., 2002). ISPs are very fast evolving sequences, therefore in most of the cases the alignments were restricted to a single species or to closely related species (see below). We searched for the occurrence of one or more ISPs representing molecular signatures for the analyzed oysters. To qualify as a molecular signature, an ISP must exhibit a unique sequence characterizing the mitogenomes of one species or multiple taxa forming a monophyletic group in a phylogenetic tree.

In the mitogenomes of O. stentina the Control Region (CoRe) was identified as the longest intergenic spacer containing AT-rich motifs, stem-and-loop, and cloverleaf secondary structures (Ghiselli et al., 2021). The capability to produce secondary structures was tested with the software RNAstructure (Reuter and Mathews, 2010). Multiple alignment of the CoRes of Ostreinae was done with the T-Coffee. The web server version, using the M-Coffee option, aligns DNA sequences by combining the output of popular aligners (Di Tommaso et al., 2011).

The total number of codons and the relative abundance of each codon family used by the 13 PCGs were computed with the MEGA X program for all Ostreinae (Kumar et al., 2018). The codon distribution was expressed as number of codons per thousand codons (CDSpT). The Relative Synonymous Codon Usage (RSCU) values were calculated with the MEGA X program (Kumar et al., 2018). First codons, as well as stop codons, complete and incomplete, were excluded from the analysis to avoid biases due to unusual putative start codons and incomplete stop codons.

The A+T/G+C content and GC-skew = (G-C)/(G+C) and AT-skew = (A-T)/(A+T) (Perna and Kocher, 1995) were used to measure the compositional biases among analyzed sequences. The base compositions were computed with MEGA X (Kumar et al., 2018). The calculations of skews were performed with Excel program (Microsoft TM). The skews were computed for the whole mitogenomes, and for PCGs, tRNAs, rRNAs and CoRes. Differences in AT-skew vs A+T content and GC-skew vs the G+C content were plotted as scatterplots in Microsoft Excel.

We also tested whether the AT- and GC-skews of various PCGs, tRNAs, rRNAs and CoRes were statistically significantly different from those computed for the whole mitogenomes. Levene’s test (Levene, 1960) revealed that the variances of AT- and GC-skews differed significantly across gene/regulatory regions (PCG, tRNA, rRNA, or CoRes), violating the assumption of homogeneity of variance required for linear models (AT-skew: F_{(23, 240)} = 2.550, p < 0.001; GC-skew: F_{(23, 240)} = 3.169, p < 0.001). For this reason, differences in AT-skew and GC-skew between the whole mitogenomes and the genes/regulatory regions were tested by bootstrapping. For each combination of mitogenome and gene/regulatory region, 50,000 bootstrapped data sets were generated. For each of these data sets, we calculated and stored the difference in AT-skew and GC-skew between the mitogenome and gene/regulatory region. To calculate the p-value of the difference in AT-skew and GC-skew between the mitogenome and genes/regulatory regions, we generated a distribution of the difference under the null hypothesis as follows: 1) for each combination of mitogenome and gene/regulatory region, we first calculated the difference in AT-skew and GC-skew in the original data set; 2) this value was added to the AT-skew or GC-skew of the gene/regulatory region, forcing the difference in AT-skew or GC-skew between the mitogenome and gene/regulatory region to be equal to zero; 3) from this new set of data, we created 50,000 bootstrapped data sets and, for each of these, calculated the difference in AT-skew or GC-skew between the mitogenome and gene/regulatory region, allowing us to obtain the distribution of the difference in AT-skew or GC-skew under the null-hypothesis. The p-value of the difference in AT-skew and GC-skew between the mitogenome and gene/regulatory region was calculated as the probability of obtaining a result equal to or more extreme than what was observed in the first bootstrapping, assuming the null hypothesis (no difference) was true.

The nucleotide substitution pattern was tracked in the stems of the secondary structures of orthologous tRNAs (Montelli et al., 2016). We looked for the occurrence of (a) hemi-conservatory base changes, (b) type I fully compensatory bases changes, (c) type II fully compensatory base changes, and (d) mismatches (Montelli et al., 2016). These patterns were identified by visual inspection of multiple sequence alignments and taking into account the predicted secondary structure of tRNAs (Montelli et al., 2016).

Given the same pair of bases in a stem, a change with respect to the background condition for the multiple alignment can involve only one of the two bases, either at the 5′ or 3′ end, without altering the secondary structure of the stem (e.g., T•G vs. T–A) (Coleman, 2003; Montelli et al., 2016). This variation is referred to as a hemi-compensatory base change (Coleman, 2003; Montelli et al., 2016). The change can also involve both bases of the pair but the secondary structure remains intact (e.g., G–C vs. A–T). This type of change is known as a fully compensatory base change because the substitution of both bases does not compromise the integrity of the secondary structure (Coleman, 2003; Montelli et al., 2016). There are two types of fully compensatory base changes: type I, which involves the substitution of a purine-pyrimidine pair with another purine-pyrimidine couple and vice versa (Montelli et al., 2016); and type II, which is characterized by a purine-pyrimidine vs. pyrimidine-purine substitution (Montelli et al., 2016). Type I occurs more easily than type II because its intermediate step is represented by a hemi-compensatory base change (Montelli et al., 2016). In contrast, type II is disfavored because its intermediate step involves a mismatch in the pair that jeopardizes the secondary structure of the stem (Montelli et al., 2016). Lastly, the change can involve a substitution pattern leading to a disruption of the secondary structure of the stem for the analyzed pair. Mismatches that do not prevent the formation of the cloverleaf structure or the tertiary structure are not uncommon in tRNAs. Various mechanisms of editing can correct mismatches in the stems, or alternatively, these mismatches can persist as unusual pairings (Cannone et al., 2002).

The level of homogeneity/heterogeneity in the substitution process in PCGs and their corresponding protein products, as well as in the two rRNAs multiple alignments, was tested with the AliGROOVE software (Kück et al., 2014). For the PCGs the AliGROOVE matrices were computed for: complete codons, the first plus the second position of each codon, single positions (first, second and third), and the translated amino acid sequences. AliGROOVE tests were performed also on the multigene concatenated data sets used in the final phylogenetic analyses (see below). In a matrix, obtained from AliGROOVE, a square ranging from brown to pink identifies a heterogeneous substitution process between the two compared sequences, while a square ranging from light to dark blue marks a homogenous substitution process (Kück et al., 2014).

The phylogenetic signal present in the different genes/multiple alignments was evaluated through two different methods: (a) the quartet puzzling analysis (Strimmer and von Haeseler, 1996) and (b) the boxplot graphics, which analyses the distribution of the pairwise distances computed according to the best-fit evolutionary model (e.g. Negrisolo et al., 2004). The quartet puzzling analysis was performed as implemented in IQ-TREE2 (Minh et al., 2020). The best fitting evolutionary models for DNA/proteins were identified with the ModelFinder program (Kalyaanamoorthy et al., 2017) implemented in IQ-TREE2 (Minh et al., 2020). IQ-TREE2 software was used also to compute the distances based on the best-fit models. The boxplots were created with the Excel software. The occurrence of phylogenetic signal was studied on the rrnS and rrnL and on each PCG. In this latter case, the analysis was done on the single positions of the codons (p1, first; p2, second; p3, third), on positions one and two (p12), overall codons (p123) and on the translated polypeptides. Finally, this analysis was extended to the concatenated data sets (see below).

For phylogenetic analyses, we created 10 concatenated data sets that are listed in Table 2 . The concatenations were done with the MEGA X software (Kumar et al., 2018).

The phylogenetic analyses were performed according to the maximum likelihood (ML) method (Felsenstein, 2004). The ML trees were computed with the program IQ-TREE2 (Minh et al., 2020). In the tree search analysis, 50 independent runs were performed for the rRNA data set, 20 runs for 13PCGpro, 13PCGpro+rRNA, 13PCGnuc, 13PCGnuc+rRNA and 10 runs for 13PCGp1, 13PCGp2, 13PCGp3, 13PCGp12 and 13PCGp12+rRNA. The optimal partitioning scheme (Chernomor et al., 2016) and the best fitting evolutionary models (Kalyaanamoorthy et al., 2017) were selected with IQ-TREE2. The data sets 13PCGnuc, 13PCGnuc+rRNA and 13PCGp12+rRNA were analyzed also with a partition scheme that considered the single positions of the codons separately.

The Ultrafast Bootstrap Test (UBT) (Minh et al., 2013) (10,000 replicates) and the approximate Likelihood-Ratio Test for branches (aLRT) (Anisimova and Gascuel, 2006) (1,000 replicates) were used for evaluating the robustness of the tree topologies obtained in the various searches. A Robinson-Foulds distance data matrix (Llabrés et al., 2021) was computed for each set of the trees generated in every phylogenetic analysis to ensure that the top-ranked topologies had a null distance and the convergence had been reached in the tree searching.

To evaluate alternative phylogenetic hypotheses, topology tests were done according to the Approximately Unbiased (AU) test (Shimodaira, 2002), the Weighted Shimodaira and Hasegawa (WSH) test (Shimodaira and Hasegawa, 1999) and the Expected Likelihood Weights (ELW) (Strimmer and Rambaut, 2002) method. Computations were performed with IQ-TREE2 (Minh et al., 2020).

Three complete mitogenomes of O. stentina were sequenced for this work ( Figure 1 ; Supplementary Table S1 ). The Illumina reads used to assemble these genomes spanned from a minimum of 12,236,243 to a maximum of 17,149,758 ( Supplementary Table S1 ). The size of the O. stentina mitogenomes varied from 16,305 to 16,314 bp ( Supplementary Table S1 ). This range was very similar to values obtained for the mitogenomes of Ostreinae and Saccostreinae, while values were higher in Crassostreinae also for the occurrence of extra genes (Magallana and Talonostrea) ( Figure 2A ; Table 1 ). The mitogenome of O. stentina contained a set of 38 genes: 13 PCGs, 23 tRNAs and 2 ribosomal RNAs ( Figure 1 ; Supplementary Table S1 ). The gene order corresponded to the typical arrangement observed in the Ostreinae subfamily (OstGO), with all genes encoded on the plus strand ( Figure 1 ). Genes were contiguous or separated by intergenic spacers (ISP) ( Figure 1 ; Supplementary Table S1 ).

The mitogenomes of Ostreinae sequenced so far exhibited the same GO ( Figure 2 ; Table 1 ). However, they varied at the microstructural level in the distribution of the ISPs (see below). OstGO exhibited the maximum level of synteny with SacGO of Saccostreinae as proved by the sharing of two large conserved gene blocks ( Figure 2A ). OstGO shared gene blocks also with CraGO, MaGO and TalGO ( Figure 2A ). Furthermore, blocks containing two or more genes were shared among all analyzed GOs ( Figure 2A ). Thus, by considering the distribution of the conserved blocks among different GOs and the reference phylogeny for the Ostreidae, we identified gene blocks that were present in the gene order of the last common ancestor of the Ostreidae (lcaO, Figure 2B ).

The mitogenomes of O. stentina were A+T rich, negatively AT-skewed, and positively GC-skewed ( Figure 3 ). This feature was shared by all oysters sequenced to date, limiting our comparison to Ostreinae ( Table 1 ). Among Ostreinae, Ostrea denselamellosa presented the most diverging values for both A+T/G+C content (60.71%) and AT-skew (-0.153) ( Figure 3 ). The range of variation of AT-skew was broader than that of GC-skew. Mitogenomes of Ostreinae exhibited a limited range of variation in A+T (60.71% - 65.41%)/G+C (34.58% - 39.29%) content and AT-/GC–skews (-0.153 − -0.128) (0.149 − 0.201) ( Figure 3 ). However, the taxon coverage was very low.

O. stentina exhibited the whole set of PCGs usually present in animal mitogenome ( Figures 1 , 2 ). All the PCGs started with standard codons (ATT, ATG, GTG, and TTG) and ended with the canonical TAA codon, except for cox3, atp6, nad5 and nad3, which ended with incomplete stop codons T(aa) or TA(a) ( Supplementary Table S1 ). None of the PCGs overlapped. When the comparison was extended to all Ostreinae sequenced to date ( Supplementary Table S2 ) it appeared that ATG was the most widespread codon followed by GTG. The genes using the most variable sets of start codons were cob (5) and nad2 (5), while the pair cox1-cox2 invariably started with ATG.

The codon distribution and the Relative Synonymous Codon Usage (RSCU) in PCGs were analyzed for the different mitogenomes of oysters of the subfamily Ostreinae. The results are summarized in Figure 4 and Supplementary Figures S16, S17 . The average number of codons for the subfamily Ostreinae was 3,710. The range of variation among the taxa was limited, spanning from 3,699 in D. sandvichensis to 3,717 codons in O. lurida and O. edulis. No intraspecific variation was observed for the multiple mitogenomes of O. stentina ( Supplementary Figure S16A ), while a very limited variability was detected in O. edulis ( Supplementary Figure S16B ). All Ostreinae exhibited a very consistent codon distribution and RSCU ( Figure 4 ; Supplementary Figures S16, S17 ). The most abundant amino acids in mitochondrial proteins determined also the richest codon families. Ser, Leu, Val and Phe were present with more than 300 residues in all taxa, accounting for ≥ 42% of the whole set of amino acids across the 13 proteins. Gly, Ala, Ile, and Met occurred with more than 200 residues (25% of the whole set), while other amino acids were less abundant, with Gln being constantly the rarest. This distribution explained the pattern observed for the codon families in Figure 4 , and Supplementary Figure S17 . The Val codon family was the most abundant, as both Leu and Ser were split into two families, with Leu2 favored over Leu1, and Ser2 more represented than Ser1. The analysis of RSCU showed that the A+T rich codons were preferred over synonymous codons with a lower content in A+T ( Figure 4 ; Supplementary Figure S17 ). This result was expected considering the compositional bias toward A+T exhibited by all PCGs (see above). However, all codons were used at least once, as the compositional bias was not so extreme to determine the elimination of rare GC-rich codons. The combined effect of A+T richness, negative AT-skew and positive GC-skews of PCGs on the codon composition was well represented by the behavior of some fourfold-degenerated codon families where the abundance of the third codon base followed this order: T, A, G, C (e.g. Pro, Val). However, this pattern was not always consistent (e.g. Ser2) suggesting that other factors played a role in the final abundance of synonymous codons in the mitogenomes of oysters. The p-distances were computed for orthologous genes, including single codon positions, and proteins of Ostreinae ( Supplementary Table S3 ). The most variable gene/protein was atp8, followed by nad2 and nad6 while the most conserved was cox1. At the intraspecific level, the variability was very low among the two genomes of O. denselamellosa ( Supplementary Table S4 ), the four mitogenomes of O. edulis ( Supplementary Table S5 ), and the three sequences of O. stentina ( Supplementary Table S6 ).

The mitogenome of O. stentina contained the full 22 tRNAs set of Metazoa plus a duplicated trnM2, a feature shared among all oysters sequenced to date ( Figure 2 ). All tRNAs exhibited the clover-leaf secondary structure ( Figure 5 ; Supplementary Figures S18–S20 ). The analysis of the multiple alignments of orthologous tRNAs ( Figure 5 ; Supplementary Figures S18–S26 ) revealed different levels of conservation among Ostreinae. Most of the variable positions were located in the single helix portions of the tRNAs, i.e. DHU loop, “extra arm and TΨC loop, which were free to vary without hampering their structure ( Figure 5 ; Supplementary Figures S18–S26 ). Some of these hyper-variable portions characterized single taxa (e.g. TΨC loop of trnC for O. stentina and O. edulis; Supplementary Figure S21 ). Base substitutions in the stems were prevalently hemi-compensatory and type I fully compensatory base changes ( Figure 5 ), as they maintained the integrity of the stems, and the molecular pathways leading to them are favored (Montelli et al., 2016). Type II fully compensatory base changes, requiring intermediate mismatches, were much rare but occurred in trnA and trnF of O. denselamellosa ( Supplementary Figures S21, S22 ), and in trnH, trnL1, trnM1, trnN and trnW of all Ostreinae ( Supplementary Figures S21–S24, S26 ). Mismatches were also present in the stems and restricted to single species (e.g. O. edulis, trnA; O. stentina, trnE) ( Supplementary Figures S21, S22 ) or common to all Ostreinae (e.g. trnD, trnN, trnQ, trnR, trnV) ( Supplementary Figures S21, S24–S26 ).

The most conserved tRNAs was trnS1 with only 5 variable positions over 70, whereas almost 50% of positions (30/66) changed in trnC, the most dynamic tRNA in Ostreinae ( Figure 5 ; Supplementary Figures S21, S25 ). The tRNAs associated to the most abundant codon families were the least variable ( Figure 4 ; Supplementary Figures S18–S20 ). The notable exception was represented by trnQ, which was associated to the least abundant amino acid ( Figure 4 ) but was among the most conserved tRNAs ( Supplementary Figures S19, S24 ). At the intraspecific level, only one base (T vs C) difference was found in the DHU loop of trnR of O. stentina. However, the three mitogenomes analyzed here were obtained from specimens of the same locality ( Supplementary Figures S27–S30 ). In contrast, variable tRNAs were present in O. denselamellosa/O. edulis ( Supplementary Figures S27–S30 ).

The rrnSs were conserved among the analyzed Ostreinae (average p-distance = 0.167 ± 0.077). The G+C content was higher than both the 5’ and 3’ halves of rrnL ( Supplementary Figure S14 ), which suggested a strong and important role of the G-C pairs in the 2D/3D structures of this molecule. The intraspecific variability was minimal for the three species with multiple sequences available (average p-distance = 0.002 ± 0.002 in O. stentina; average p-distance = 0.006 ± 0.007 in O. edulis; average p-distance = 0.001 in O. denselamellosa).

As secondary structure models existed for rrnL of the phylum Mollusca (Lydeard et al., 2000) and for the family Ostreidae (Milbury et al., 2010), we used these templates to infer the secondary structure of the rrnLs of Ostreinae.

The overall structure of rrnL of O. stentina is presented in Figure 6 , while a detailed representation of the 2D structure is available in Supplementary Figure S31 . The structure mirrored those available for other oyster species (Milbury et al., 2010) and more in general molluscs (Lydeard et al., 2000), with domain I and II located in the 5’ half, domain III lacking, and domain IV-VI located in the 3’ half. The rrnL structures inferred for other Ostreinae overlapped with that presented here for O. stentina ( Figure 6 ). Among the Ostreinae, most of the variable positions in rrnL were located in the 5’ half (average p-distance = 0.281 ± 0.112) while the 3’ half was much more conserved (average p-distance = 0.148 ± 0.072). This higher level of conservation reflects the prominent structural role of the 3’ half for the functioning of the whole rrnL molecule (Lydeard et al., 2000; Milbury et al., 2010). Furthermore, the 3’ half was markedly GC-richer than the 5’ half ( Supplementary Figure S14 ), and the stability of the stems in its highly conserved domains IV-V was often guaranteed by the pairs G-C and C-G ( Figure 6 ; Supplementary Figure S31 ). Intraspecific behavior mirrored that observed in the comparisons among different species of Ostreinae. The 5’ half was more variable than the 3’ segment in all tree species ( Supplementary Figures S32–S34 ). However, the level of variation was very limited or non-existent as in case of the 3’ halves of O. denselamellosa ( Supplementary Figure S34B ).

In the mitogenome of Ostreinae, the control region was located between trnD and cox1 ( Figure 2 ) and, for the first time, was characterized in detail in this paper for this subfamily of oysters. Its size varied from 688 bp (O. denselamellosa ON964460) to 742 bp (N. pinnicola). As stated above, the CoRe was extremely AT-rich ( Supplementary Figure S15 ). The sequences of CoRe were highly diverging, as proved by the very low number of fully conserved positions in their multiple alignment ( Supplementary Figure S35 ). The alignment showed large portions that were difficult to align with high accuracy (Chang et al., 2014), even when using a highly sophisticated software such as T-coffee (Di Tommaso et al., 2011). Despite the high variability, two segments appeared rather conserved in the alignment: one spanning positions 70 to 130, and the other ranging from positions 480 to 570 ( Supplementary Figure S35 ). Two fully conserved motifs were present in the CoRes of Ostreinae. The first one (AAAGGGG) started at position 171 of the alignment ( Supplementary Figure S35 ). This motif was present also in rrnS. It occurred in the mitogenome of other Osteidae (Saccostrea, Magallana and Talonostrea), but not in their CoRe. A second fully conserved motif of 11 nucleotides (CTATGTAAATA) extended from position 552 to position 562 ( Supplementary Figure S35 ). This motif was exclusive of the CoRe of Ostreinae sequenced to date and did not occur in the mitogenomes of other oysters. The CoRe of O. denselamellosa presented a second copy of this motif (positions: 132-142, Supplementary Figure S35 ).

The CoRes of O. stentina ranged from 701 bp to 711 bp ( Supplementary Table S1 , Supplementary Figure S35 ) and their average p-distance was 0.010 ± 0.008. The four CoRes of O. edulis ranged from 695 bp to 700 bp, and two were identical (JF274008 and CM063324). Their average p-distance was 0.022 ± 0.023. Additionally, the two CoRes of O. denselamellosa differed in length by one nucleotide (688 vs. 689) ( Supplementary Figure S35 ) and their p-distance was 0.017. A very peculiar case was observed for D. sandvichensis and P. pestigris, where the available mitogenomes exhibited identical CoRes.

Stretches of polyA, polyT and polyG, as well as polyAT, characterized the CoRes of Ostreinae ( Supplementary Figure S35 ). These features are peculiar of, and specific to, the control regions of molluscs and, more in general, animals (Ghiselli et al., 2021). Finally, all CoRes of Ostreinae were able to form stem-and-loop secondary structures ( Figure 7 ), and these structures were located in a highly variable portion of their multiple alignment. These structures are considered important for the replication and transcription of mitogenomes (Ghiselli et al., 2021).

The level of compositional heterogeneity occurring among orthologous sequences was evaluated for all 13 PCGs and rDNAs genes with the software AliGROOVE (Kück et al., 2014) ( Supplementary Figures S36–S43 ). The third codon positions of all PCGs exhibited high heterogeneous substitution patterns, while the most homogenous single positions were the second positions of several PCGs (i.e. cob, cox1-cox2, nad1, nad3-nad5). Amino acid sequences were very homogeneous in their substitution patterns with rare exceptions observed in nad2 and nad6 ( Supplementary Figures S39, S42 ). Substitution patterns for rrnL and rrnS were homogeneous among Ostreinae species but were heterogeneous compared with the sequences of other subfamilies ( Supplementary Figure S42 ). Additionally, when the 13PCG data sets were considered ( Table 2 ), the substitution process was homogenous for amino acids, the first (mostly) and second positions of codons, as well as for first + second positions, plus a large part of whole codons. On the contrary, the substitution pattern was heterogeneous for third positions, except in intraspecific comparisons ( Supplementary Figure S43 ).

The phylogenetic signal for single PCG was the highest for amino acids and first + second codon positions, while third positions appeared highly saturated, as evidenced by the maximum likelihood distances, largely exceeding 1, and the lowest percentage of fully resolved quartets ( Supplementary Figures S44–S46 ). Both rrnL and rrnS exhibited a good phylogenetic signal. Similarly, the best signal among the concatenated alignments was observed for 13PCGpro data set and 13PCGp12 data set ( Table 2 ; Supplementary Figure S46 ), which contained respectively the amino acid sequences and the first + second position of the 13 PCGs.

13PCGpro exhibited the best signal and the most homogeneous substitution pattern among the analyzed data sets (see above) ( Table 2 ). The ML tree (hereafter Tree 1) obtained from this set is provided in Figure 8 ( Supplementary Table S7 ). Most of the nodes and branches received very strong statistical corroboration. Within Ostreinae, the genus Ostrea appeared polyphyletic, with O. stentina sister species of O. lurida and O. edulis + O. denselamellosa forming a separated group nested within a second clade encompassing N. pinnicola + P. plestigris, their sister taxon, and D. sandvichensis. UBT/aLRT values strongly support this clade. 13PCGpro+rRNAs, 13PCGp2, 13PCGp12.a/b, and 13PCGp12+rRNAs.a/b produced also Tree 1 ( Supplementary Figures S47–S52 , and Supplementary Table S7 ).

Some data sets listed in Table 2 generated trees that differed from tree 1( Supplementary Figures S53–S59 , and Supplementary Table S7 ). However, these alternative topologies lacked strong statistical support. We will analyze Tree 3, 4, and 6 ( Supplementary Figures S54, S55, S57 ) in more detail here, as the Ostrea genus resulted monophyletic, although the most basal node, the critical one, did not receive strong statistical support. All these trees were the product of the analyses performed on data sets including third positions of codons and/or rRNAs ( Table 2 ; Supplementary Table S7 ). The third positions of codons exhibited a highly heterogeneous substitution pattern and their phylogenetic signal was mostly/completely lost, two factors that are highly detrimental for phylogenetic analyses (Negrisolo et al., 2004; Kück et al., 2014). Furthermore, they failed a test of stationarity or homogeneity (p < 0.05) computed with IQ-TREE2, raising serious concerns about the reliability of the trees derived from their analyses (Naser-Khdour et al., 2019). Tree 6 was the phylogenetic output of the rRNAs data set. The ribosomal markers exhibited good phylogenetic signals. However, the substitution pattern was not homogeneous between ingroup and outgroup sequences ( Supplementary Figure S42 ), a factor that can influence the phylogenetic outputs (Kück et al., 2014).

The results of alternative topologies tests performed on the data sets analyzed in the present paper are summarized in Supplementary Table S8 . Tree 1, our reference topology, was rejected only by 13PCGp3, the least reliable analyzed data set. Conversely, the most robust data set, i.e. 13PCGpro, rejected nearly all alternative topologies, except for Tree 4 in the AU test and Trees 4 and 5 in the highly conservative WSH test (Shimodaira, 2002).

The mitogenomes of O. stentina contained 28 ISPs, ranging from 1 (ISP trnL2-trnP) to 117-118 (ISP trnG-cox3) nucleotides ( Figure 1 ; Supplementary Table S1 ). No genes overlapped. Similar patterns characterized the mitogenomes of Ostreinae sequenced to date ( Figure 9 ). In several pair of consecutive genes, the behavior is fixed: they were either separated by an ISP (e.g. cox1 and trnG, trnK and trnL1), or adjacent (e.g. trnH and nad4, nad3 and trnK). In other cases, the pattern changed in different species. In particular, trnL1 and trnF were spaced in O. edulis, O. denselamellosa, N. pinnicola, and P. plestigris, and adjacent in other oysters ( Figure 9 ; Supplementary Figure S60G ).

We analyzed the sequences of the ISPs by aligning them manually or with MAFFT. Notably, in 19 ISPs we identified sequences ( Figure 8 ) that were exclusive to and characterized a single species or group of species located downstream to a well-supported node of the reference Tree 1 (e.g. ISP trnG-cox3; O. lurida + O. stentina) ( Figure 8 ; Supplementary Figure S60 ).

Phylogenetic relationships among the four subfamilies of Ostreidae are well defined ( Figure 2B ) (Salvi and Mariottini, 2017; Li et al., 2021; Salvi and Mariottini, 2021). We identified gene blocks shared among the different subfamilies and, considering their phylogenetic relationships, we were able to partly infer the gene order arrangement of the mitogenome of the last common ancestor (lcaO) of all Ostreidae ( Figure 2B ). Saccostreinae and Ostreinae are not sister taxa, but share two large gene blocks, which represent a plesiomorphic condition for the family Ostreidae ( Figure 2B ). Similarly, the four blocks shared among the three subfamilies Crassostreinae, Ostreinae and Saccostreinae represent further plesiomorphies ( Figure 2B ). The complete transformational pathway that lead to the diversity of GOs observed today in oysters, particularly in Crassostreinae, remains to be fully understood. Sequencing the mitogenomes of Striostreinae, the sister group of Crassostreinae ( Figure 2B ), and the only subfamily without available mitochondrial sequences, is a priority to properly address this issue.

Standard evolutionary models used in phylogenetic analyses assume that the substitution process among orthologous sequences in the multiple alignments is homogeneous and the violation of this assumption may generate misleading phylogenetic outputs (Kück et al., 2014). The lack of phylogenetic signal is another important source of distortive effects on phylogenetic results. We performed the quartet puzzling analysis (Strimmer and von Haeseler, 1996) and analyzed the distribution of the pairwise distances computed according to the best-fit evolutionary model to test this amount of phylogenetic signal (e.g. Negrisolo et al., 2004). It is well known that when the distribution of these distances is considerably greater than one, there is a substantial loss of phylogenetic signal in the analyzed dataset (e.g. Negrisolo et al., 2004). The best phylogenetic markers, when working at the taxonomic level of the family, proved to be the proteins. The 13PCGpro multiple alignment exhibited the best signal and the most homogeneous substitution pattern among the analyzed data sets. In contrast, third positions of codons showed very heterogeneous substitution patterns and substantial lack of phylogenetic signal.

The genus Ostrea was polyphyletic in our reference tree ( Figure 8 ). We obtained also alternative topologies implying the monophyly of this taxon, but these were obtained from data sets that proved to be unreliable markers (third positions) or difficult to manage (ribosomal genes) due to non-homogeneous substitution patterns occurring between ingroup and outgroup (Kück et al., 2014). Ostrea resulted para/polyphyletic also in previous phylogenetic analyses based on both nuclear and mitochondrial genes (Salvi et al., 2014; Guo et al., 2018; Li et al., 2021; Salvi and Mariottini, 2021), suggesting that homoplasy characterizes the morphological evolution of the genus. Our results strongly support a polyphyletic nature for Ostrea. However, we worked with a limited taxon sampling. Therefore, a wider species coverage is necessary for corroborating this point.

The atp8 was the most variable PCG, followed by nad2 and nad6, while the most conserved was cox1. This result is very interesting and supports the hypothesis that the mitogenomes of oysters contain multiple PCGs that can be used for molecular identification of the species outside of cox1 and further corroborates earlier findings (e.g. Xiao et al., 2015). At the intraspecific level, variability was limited in O. stentina. However, all mitogenomes were obtained from specimens collected in the Venice Lagoon, thus they do not represent the global diversity of the species.

The mitochondrial tRNAs harbor a considerable amount of taxonomic and evolutionary information that fully stands out when their secondary structure is considered (e.g. Simonato et al., 2013; Montelli et al., 2016). Unfortunately, these markers are often given only a cursory treatment. In this study, we analyzed in details the substitution process characterizing the multiple alignments of orthologous tRNAs. Particularly interesting are the base changes occurring in the stems of tRNAs (Coleman, 2003; Montelli et al., 2016). In our study, the tRNAs associated with the most abundant codon families were the least variable. This pattern of conservation supports the hypothesis that these tRNAs have a more constrained nucleotide substitution pattern, associated to their high frequency of usage in the protein synthesis.

Some tRNAs exhibited AT-/GC- skews values that differed greatly from those of the strand encoding them. This is not unique to oysters. A similar pattern was observed in the tRNAs of Cetacea (Montelli et al., 2016). It was not possible to identify a single cause (e.g. tRNAs associated to abundant amino acids) that explained this result. The short length of tRNAs likely played a role, as even a small number of substitutions can have a strong impact on their skew values.

The hyper-variable portions of DHU loop, “extra arm and TΨC loop of several mitochondrial tRNAs exhibited sequence motifs that characterized single/group of species of oysters. Fully compensatory base changes, as well as mismatches, were also present in the stems of tRNAs, either restricted to single oyster or, conversely, exclusive to the entire subfamily Ostreinae. Our taxon coverage is very sparse, but despite this limitation, these tRNAs features could serve as additional taxonomic tools for the family Ostreidae, where identification of species and taxa relationships are problematic (Harry, 1985), as observed in other groups of invertebrates (e.g. Simonato et al., 2013).

For rrnS, a secondary structure model did not exist for Ostreidae and we did not attempt to develop a new one. In contrast, we used the secondary structure models of rrnL available for the phylum Mollusca (Lydeard et al., 2000) and for the family Ostreidae (Milbury et al., 2010) to infer the secondary structure of the rrnL of O. stentina. The analyses of compositional biases and AT-/GC-skews of rrnLs and rrnSs suggests that structural constraints played a key role in shaping these features.

The CoRe of all analyzed mitogenomes contained the peculiar sequence motif CTATGTAAATA. If this motif is found to be exclusive to all Ostreinae, it might become a very useful marker to unambiguously identify this genomic portion, similar to other motifs identified in various animal groups (e.g. Lepidoptera; Salvato et al., 2008). A very peculiar case was observed for D. sandvichensis and P. pestigris, where the available mitogenomes exhibited identical CoRes. These sequences were produced by the same research group at different times. As shown above, CoRes are variable at the intraspecific level. Furthermore, D. sandvichensis and P. pestigris are not sister species ( Figure 8 ). Therefore, the occurrence of an identical control region in their mitogenomes requires independent confirmation.

Many intergenic spacers located throughout the mitogenome ( Figure 9 ) contain sequences characteristic of a single species or clade. These sequences are mito-signatures (Liu et al., 2022), i.e. molecular markers useful to define/identify taxa in a phylogenetic context, but cannot be considered true synapomorphyes. Uniqueness is the hallmark of a true apomorphy (Page and Holmes, 2009). However, it is very unlikely that an often short sequence of ISP could fulfill this stringent requirement. A mito-signature can be very useful to identify a species, a group of species, or even bigger taxa, within a well-established phylogenetic framework. This is particularly relevant in animals like oysters, as they are difficult to identify on a morphological basis (e.g. Harry, 1985). In the past, the value of ISP as intraspecific phylogenetic markers has been shown in the Crassostreinae (Ren et al., 2016). Our findings further corroborate this point and extend, at the interspecific level, the taxonomic/phylogenetic value of these short sequences for oysters, as already known in other groups of animals (e.g. Simonato et al., 2013; Basso et al., 2017; Liu et al., 2022).