Chromosome-Level Genome Assembly of the Blue Mussel Mytilus chilensis Reveals Molecular Signatures Facing the Marine Environment

The blue mussel Mytilus chilensis is an endemic and key socioeconomic species inhabiting the southern coast of Chile. This bivalve species supports a booming aquaculture industry, which entirely relies on artificially collected seeds from natural beds that are translocated to diverse physical–chemical ocean farming conditions. Furthermore, mussel production is threatened by a broad range of microorganisms, pollution, and environmental stressors that eventually impact its survival and growth. Herein, understanding the genomic basis of the local adaption is pivotal to developing sustainable shellfish aquaculture. We present a high-quality reference genome of M. chilensis, which is the first chromosome-level genome for a Mytilidae member in South America. The assembled genome size was 1.93 Gb, with a contig N50 of 134 Mb. Through Hi-C proximity ligation, 11,868 contigs were clustered, ordered, and assembled into 14 chromosomes in congruence with the karyological evidence. The M. chilensis genome comprises 34,530 genes and 4795 non-coding RNAs. A total of 57% of the genome contains repetitive sequences with predominancy of LTR-retrotransposons and unknown elements. Comparative genome analysis of M. chilensis and M. coruscus was conducted, revealing genic rearrangements distributed into the whole genome. Notably, transposable Steamer-like elements associated with horizontal transmissible cancer were explored in reference genomes, suggesting putative relationships at the chromosome level in Bivalvia. Genome expression analysis was also conducted, showing putative genomic differences between two ecologically different mussel populations. The evidence suggests that local genome adaptation and physiological plasticity can be analyzed to develop sustainable mussel production. The genome of M. chilensis provides pivotal molecular knowledge for the Mytilus complex.


Introduction
The blue mussel M. chilensis (Hupé, 1854) is Chile's endemic, ecological, and socioeconomic key species that leads the national shellfish aquaculture. The farming of M. chilensis Moreover, RNA libraries were constructed from hemocytes, digestive gland, gill, and mantle tissues for transcriptome sequencing to obtain whole-transcriptome profiling from the same mussels used for genome DNA sequencing. Additionally, twelve available Sequence Read Archive (SRA) transcriptomic data (GenBank accession number SRP261955), representing gills and mantle tissues collected from individuals of Cochamó and Yaldad mussel populations [37,38], were incorporated to analyze population-specific transcriptome profiles. Total RNA from three biological replicates (five total RNA extractions each) from each mussel population was extracted by the Trizol reagent method (Invitrogen, Waltham, MA, USA). The quality and integrity of extracted RNAs were measured in a Tape Station 2200 instrument (Agilent, Santa Clara, CA, USA), using the R6K Reagent Kit based on the manufacturer's instructions. RNA samples with RNA integrity numbers (RIN) >9 were selected for the preparation of high-quality libraries using a TrueSeq Stranded mRNA LT Sample Prep Kit and sequenced in a HiSeq 4000 (Illumina, San Diego, CA, USA).

De Novo Genome Assembly and Hi-C Scaffolding of M. chilensis
Two HiFi single-molecule real-time cells in the PacBio Sequel platform yielded 53.8 Gb of high-quality DNA genome information. These long reads were assembled with the Hifiasm (v.0.19.3) package using default parameters [45]. For Hi-C scaffolding, reads were aligned to the primary draft assembly, following the manufacturer's instructions [46]. Briefly, reads were aligned using BWA-MEM (v.0.7.17) [47] with the -5SP and -t 8 options specified and all other options default. The package SAMBLASTER (v.0.1.26) [48] was used to flag duplicates excluded from further analysis. Sequence alignments were filtered with SAMtools (v. 1.17) [49,50] using the -F 2304 filtering flag to remove non-primary and secondary alignments. This step was conducted to remove alignment errors, low-quality alignments, and other alignment noise due to repetitiveness, heterozygosity, and other ambiguous assembled sequences. Finally, Phase Genomics' Proximo Hi-C genome-scaffolding platform was used to create chromosome-scale scaffolds from the corrected assembly, according to Bickhart et al. [51]. Moreover, RNA libraries were constructed from hemocytes, digestive gland, gill, and mantle tissues for transcriptome sequencing to obtain whole-transcriptome profiling from the same mussels used for genome DNA sequencing. Additionally, twelve available Sequence Read Archive (SRA) transcriptomic data (GenBank accession number SRP261955), representing gills and mantle tissues collected from individuals of Cochamó and Yaldad mussel populations [37,38], were incorporated to analyze population-specific transcriptome profiles. Total RNA from three biological replicates (five total RNA extractions each) from each mussel population was extracted by the Trizol reagent method (Invitrogen, Waltham, MA, USA). The quality and integrity of extracted RNAs were measured in a Tape Station 2200 instrument (Agilent, Santa Clara, CA, USA), using the R6K Reagent Kit based on the manufacturer's instructions. RNA samples with RNA integrity numbers (RIN) >9 were selected for the preparation of high-quality libraries using a TrueSeq Stranded mRNA LT Sample Prep Kit and sequenced in a HiSeq 4000 (Illumina, San Diego, CA, USA).

De Novo Genome Assembly and Hi-C Scaffolding of M. chilensis
Two HiFi single-molecule real-time cells in the PacBio Sequel platform yielded 53.8 Gb of high-quality DNA genome information. These long reads were assembled with the Hifiasm (v.0.19.3) package using default parameters [45]. For Hi-C scaffolding, reads were aligned to the primary draft assembly, following the manufacturer's instructions [46]. Briefly, reads were aligned using BWA-MEM (v.0.7.17) [47] with the -5SP and -t 8 options specified and all other options default. The package SAMBLASTER (v.0.1.26) [48] was used to flag duplicates excluded from further analysis. Sequence alignments were filtered with SAMtools (v. 1.17) [49,50] using the -F 2304 filtering flag to remove non-primary and secondary alignments. This step was conducted to remove alignment errors, lowquality alignments, and other alignment noise due to repetitiveness, heterozygosity, and other ambiguous assembled sequences. Finally, Phase Genomics' Proximo Hi-C genomescaffolding platform was used to create chromosome-scale scaffolds from the corrected assembly, according to Bickhart et al. [51].

Karyotype of M. chilensis
Metaphase plates of 24-h-post-fertilization larvae were used to obtain chromosomes from M. chilensis, according to Gallardo-Escárate et al. [52]. Briefly, antimitotic treat-Genes 2023, 14, 876 5 of 27 ment with colchicine 0.05% solution was applied for 4 h. Then, the larvae were rinsed in clean seawater and immersed in a hypotonic solution (seawater:distilled water, 1:1) for 30 min. Finally, the larvae were fixed in modified Carnoy solution (methanol:acetic acid, 3:1). Chromosome spreads were obtained by dissociating larva tissue in acetic acid (50%), pipetting suspension drops onto slides preheated at 43 • C and air-dried according to Amar et al. [53]. A FISH experiment was performed to validate the physical localization of specific genes. Here, 28S rDNA was labeled following methods previously published [21]. Briefly, metaphase preparations were denatured at 69 • C for 2 min and hybridized overnight at 37 • C. Signal detection was performed using fluorescein avidin and biotinylated antiavidin for the biotinylated probes, and mouse anti-digoxigenin, goat anti-mouse rhodamine, and rabbit anti-goat rhodamine for the digoxigenin-labeled probes. Fluorescent staining was carried out with 4,6-diamidino-2-phenylindole (DAPI) and mounted with Vectacshield antifading solution. Chromosome spreads were observed using an epifluorescent microscope Nikon Eclipse 80i (Minato-ku, Tokyo, Japan) equipped with a digital camera DS-5Mc.

Genome Annotation of M. chilensis
Our repeat annotation pipeline applied a combined homology alignment strategy, and de novo search to identify the whole-genome repeats. Tandem repeat was extracted using TRF (v.4.09.1) (http://tandem.bu.edu/trf/trf.html (accessed on 3 October 2021)) by ab initio prediction.
The structural annotation approach was applied to incorporate de novo, homolog prediction, and RNA-Seq-assisted predictions to annotate gene models. For gene prediction based on de novo, Augustus (v3.2.3), Geneid (v1.4), Genescan (v1.0), GlimmerHMM (v3.04), and SNAP (29 November 2013) were used in our automated gene prediction pipeline. Sequences of homologous proteins were downloaded from Ensembl/NCBI/others for homolog prediction. Protein sequences were aligned to the genome using TblastN (v2.2.26; E-value ≤ 1 × 10 −5 ), and then the matching proteins were aligned to the homologous genome sequences for accurate spliced alignments with GeneWise (v2.4.1) software to predict the gene structure contained in each protein region. Finally, for RNA-seq data, transcriptome reads assemblies were generated with Trinity (v2.1.1) for the genome annotation. For the genome annotation optimization, the RNA-Seq reads from different tissues were aligned to genome fasta using Hisat (v2.0.4)/TopHat (v2.0.11) with default parameters to identify exons region and splice positions. The alignment results were inputted into Stringtie (v1.3.3)/Cufflinks (v2.2.1) with default parameters for genome-based transcript assembly. The non-redundant reference gene set was generated by merging genes predicted by three methods with EvidenceModeler (EVM v1.1.1) using PASA (Program to Assemble Spliced Alignment) terminal exon support and including masked transposable elements as input into gene prediction. Individual families of interest were selected for further manual curation.
Gene functions were assigned according to the best match by aligning the protein sequences to the Swiss-Prot database using Blastp (with a threshold of E-value ≤ 1 × 10 −5 ). The motifs and domains were annotated using InterProScan70 (v5.31) by searching against publicly available databases, including ProDom, PRINTS, Pfam, SMRT, PANTHER, and PROSITE. Each gene's gene ontology (GO) IDs were assigned according to the corresponding InterPro entry. We predicted the protein function by transferring annotations from the closest BLAST hit (E-value < 1 × 10 −5 ) in the Swiss-Prot database and DIAMOND (v0.8.22)/BLAST hit (E-value < 1 × 10 −5 ) in the NR database. We also mapped the gene set to a KEGG pathway and identified the best match for each gene.

Comparative Genomics between M. chilensis and M. coruscus
Syntenic relationships were explored among mussel species for which chromosomelevel reference genomes are publicly available. Here, the analysis was performed between the two congeneric species M. chilensis (this study) and M. coruscus [28], where gene annotations were explored by MCScanX (v. 1.0) [56] implemented in the TBtools (v.1.115) package [57]. This approach detects groups of orthologous genes and compares their arrangement to identify colinear segments in the compared genomes. MCScanX was used to discover microsyntenic relationships, focusing on the local arrangement of genes near the syntenic blocks. The microsyntenic arrangement of genes identified by MCScanX was evaluated through GO analysis to identify the primary molecular function and biological processes enrichened for each genomic region where macromutations or chromosome rearrangements were detected.
Disseminated neoplasia is a disease horizontally transmitted by clonal cancer cells, which causes leukemia in mollusk bivalves [58,59]. The neoplastic cells gradually replace normal hemocytes leading to increased mortality, and the disease has been detected in 15 species of marine bivalve mollusks worldwide [60]. Notably, disseminated neoplasia has been observed among mussel species with varying epizootic prevalences. For instance, M. trossulus has shown high prevalences in some areas, whereas in Mytilus edulis, the prevalences are generally lower. Furthermore, M. galloprovincialis has been suggested as a species resistant to the disease in Spanish and Italian mussel populations [61]. This observation extends the relevance of exploring mussel species' genetic features associated with disseminated neoplasia. Herein, the molecular characterization of Steamer-like elements in M. chilensis was conducted by cloning and the walking primer method according to Arriagada et al. [62]. The putative M. chilensis Steamer-like was scanned through twelve reference genomes assembled at chromosome level for Bivalvia: Mytilus coruscus (GCA_017311375.1), Mytilus edulis (GCA_019925275.1), Dreissena polymorpha (GCA_020536995.1), Mercenaria mercenaria (GCA_014805675.2), Solen grandis (GCA_021229015.1), Ruditapes philippinarum (GCA_009026015.1), Pecten maximus (GCA_902652985.1), Pinctada imbricata (GCA_002216045.1), Crassostrea gigas (GCA_902806645.1), Crassostrea ariakensis (GCA_020458035.1), and Crassostrea virginica (GCA_002022765.4). The putative long terminal repeat (LTR) sequences were identified using BLAST search, where open reading frames (ORFs) between flanking LTRs were detected. The identified Steamer-like elements were aligned using ClustalW and annotated based on a search of the NCBI Conserved Domain database. Amino acid sequences for the full-length Gag-Pol polyprotein region were aligned among the studied bivalve species. The Steamer element was reported for Mya arenaria (Accession AIE48224.1) and M. chilensis (this study). DNA sequence genealogy analysis was conducted to investigate horizontal transmission events among bivalve species. The maximum likelihood (ML) method was conducted on the SLEs loci localized in all the publicly available bivalve genomes assembled at the chromosome level.

Whole-Genome Transcript Expression Analysis in Two M. chilensis Populations
The transcriptomes of mussels collected from the Yaldad and Cochamó populations were analyzed using a hierarchical clustering approach to detect transcriptional similarities among tissues/populations. The differentially expressed transcripts compared to normalized expression values were visualized in a clustering heatmap and selected according to the identified cluster. For an optimal comparison of the results, k-means clustering was performed to identify candidate genes involved in specific gene expression patterns. The distance metric was calculated with the Manhattan method, where the mean expression level in 5-6 rounds of k-means clustering was subtracted.
Moreover, raw data from mussel tissues collected from the Yaldad and Cochamó populations were trimmed and mapped to the M. chilensis genome using CLC Genomics Workbench v22 software (Qiagen Bioinformatics, Aarhus, Denmark). Threshold values for transcripts were calculated from the coverage analysis using the Graph Threshold Areas tool in CLC Genomics Workbench v22 software. Here, an index denoted as chromosome genome expression (CGE) was applied to explore the whole-genome transcript expression profiling, according to Valenzuela et al. [63]. The CGE calculates the mean coverage of transcripts mapped into a specific chromosome region, comparing mussel populations and tissues. Specifically, the CGE index represents the percentage of the transcriptional variation between two or more RNA-seq data for the same locus. The transcript coverage values for each dataset were calculated using a threshold of 20,000 to 150,000 reads. A window size of 10 positions was set to calculate and identify chromosome regions differentially transcribed. This approach was used to visualize actively transcribed chromosome regions, identify differentially expressed genes, and observe tissue-specific patterns in the evaluated mussel populations. Finally, the threshold values for each dataset and the CGE index were visualized in Circos plots [64].
RNA-seq data analyses were carried out using the raw sequencing reads and mapped on the assembled genome by CLC Genomics Workbench v22 software (Qiagen Bioinformatics, Aarhus, Denmark) separately for each tissue/population. In parallel, de novo assembling was performed to evaluate PAVs and dispensable genes affecting the in-silico transcription analysis. The assembly was performed with overlap criteria of 70% and a similarity of 0.9 to exclude paralogous sequence variants. The settings were set as mismatch cost = 2, deletion cost = 3, insert cost = 3, minimum contig length = 200 base pairs, and trimming quality score = 0.05 using CLC Genomics Workbench v22. After assembly, the contigs generated for each data set were mapped on the genes annotated in the reference genome to evaluate genome coverage and detect PAV features. The analysis did not show bias putatively associated with PAVs between the analyzed mussel populations. Then, mRNA sequences annotated for the M. chilensis genome were used to evaluate the transcription level between mussel populations, where differential expression analysis was set with a minimum length fraction = 0.6 and a minimum similarity fraction (long reads) = 0.5. The obtained genes from each tissue/population were blasted to CGE regions to enrich the number of transcripts evaluated by RNA-Seq analysis. In addition, sequences were extracted near the threshold areas in a window of 10 kb for each transcriptome. The expression value was set as transcripts per million model (TPM). The distance metric was calculated with the Manhattan method, with the mean expression level in 5-6 rounds of k-means clustering subtracted. Finally, the Generalized Linear Model (GLM) available in the CLC software was used for statistical analyses and to compare gene expression levels regarding the log 2 fold change (p = 0.005; FDR corrected).
Moreover, innate immunity in marine invertebrates may play an important role in speciation and environmental adaptation [65,66]. Herein, we investigate the immunerelated genes associated with the Toll-like receptor (TLR) and apoptosis signaling pathways given that the functional annotation revealed that they were mainly enriched between the mussel populations analyzed. In addition, bioinformatic analyses were carried out using the CLC Genomics Workbench software to identify single nucleotide variants (SNV) from the transcriptomes sequenced for Yaldad and Cochamó. Candidate SNVs were called with the following settings: window length = 11, maximum gap and mismatch count = 2, minimum average quality of surrounding bases = 15, minimum quality of central base = 20, maximum coverage = 100, minimum coverage = 8, minimum variant frequency (%) = 35.0, and maximum expected variations (ploidy) = 2. Furthermore, the genotypes of DEGs were also identified for detecting putative genetic variations between mussel populations. Here, singleton, dispersed, tandem, proximal, and whole-genome duplication (WGD) gene events were evaluated using MCScanX. The amino acid changes and the zygosity proportions were also estimated in DEGs between the Yaldad and Cochamó populations.

GO Enrichment Analysis
Differentially expressed mRNA were annotated through BlastX analysis using a custom protein database constructed from GeneBank, KEGG, GO, and UniProtKB/Swiss-Prot. The cutoff E-value was set at 1 × 10 −10 . Transcripts were subjected to gene ontology (GO) analysis using the Blast2GO plugins included in the CLC Genomics Workbench v22 software (Qiagen Bioinformatics, Aarhus, Denmark). The results were plotted using the Profiler R package [67]. GO enrichment analysis was conducted to identify the most represented biological processes among protein-coding genes proximally located to the CGE regions. The enrichment of biological processes was identified using Fisher's exact test tool of Blast2GO among the different tissues and mussel populations.

Chromosome Genome Assembly of M. chilensis Using Proximity Ligation
With two HiFi single-molecule real-time cells in the PacBio Sequel platform, we generated 53.8 Gb of high-quality DNA genome information.
This data comprised 63 million reads with a total length of 882 Gbp (Table 1). These long reads were assembled with the Hifiasm package using default parameters [45], yielding a primary assembly of 13,762 contigs equivalent to 2.19 Gb, with an N50 of 206 Mb. The genome size assembly made by Hifiasm was comparable with the previous genome size described for closely related species: 1.28 Gb for M. galloprovincialis [32], 1.57 Gb for M. coruscus [28], and 1.79 Gb for Dreissena polymorpha [27]. Table 1. Statistics of whole-genome sequence assembly and transcriptome analysis of the blue mussel Mytilus chilensis using Illumina, PacBio, and Hi-C. In vivo Hi-C is a technique that maps physical DNA-DNA proximity across the entire genome [68,69]. The method was introduced as a genome-wide version of its predecessor, 3C (chromosome conformation capture). It has been a powerful tool in chromosome-scale genome assembly of many animals in recent years [70,71]. In this study, Hi-C experiments and data analysis of hemocyte cells were used for the chromosome assembly of the blue mussel M. chilensis. Here, Phase Genomics (Seattle, WA, USA) prepared and sequenced two Hi-C libraries, resulting in~20× coverage and~253 million 150 bp paired-end reads ( Table 1). The Hi-C analysis evidenced that 44.68% of high-quality reads showed intercontig signals or Cis-close position (<10 kbp on the same contig), and an additional 4.09% of sequence reads revealed a Cis-far conformation (>10 kbp on the same contig) ( Table 2). Hi-C reads were aligned using Bowtie version 1.3.1 [72] to order and orient the 13,762 contigs, and scaffolding was performed using Proximo (Phase Genomics, Seattle, WA, USA). We then applied Juicebox for visual inspection and manual correction [73]. We also manually removed 1894 scaffolds that were microbe-sized and disconnected from the rest of the assembly. Then, 11,868 contigs were used for the first chromosome-level high-quality M. chilensis assembly (Table 3). The N50 and total genome length were calculated in 134 Mbp and 1938 Gbps, respectively. The M. chilensis genome provides a valuable genomic resource for research in mussel biology and for developing novel sustainable strategies in mussel aquaculture. The Hi-C data generated 14 chromosomes assembled with HiFi consensus long DNA reads ( Figure 1B). The cytogenetic analysis performed for M. chilensis revealed a conservative karyotype for the Mytilus genus composed of 2n = 14 [21]. Physical localization of 28S-rRNA revealed two loci mapped in different submetacentric/subtelocentric chromosome pairs ( Figure 1C), confirming the presence of major rDNA clusters subterminal to the long arms of two chromosome pairs reported in M. edulis and M. galloprovincialis [74]. Concerning genome assembly, the largest scaffold was assembled from 998 contigs with a total size of 173.3 Mb. Meanwhile, the smallest scaffold was 117.3 Mb, consisting of 744 contigs (Table 3). Notably, the number of contigs in the scaffolds was 11,868 (100% of all contigs in chromosome clusters, 86.24% of all contigs) and accounted for 1.93 Gbps of genome size (100% of all length in chromosome clusters, 88.43% of all sequence length). The completeness of genome assembly was assessed by the single-copy ortholog set (BUSCO, V5.3.2) [75]. The results showed the following BUSCO scores: (i) Eukaryota Odb10; C:94.1% (S:72

Genome Annotation of M. chilensis
The genome assembly was annotated using de novo and protein-and transcriptguided methods (Figure 2A). The first step of the annotation process was to identify the DNA repeats through the M. chilensis genome. Repetitive elements and non-coding genes in the blue mussel genome were annotated by homologous comparison and ab initio prediction. RepeatMasker [76] was used for homologous comparison by searching against the Repbase database [77] and RepeatModeler [78]. According to these analyses, about 1.1 Gbps of repeat sequences were annotated, which accounted for 56.73% of the whole genome. Herein, DNA transposons, LINE, and LTR transposable elements were identified (Table 4). Useful genome information for population genetic studies is the identification of simple sequence repeats (SSRs) or microsatellites. The mining of SSRs revealed that the M. chilensis genome has 548,360 SSR sequences, where 9% and 6% of the SSR loci were annotated for each mussel chromosome ( Figure S1). The most frequent SSR motif was the tetranucleotide, followed by the dinucleotides, accounting for 206,103 and 197,700 repeats, respectively. The entire SSR sequences accounted for 0.35% of the whole genome. The development of SSR markers offers a shortcut to assessing genetic diversity, which can potentially be applied in food authentication and genetic traceability for mussel species [79][80][81].
pairs ( Figure 1C), confirming the presence of major rDNA clusters subterminal to the long arms of two chromosome pairs reported in M. edulis and M. galloprovincialis [74]. Concerning genome assembly, the largest scaffold was assembled from 998 contigs with a total size of 173.3 Mb. Meanwhile, the smallest scaffold was 117.3 Mb, consisting of 744 contigs (Table 3). Notably, the number of contigs in the scaffolds was 11,868 (100% of all contigs in chromosome clusters, 86.24% of all contigs) and accounted for 1.93 Gbps of genome size (100% of all length in chromosome clusters, 88.43% of all sequence length). The completeness of genome assembly was assessed by the single-copy ortholog set (BUSCO, V5.3.2) [75]. The results showed the following BUSCO scores: (i) Eukaryota Odb10; C:94.1% (S:72

Genome Annotation of M. chilensis
The genome assembly was annotated using de novo and protein-and transcriptguided methods (Figure 2A). The first step of the annotation process was to identify the DNA repeats through the M. chilensis genome. Repetitive elements and non-coding genes in the blue mussel genome were annotated by homologous comparison and ab initio prediction. RepeatMasker [76] was used for homologous comparison by searching against the Repbase database [77] and RepeatModeler [78]. According to these analyses, about 1.1 Gbps of repeat sequences were annotated, which accounted for 56.73% of the whole genome. Herein, DNA transposons, LINE, and LTR transposable elements were identified (Table 4). Useful genome information for population genetic studies is the identification of simple sequence repeats (SSRs) or microsatellites. The mining of SSRs revealed that the M. chilensis genome has 548,360 SSR sequences, where 9% and 6% of the SSR loci were annotated for each mussel chromosome ( Figure S1). The most frequent SSR motif was the tetranucleotide, followed by the dinucleotides, accounting for 206,103 and 197,700 repeats, respectively. The entire SSR sequences accounted for 0.35% of the whole genome. The development of SSR markers offers a shortcut to assessing genetic diversity, which can potentially be applied in food authentication and genetic traceability for mussel species [79][80][81].

Protein-Coding Gene Prediction and Functional Annotation in the M. chilensis Genome
For the identification of protein-coding genes, de novo, homolog prediction, and RNAseq evidence were used as the training set ( Figure 2A). For homologous predictions, the protein sequences from Crassostrea gigas, Mytilus galloprovincialis, M. coruscus, and Dreissena polymorpha genomes were extracted using the respectively published references and aligned against the blue mussel genome using TBLASTN (E-value < 1 × 10 −5 ) ( Table 5). The gene sequence structure of each candidate gene and previously mentioned tools were used to predict protein-coding genes. Finally, a non-redundant reference gene set was generated using the EvidenceModeler (v.2.0) (EVM) and PASA2 tools (v.2.5.2) (Figure 2A). Taken together, 34,530 protein-coding genes were identified with a 6531 bp average transcript length, 1377 bp average CDS length, 4.92 average of exons per gene, and 1377 and 1316 average length of exons and introns, respectively (Table 6). Additionally, 516 tRNAs were predicted using tRNAscan-SE, and 143 rRNA genes were annotated using RNAmmer. For non-coding RNAs with putative regulatory roles, 1365 miRNAs and 43,011 long non-coding RNAs were identified and annotated within the M. chilensis genome ( Table 7). For functional annotation, the predicted proteins within the blue mussel genome were searched by homology against seven databases: Swiss-Prot, Nr, Nt, KEGG, eggnog, GO, and Pfam ( Figure 2A). Overall, 70.45%, 73.01%, 8.98%, 64.94%, 80.57%, 33.61%, and 96.33% of genes matched entries in these databases, respectively. A total of 34,530 genes (100%) were successfully annotated by gene function and conserved protein motifs ( Table 8). The genomic features annotated for the native blue mussel M. chilensis were displayed using a Circos plot [64]. Herein, this graphical representation shows the primary genomic features for the 14 chromosomes. Specifically, gene density, repeat density, GC content, rRNA localization, and ncRNAs were plotted. The transcriptome expression profiles for the mantle, gills, hemocytes, and digestive gland tissues were also displayed in connection with the syntenic blocks ( Figure 2B).

Comparative Genomics
Smooth-shelled blue mussels of the genus Mytilus represent a model group because of their cosmopolitan distribution, socioecological importance, and intriguing evolutionary history. This taxon provides new insights into the process of speciation and how hybridization and introgression can be one of the biggest threats to global mussel biodiversity [82]. A survey of single nucleotide polymorphisms (SNPs) on southern hemisphere blue mussels has provided a new layer for understanding their biology, taxonomy, and phylogeography [83,84]. However, SNP markers cannot be applied as a single tool to evidence chromosome rearrangement events during the Mytilus evolution. Here, whole-genome sequencing in smooth-shelled blue mussels and relative bivalve species is a priority for global mussel aquaculture, biosecurity, and conservation.
With the aim of exploring genomic rearrangements in Mytilus, the reported reference genomes for M. coruscus and M. chilensis were analyzed. Of the 34,530 predicted genes from the M. chilensis genome, 18,758 (54.32%) were found in syntenic collinear blocks after being compared with the M. coruscus genome ( Figure 3A). These syntenic blocks consisted of 671 alignments with a minimum of 5 genes per block. The number of alignments per chromosome ranged from 27 on chromosome 13 to 69 on chromosome 3. Chromosomes with higher genes in collinear blocks were chromosomes 1, 4, and 6, with 1227, 1091, and 1088 genes, respectively. Blocks with less than five genes or E-value < 1 × 10 −5 were discarded from this analysis. Most collinear blocks were located at the same pair of chromosomes between the two genomes. For example, M. chilensis Chr1 had only syntenic blocks with LG01 from M. coruscus in the same order. However, chromosomes 6 and 10 from M. chilensis had collinearity with chromosomes LG09 and LG02 in M. coruscus but were orientated as two inversed blocks per pair of chromosomes (red lines in Figure 3A and Figure S2). The genes in these alignments from inversed blocks were extracted, blasted, and gene ontology terms were identified. Enrichment analyses from GO terms were obtained from Chr10, and LG09 inversed blocks and Chr6 and LG02 pair of chromosomes ( Figure 3B,C). Most molecular function-enriched GO terms in the Chr10 and LG09 pair were associated with heat shock protein (HSP) binding. By contrast, in the Chr6 and LG02 pair, most of the enriched GO terms were associated with the mitochondria and biological processes related to autophagy or regulation of gene expression by epigenetic changes. Notably, chromosome rearrangements have been associated with adaptative genetic traits in marine organisms [85], where specific architectural proteins such as HSPs may have distinct roles in establishing 3D genome organization [86].  LG02, and Chr10 vs.

Comparative Analysis of Steamer-like Elements in Bivalvia
To explore the gene expansion of retrotransposon elements among representative species from Bivalvia, we primarily characterized the Steamer-like elements (SLEs) in M. chilensis using the approach described by Arriagada et al. [62]. The analysis evidenced that the genome of M. chilensis contains five copies of SLEs distributed in chromosomes 1, 6, 7, 10, and 11. The alignment showed that all SLE copies are flanked by two LTRs (5′ and 3′) containing the Gag-Pol ORFs and the domains annotated to protease, reverse transcriptase, RNAaseH, and integrase. Notably, an insertion composed of 12 nucleotides at posi- LG02, and Chr10 vs.
LG09. The x-axis indicates the negative log 10 (p-value).

Comparative Analysis of Steamer-like Elements in Bivalvia
To explore the gene expansion of retrotransposon elements among representative species from Bivalvia, we primarily characterized the Steamer-like elements (SLEs) in M. chilensis using the approach described by Arriagada et al. [62]. The analysis evidenced that the genome of M. chilensis contains five copies of SLEs distributed in chromosomes 1, 6, 7, 10, and 11. The alignment showed that all SLE copies are flanked by two LTRs (5 and 3 ) containing the Gag-Pol ORFs and the domains annotated to protease, reverse transcriptase, RNAaseH, and integrase. Notably, an insertion composed of 12 nucleotides at position 933ˆ934 was exclusively found in chromosomes 7 and 11. The translation for the inserted nucleotides suggests four amino acids, K, T, S, and H, in a positive orientation. However, the translation evaluated in the reading frame (−1) evidenced a methionine localized before the RNAaseH coding gene ( Figure 4A).   Furthermore, the phylogenetic analysis using publicly available reference genomes assembled at chromosome level for eleven bivalve species using maximum likelihood (ML) revealed a six-chromosome cluster composed of bivalves belonging to the families Veneridae, Solenidae, Pectinidae, Ostreidae, Pteriidae, and Mytilidae ( Figure 4B). The phylogenetic-reconstruction-rooted SLEs were found in three chromosomes (2, 4, and 17) from R. philippinarum. The other Veneridae member, M. mercenaria showed a cluster of four chromosomes (10, 12, 13, and 16) and related to two chromosomes of S. grandis (10 and 16). This last species formed a unique cluster composed of three chromosomes (8, 15, and 17), similar to P. maximus, with three chromosomes. Concerning the mussel and oyster genomes assembled at the chromosome level, the phylogenic analysis revealed two main clusters composed of species belonging to Ostreida and Mytilidae, where the first taxon was comprised of the Ostreidae and Pteriidae families. Herein, one cluster was rooted with three SLE sequences from C. virginica, C. gigas, and C. ariakensis located on chromosomes 9, 2, and 5, respectively. The second major cluster was composed of SLEs annotated in chromosomes from Ostreidae and Pteriidae, where C. virginica chromosomes were closely related to P. imbricata. The third cluster was observed containing three SLE sequences from C. virginica and C. ariakensis: chromosomes 1, 2, 8, and 1, 2, and 6, respectively ( Figure 4B). The analysis of the Mytilidae family revealed two primary clusters comprising SLEs lo-cated in chromosomes from M. edulis and D. polymorpha, and M. coruscus, M. chilensis, respectively ( Figure 4B). This last cluster grouped five chromosomes from M. coruscus (Chr. 1, 4, 5, 7, and 11), and two from M. edulis (Chr. 4 and 6). The Steamer-like sequence characterized for M. chilensis was also observed in this cluster. Finally, a detailed analysis of the three mussel species reported with genome assemblies at the chromosome level was conducted ( Figure 4C). Notably, a rooted cluster comprising chromosomes 7, 11, 9 for M. couscous, and 4 and 6 for M. edulis were closely related. Herein, two primary clusters of SLEs located in chromosomes from M. chilensis, M. edulis, and M. coruscus were observed. The analysis suggested that the SLEs identified on M. chilensis chromosomes are closely related to the SLEs annotated on chromosomes 9 and 4 in M. edulis; meanwhile, the SLEs located in chromosomes 1, 5, and 11 in M. coruscus were also identified in the same chromosome cluster. The second main cluster observed comprised exclusively SLEs annotated in D. polymorpha chromosomes, except the SLE copies identified in chromosomes 9 and 10 of M. edulis. Interestingly, the SLEs annotated in chromosome 9 from M. edulis are shared among the three primary clusters analyzed, suggesting putative translocation gene events in Mytilidae.
Overall, the phylogenic relationships of SLEs revealed that the reported bivalve genomes comprise between 3 and 6 loci. A lower number of SLEs was found in Solenidae, Pectinidae, and Veneridae, followed by Mytilidae. A higher number of SLE loci was observed in genomes belonging to the Ostreida order. As far as we know, the evolution of the bivalve chromosomes has mainly been studied using cytogenetic techniques combining molecular probes on candidate genes to detect genome rearrangements that drive the speciation process [87][88][89]. However, the availability of reference genomes assembled at the chromosome level opens new perspectives for exploring molecular evolution in several taxonomic orders through gene collinearity analysis. The study by Yang [28] highlighted putative chromosome rearrangements among the king scallop Pecten maximus, the blood clam Scapharca broughtonii, the hard-shelled mussel Mytilus coruscus, the pearl oyster Pinctada martensii, and the Pacific oyster Crassostrea gigas genomes. Notably, the chromosome synteny illustrated that large-scale rearrangements are common events between the scallop and oyster but scarce between the scallop and mussel genomes. The reported evidence suggested that almost all the chromosome rearrangements between the mussel and oyster genomes are different, implicating independent chromosome fusion events. The SLE loci identified in all the genomes analyzed in the current study suggest that SLEs are relatively conserved in chromosome position for some taxa. For instance, the SLE loci in Veneridae, Pectinidae, and Solenidae appear to be associated with chromosomes 10, 13, 12, and 16. This sharing characteristic can reflect common genetic events during the evolution of these taxonomical groups. Similarly, the Ostreidae and Mytilidae families share SLE loci annotated to chromosomes 1, 2, 8, and 10. The detailed analysis of SLEs in Mytilidae shows that the transposon identified in M. chilensis was shared between M. edulis and M. coruscus, where SLEs in D. polymorpha appear to be more phylogenetically distant than Mytilus species. Interestingly, the mutation identified on the SLEs localized in the M. chilensis genome (insertion of twelve nucleotides), specifically on chromosomes 7 and 11, was shared with the SLE annotated on chromosome 9 in M. edulis. This cumulative evidence reveals diverse chromosome rearrangements, reflecting a complex evolutionary history of bivalve chromosomes.

The Marine Environment of M. chilensis Populations
Temporal and spatial variability of sea surface temperature (SST) around Chiloé island and at Yaldad and Cochamó were analyzed over the past two decades. The oceanographic variability for the location studied was analyzed from remote sensing data and in situ measurements ( Figure 5A-D). Here, the daily time series of SST extracted from satellite-derived data for both sites evidenced high surface temperature variability between Yaldad and Cochamó, where this last location was constantly higher throughout the year ( Figure S3). Notably, the monthly medians computed from the SST time series showed that the main differences were observed during the austral summer from December to March. During the winter, the oceanographic variability was less pronounced, showing temperatures between 13 • C and 10 • C from April to July ( Figure 5C).
Temporal and spatial variability of sea surface temperature (SST) around Chiloé island and at Yaldad and Cochamó were analyzed over the past two decades. The oceanographic variability for the location studied was analyzed from remote sensing data and in situ measurements ( Figure 5A-D). Here, the daily time series of SST extracted from satellite-derived data for both sites evidenced high surface temperature variability between Yaldad and Cochamó, where this last location was constantly higher throughout the year ( Figure S3). Notably, the monthly medians computed from the SST time series showed that the main differences were observed during the austral summer from December to March. During the winter, the oceanographic variability was less pronounced, showing temperatures between 13 °C and 10 °C from April to July ( Figure 5C).  Furthermore, in situ data were collected from June 2017 to May 2018, exhibiting significant differences in both locations for temperature and salinity between 0 and 20 m of depth. Interestingly, the oceanographic survey revealed a pronounced vertical stratification with higher temperatures and lower salinity in Cochamó compared with Yaldad ( Figure 5D). These observations support the idea that two oceanographically different zones exist in the inner sea of Chiloé Island. In this northern area, mussels from Cochamó and Yaldad were sampled for the current study. Taken together, we can hypothesize that the mussels inhabiting Cochamó are significantly more exposed to environmental stress than the Yaldad mussel population. To date, there are few studies exploring how population genetic variation is related to, or caused by, the marine environmental variation in mussel populations. Notably, a study conducted by Wenne et al. [90] examined the genetic differentiation of native populations of M. galloprovincialis throughout its entire geographic range in the Mediterranean Sea, the Black Sea, and the Sea of Azov using 53 SNP loci. The results indicated that 7 of the 13 environmental variables explained significant variation in population-specific SNP locus allele frequencies. These seven variables explained 75% of the variation in the SNP dataset, suggesting that a complex mix of environmental variables contributes to the genetic variation in M. galloprovincialis populations in the Mediterranean Sea.

Whole-Genome Transcript Expression Analysis in Two M. chilensis Populations
The transcriptome profiling among mussels collected during the austral summer in 2019 from Yaldad and Cochamó evidenced three primary transcriptional clusters. Herein, gene cluster 1 was highly expressed in the gills of mussels exposed to the Yaldad marine conditions; meanwhile, gene clusters 2 and 3 were highly expressed in individuals collected in Cochamó or mussels exposed to estuarine conditions ( Figure 6A). Notably, the RNAseq from individuals collected as Cochamo1 (replicate) showed a highly expressed gene cluster, indicating a wide transcriptome variation among mussels from this population. The RNA-seq analysis was performed with the mRNA sequences annotated on the M. chilensis genome. Herein, it is essential to note that in mussel species, specifically in M. galloprovincialis, the phenomenon of presence-absence variation (PAV) has been described. This fact means that PAVs can bias the analyses of transcriptome profiles in the studied mussel populations. We previously conducted a de novo assembling for the RNA-data sets sequenced from Yaldad and Cochamó populations. The results showed that the number of genes with expression values >1 (total gene reads) did not show statistical differences between both mussel populations. For instance, Yaldad and Cochamó mussels showed 25,086 ± 215 and 25,344 ± 212 (three replicates per population), respectively, of expressed genes in gill tissue (p-value = 0.98). Collectively, between 72.6% and 73.3% of the annotated genes in the M. chilensis genome were transcriptionally active in gills independently of the population analyzed.
Genes 2023, 14, x FOR PEER REVIEW  19 of 29 processes, defense response, cell differentiation, and anatomical structure development ( Figure 6E). Clusters 2 and 3 were less enriched, revealing significant GO terms for transmembrane transport, reproductive processes, protein-containing complex assembly, microtube-based movement, cytoskeleton organization, chromatin organization, and metabolic process ( Figure 6E). The cluster gene expression analysis was used to identify genetic polymorphisms annotated in differentially expressed genes (DEGs) between the Yaldad and Cochamó mussel populations. The DEGs were evaluated by cluster transcriptome analysis displayed using a Circos plot to visualize specific loci where DEGs were highly transcribed. The foldchange values calculated showed high transcription levels in clusters 1 and 2 through all chromosomes scanned (see red dots in Figure 7A). Congruently with the previous RNAseq results in this study, the highest fold-change values were observed in DEGs annotated in cluster 1 (Yaldad population). By contrast, cluster 3 showed a small number of DEGs with high fold-change values. Notably, the physical mapping of DEGs on chromosomes evidenced specific transcriptome patterns, revealing genes differentially expressed through the mussel genome exposed to the marine environment. The synteny analysis for DEGs showed a marked pattern among chromosomes 5, 7, and 12 for cluster 1; mean- The evaluation of differentially expressed genes (DEGs) showed that the main factor of differences in the number of DEGs was the population rather than the replicates assessed ( Figure 6B). The proportion of DEGs evaluated among the gene clusters revealed that cluster 1, highly expressed in Yaldad, accounted the 78.85% of the total DEGs analyzed. Clusters 2 and 3 are primarily characterized by high transcription values in the Cochamó population, evidenced by 7.32% and 13.82% of DEGs, respectively. The total number of DEGs analyzed was 1570 ( Figure 6C). Notably, the fold-change values estimated among the replicates and populations revealed high values in gene transcriptional cluster 1, compared with clusters 2 and 3 where the fold-change values were significantly lower ( Figure 6D). The functional analysis showed that cluster 1 was enriched by GO terms related to protein modification processes, programmed cell death, immune system processes, defense response, cell differentiation, and anatomical structure development ( Figure 6E). Clusters 2 and 3 were less enriched, revealing significant GO terms for transmembrane transport, reproductive processes, protein-containing complex assembly, microtube-based movement, cytoskeleton organization, chromatin organization, and metabolic process ( Figure 6E).
The cluster gene expression analysis was used to identify genetic polymorphisms annotated in differentially expressed genes (DEGs) between the Yaldad and Cochamó mussel populations. The DEGs were evaluated by cluster transcriptome analysis displayed using a Circos plot to visualize specific loci where DEGs were highly transcribed. The fold-change values calculated showed high transcription levels in clusters 1 and 2 through all chromosomes scanned (see red dots in Figure 7A). Congruently with the previous RNAseq results in this study, the highest fold-change values were observed in DEGs annotated in cluster 1 (Yaldad population). By contrast, cluster 3 showed a small number of DEGs with high fold-change values. Notably, the physical mapping of DEGs on chromosomes evidenced specific transcriptome patterns, revealing genes differentially expressed through the mussel genome exposed to the marine environment. The synteny analysis for DEGs showed a marked pattern among chromosomes 5, 7, and 12 for cluster 1; meanwhile, the synteny observed for the DEGs annotated in clusters 2 and 3 revealed a wide distribution along the M. chilensis genome ( Figure 7A). Interestingly, the analysis carried out to detect macro-genome mutation in gene families between the Yaldad and Cochamó populations evidenced a similar number of dispersed genes, suggesting that those might arise from transposition. Tandem or repeatedly duplicated genes were observed with a low proportion in cluster 3 (Cochamó); meanwhile, the proximal genes showed a similar proportion to cluster 1 (Yaldad). These results might suggest small-scale transposition or duplication/insertion events. An interesting finding was observed for whole-genome duplication (WGD). The primary proportion was evidenced in cluster 3 (Cochamó), compared with clusters 1 and 2 from the Yaldad population ( Figure 7B). Furthermore, the bioinformatic analysis conducted for detecting amino acid changes (AAC) in DEGs showed that 38% of non-synonymous AAC were identified in mussels collected from Yaldad. By contrast, the main proportion of synonymous AAC was detected in mussels exposed to Cochamó's estuarine conditions ( Figure 7C). Notably, the analysis performed for DEGs annotated in cluster 2 did not show non-synonymous and synonymous AAC in mussels collected from Yaldad. Finally, evaluating the zygosity proportion estimated for each mussel population evidenced an inverse pattern between both populations. The Yaldad cluster was higher in the homozygous proportion than Cochamó, where heterozygous AAC were detected in a higher proportion ( Figure 7D).
To explore the transcriptome signatures between the Yaldad and Cochamó mussel populations, we applied the genome chromosome expression (CGE) approach to test differences among tissues and individuals through the M. chilensis genome. The CGE analysis revealed high differences among chromosome regions, where the gill tissue was more modulated than the mantle tissue ( Figure 8A). Interestingly, there are some levels of congruence among the CGE annotated for both mussel populations. We conducted a gene ontology enrichment analysis using this finding from genes identified by CGE analysis. The results evidenced that gill transcriptomes displayed functional processes associated with transmembrane transport, protein catabolic, nervous system, and metal homeostasis. Notably, immune system process GO terms were highly enriched in gills. Moreover, the chromosome region differentially expressed in mantle tissue revealed that the reproductive process, protein modification process, cell differentiation, anatomical structure development, and gene silencing by RNA were mainly annotated ( Figure 8B). Taken together, the results reported in this study are highly congruent with the previous study conducted by Yévenes et al. [38], through the transcriptome responses of M. chilensis collected in ecologically different farm-impacted seedbeds. pared with clusters 1 and 2 from the Yaldad population ( Figure 7B). Furthermore, the bioinformatic analysis conducted for detecting amino acid changes (AAC) in DEGs showed that 38% of non-synonymous AAC were identified in mussels collected from Yaldad. By contrast, the main proportion of synonymous AAC was detected in mussels exposed to Cochamó's estuarine conditions ( Figure 7C). Notably, the analysis performed for DEGs annotated in cluster 2 did not show non-synonymous and synonymous AAC in mussels collected from Yaldad. Finally, evaluating the zygosity proportion estimated for each mussel population evidenced an inverse pattern between both populations. The Yaldad cluster was higher in the homozygous proportion than Cochamó, where heterozygous AAC were detected in a higher proportion ( Figure 7D).  Singleton means that the gene is single-copy, which should not be the type of members of gene families. Dispersed means that the gene might arise from transposition. Tandem means that genes were repeatedly duplicated. Proximal means that the gene might arise from small-scale transposition or from tandem duplication and insertion of some other genes. Wholegenome duplication (WGD) means the gene might arise from a chromosome duplication region. The analysis was carried out using MCScanX. (C) Amino acid change proportions (%) between the Yaldad and Cochamó populations. The non-synonymous and synonymous were annotated for the DEGs selected for each mussel population according to the cluster analysis. (D) Zygosity proportion (%) estimated for each mussel population. Cluster 1 (Yaldad) is represented by blue bars, and clusters 2 and 3 (Cochamó) are displayed in brown bars.
The cumulative findings of this study suggest that the immune system was primarily modulated between mussels exposed to the Yaldad and Cochamó environmental conditions. With the aim of exploring the transcription profiling of immune-related genes, we selected two KEGG pathways annotated in the M. chilensis genome (Figure 9). Herein, Toll-like receptor signaling pathway and apoptosis were analyzed in terms of transcription activity and single nucleotide variation (SNV) between mussel populations. Notably, a nonsynonymous SNV was detected on the TLR2 gene (28T>G) in individuals collected from the Yaldad population. The translation evidenced an amino acid change from phenylalanine to valine at position 10 in the ORF (Phe10Val) ( Figure 9A). The analysis also evidenced SNV on genes such as AKT and TAB1, where no amino acid changes were detected. The transcriptome profiling for the TLR pathway evidenced a high modulation of genes such as TLR3, AKT, TRAFF6, FADD, IRAK4, and RAC1 in mussels collected from Yaldad. Interestingly, mitogen-activated protein kinases (MAP2K and MAPK1) and c-Jun N-terminal kinase (JNK) were differentially expressed, suggesting putative roles related to stress signaling pathways ( Figure 9B). Furthermore, the apoptosis pathway revealed two SNV localized in eukaryotic initiation factor 2 α (EIF2α) and inhibitor of apoptosis (IAP) in mussels sampled from the Yaldad and Cochamó populations, respectively ( Figure 9C). The 2613delG in the EIF2α gene produces a frameshift at the Thr872; meanwhile, the 968_970delCTC localized in the IAP gene produces a deletion of proline at position 323. The transcriptome profiling of apoptosis-related genes showed a conspicuous differentiation between gills and mantle tissue, where three primary gene expression clusters were identified ( Figure 9D).
To explore the transcriptome signatures between the Yaldad and Cochamó mussel populations, we applied the genome chromosome expression (CGE) approach to test differences among tissues and individuals through the M. chilensis genome. The CGE analysis revealed high differences among chromosome regions, where the gill tissue was more modulated than the mantle tissue ( Figure 8A). Interestingly, there are some levels of congruence among the CGE annotated for both mussel populations. We conducted a gene ontology enrichment analysis using this finding from genes identified by CGE analysis. The results evidenced that gill transcriptomes displayed functional processes associated with transmembrane transport, protein catabolic, nervous system, and metal homeostasis. Notably, immune system process GO terms were highly enriched in gills. Moreover, the chromosome region differentially expressed in mantle tissue revealed that the reproductive process, protein modification process, cell differentiation, anatomical structure development, and gene silencing by RNA were mainly annotated ( Figure 8B). Taken together, the results reported in this study are highly congruent with the previous study conducted by Yévenes et al. [38], through the transcriptome responses of M. chilensis collected in ecologically different farm-impacted seedbeds. The cumulative findings of this study suggest that the immune system was primarily modulated between mussels exposed to the Yaldad and Cochamó environmental conditions. With the aim of exploring the transcription profiling of immune-related genes, we selected two KEGG pathways annotated in the M. chilensis genome (Figure 9). Herein, Toll-like receptor signaling pathway and apoptosis were analyzed in terms of transcription activity and single nucleotide variation (SNV) between mussel populations. Notably, a non-synonymous SNV was detected on the TLR2 gene (28T>G) in individuals collected Notably, genes such as P53, ERK, TP53, PARP2, and JNK were highly expressed in gill tissue. The gene expression analysis in mantle tissue evidenced high transcriptional activity in genes related to the intrinsic (mitochondria-mediated) pathway, such as the B- Notably, genes such as P53, ERK, TP53, PARP2, and JNK were highly expressed in gill tissue. The gene expression analysis in mantle tissue evidenced high transcriptional activity in genes related to the intrinsic (mitochondria-mediated) pathway, such as the B-cell lymphoma (BCL) gene and the second mitochondria-derived activator of caspase (DIABLO) gene. Concerning the mitogenome of M. chilensis, it was previously reported by Gaitán-Espitia et al. [34] andŚmietanka and Burzyński [91], evidencing a genome size of 16,748bp and structurally identical to the northern hemisphere M. edulis and M. galloprovincialis mitogenomes. Furthermore, the putative adaptive contribution of the mitochondrial genes was recently reported by Yevenes et al. [37]. The RNA-Seq analysis detected differences in the number of upregulated mitogens between individuals from Cochamó and Yaldad, some being tissue-specific (ND4L and COX2). Several monomorphic location-specific mitochondrial genetic variants were detected in samples from Cochamó and Yaldad, representing standing genetic variability to optimize mitochondrial functioning under local habitats. Overall, these mitochondrial transcriptomic differences reflect the impact of environmental conditions on the mitochondrial genome functioning and offer new markers to assess the effects of habitat translocations on mussel fitness, a routine industry practice. Likewise, these mitochondrial markers should help monitor and maintain population differences in this keystone and heavily exploited native species.
Further functional studies will be conducted to validate the association between single nucleotide polymorphisms and the fitness traits observed or how the translocation process associated with aquaculture activity can evolve with the loss of locally adapted alleles. Interestingly, a recent transplant experiment reported by Jahnsen-Guzmán et al. [92] demonstrated that M. chilensis individuals are adapted to the subtidal environment (4 m depth), as they exhibit significantly higher fitness (growth and calcification rates) than those transferred to the intertidal environment (1 m depth), which showed increased metabolic stress. Herein, the mussel lives in extreme environmental variability, where their ability to cope with perturbations, and build plasticity and adaptive responses, seems based on the genome architecture.

Conclusions
This study reports the first chromosome-level genome assembly of the native blue mussel Mytilus chilensis. This genomic resource was used to identify genome signatures related to the phenotypic plasticity in the mussel population inhabiting contrasting marine environments. Collectively, the putative mutations associated with immune-and metabolism-related genes suggest molecular regulatory mechanisms to reduce the number of genes and their transcriptional activity under stress conditions. This evolutionary strategy can suggest that the expression of those genes has evolved a degree of "frontloading" that potentially pre-adapts the mussel populations to frequent heat and salinity stress, contributing to their physiological tolerance and fitness. We believe that the generated genomic resource must be instrumental for future research on population genomics informing management and sustainable strategies for the Chilean mussel aquaculture.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes14040876/s1, Figure S1: Distribution of Simple Sequence Repeats (SSR) identified by SSR Finder in the M. chilensis genome; Figure S2: Dot-blot analysis between M. chilensis and M. coruscus chromosomes, and macro-syntenic relationships; Figure S3: Temporal and spatial variability of Sea Surface Temperature (SST) around Chiloé island, and at sites Yaldad and Cochamo, over the past two decades. Funding: This study was funded by FONDAP #1522A0004, FONDECYT #1210852, and #1130852, granted by ANID-Chile.

Institutional Review Board Statement:
The Committee of Ethics, Bioethics, and Biosafety, University of Concepción, Chile, approved this project (CEB324-2022, September 2021).

Informed Consent Statement: Not applicable.
Data Availability Statement: The Mytilus chilensis whole-genome sequencing data supporting this study's findings are available from NCBI under BioProject PRJNA861856. The sequencing data supporting this study's findings are available in SRA at SRR20966976, SRR20593343, and SRP261955. The benefits from this study accrue from sharing our data and results on public databases as described above. The assembled genome and the genome annotation results were deposited in the Figshare database [93].

Conflicts of Interest:
The authors declare no conflict of interest.