The diversity and evolution of retroviruses: Perspectives from viral “fossils”

Retroviruses exclusively infect vertebrates, causing a variety of diseases. The replication of retroviruses requires reverse transcription and integration into host genomes. When infecting germline cells, retroviruses become inherited vertically, forming endogenous retroviruses (ERVs). ERVs document past viral infections, providing molecular fossils for studying the evolutionary history of retroviruses. In this review, we summarize the recent advances in understanding the diversity and evolution of retroviruses from the perspectives of viral fossils, and discuss the effects of ERVs on the evolution of host biology.

Introduction
Retroviruses (the Retroviridae family) exclusively infect vertebrates and cause a variety of diseases, including acquired immune deficiency syndrome (AIDS) and cancers (Goff, 2007;Coffin et al., 1997;Telesnitsky, 2010). Unique among RNA viruses, retrovirus replication requires reverse transcription and integration into host genomes. Retroviruses usually infect somatic cells. But when retroviruses infect germline cells, the integrated retroviruses are passed on vertically to the next generation, forming endogenous retroviruses (ERVs). ERVs were first discovered in 1960s (Weiss, 2013). In the pre-genomic era, many ERVs had been identified by Southern blot and/or PCR approaches (e.g. Repaske et al., 1985;Tristem et al., 1995;Tristem et al., 1996). Thanks to the recent development of next-generation sequencing techniques, the genomes of many vertebrates have been sequenced, providing important resources for studying the extent and distribution of ERVs with the vertebrate genomes. ERVs are widespread and highly abundant in vertebrate genomes; for example, ERVs account for~8% of the human genome (Hayward et al. 2013(Hayward et al. , 2015Xu et al., 2018).
It is notoriously hard to study the deep evolution of viruses over geographic time scale due to high mutation rates of viruses (Duffy et al., 2008;Gago et al., 2009) and no physical fossils preserved. Retroviruses are no exception. Reverse transcriptase lacks the ability of proofreading, resulting in high mutation rates in retroviruses. However, ERVs evolve with their hosts as genomic loci, and thus have much lower evolutionary rates (Duffy et al., 2008;Gago et al., 2009). ERVs provide "molecular fossils" for studying the deep history of retroviruses and the ancient evolution of host-retrovirus interaction. More generally, many viruses, both RNA and DNA viruses, can integrate into host genomes, forming endogenous viral elements (EVEs) (Holmes, 2011;Feschotte and Gilbert, 2012). The discovery of EVEs in a wide range of eukaryotes lays the foundation of an emerging field, namely Paleovirology. Paleovirology is the study of the ancient evolution of viruses through analyzing endogenous viral elements or evolutionary signatures of host-virus conflicts left in the host genomes (Emerman and Malik, 2010;Katzourakis, 2013). In this review, we introduce methods to analyze ERVs, summarize the recent advances in understanding the diversity and evolution of retroviruses from the perspectives of viral fossils, and discuss the effects of ERVs on the evolution of host biology.

Retroviruses and endogenous retroviruses
Retroviruses are enveloped viruses with RNA genomes~7-10 kb in length. All the retroviruses encode three common genes, namely gag, pol, and env. The gag gene encodes matrix (MA), capsid (CA), and nucleocapsid (NC) proteins, the pol gene encodes protease (PR), reverse transcriptase (RT), RNase H (RH), and integrase (IN) proteins, and the env gene encodes transmembrane (TM) and surface (SU) glycoproteins. Some complex retroviruses can also encode several accessory genes (Goff, 2007;Coffin et al., 1997;Telesnitsky, 2010). Unlike other RNA viruses, the replication of retroviruses requires reverse transcription of viral RNA genomes into double stranded DNA, and integration of the newly synthesized DNA into host genomes. Retrovirus infection typically occurs in the somatic cells of hosts. On occasion, retroviruses infect host germline cells, and the integrated retroviruses flanked by two long terminal repeats (LTRs) become vertically inherited to next generation as part of the host chromosomes, forming endogenous retroviruses (ERVs) (Stoye, 2012;Johnson, 2015) (Fig. 1).
Like new mutations, the fate of novel retrovirus insertions within the host population is governed by population genetics processes. As new mutations, most of new ERV insertions are deleterious and removed from host population by purifying selection (Wang and Han, 2020). The integration of exogenous retroviruses occurs preferentially near or within host genes, but ERVs are found mainly outside host genes, suggesting that a majority of ERVs have been removed by negative selection (Brady et al., 2009;Emerman and Malik, 2010). The fates of the ERVs that have little effect on host fitness and are selectively neutral are governed by random genetic drift. The probability that a newly retroviral insertion is fixed in the host population is 1/2N, where N is the effective population size of hosts (Hartl and Clark, 2006). These neutral ERVs might accumulate disruptive mutations (premature stop codons, frameshift mutations, and insertions or deletions), which eventually lead to the degradation of these ERVs. In extreme but frequent cases, recombination occurs between 5 0and 3-LTRs, generating a solo-LTR. In the rare cases that an ERV insertion is beneficial, it will be fixed in the host population by positive selection (Fig. 2). In the meantime, ERVs retain the ability to proliferate within the host genomes by retrotransposition in cis (ERVs use their own proteins for amplification), complementation in trans (ERVs hitchhike proteins produced by other ERVs or transposable elements within the same cell), or reinfection (ERVs produce virus particles to infect other germline cells), increasing their copy numbers within the host genomes (Belshaw et al., 2005;Bannert and Kurth, 2006;Magiorkinis et al., 2012) (Fig. 3). Together, these processes generate complex insertion polymorphism pattern of ERVs within the host populations.
ERV insertion polymorphisms have been found in the populations of many vertebrates (Zhang et al., 2008;Turner et al., 2001;Roca et al., 2004;Polani et al., 2010;Ngo et al., 2019). An ongoing fixation of endogenous koala retrovirus (KoRV) had been observed in the wild koala (Phascolarctos cinereus) population: all the individuals in the north of Australia had endogenous KoRV, while KoRV was not detected in some individuals in the south of Australia (Tarlinton et al., 2006). Consistently, southern koala individuals carried fewer KoRV loci than northern ones (Ishida et al., 2015). The KoRV endogenization started before the late 1800s, but no more than 22,200-49,990 years ago ( Avila-Arcos et al., 2013;Ishida et al., 2015). Tens of HERV-K elements were found to be polymorphic in the human population (Wildschutte et al., 2016), and several HERV-K elements were absent from modern human genomes and present in archaic (Denisovan and Neanderthal) hominids (Lee et al., 2014), suggesting the recent activity of the HERV-K elements in the human population.

Mining ERVs
ERVs often accumulate a large number of mutations, making them difficult to identify. Several main strategies have been commonly used to mine ERVs within the vertebrate genomes. (i) ERVs can be identified using similarity search with representative retroviruses as queries. However, representative retroviruses might not cover the major diversity of retroviruses. (ii) Because ERVs are flanked by two LTRs, ERVs can be identified through first detecting candidate LTRs and then detecting and reconstructing retroviral proteins, as exemplified by RetroTector (Sperber et al., 2007). However, ERVs frequently accumulate large-scale deletions and become fragmented, and thus not all the ERVs are flanked by two LTRs. (iii) ERVs can also be identified through first searching homologs of retrovirus hallmark RT proteins. Because retrovirus RT proteins share detectable similarity with retrotransposon RT proteins, phylogenetic analyses are then performed to identify ERVs that cluster with representative retroviruses (Xu et al., 2018). Once again, fragmented ERVs that lack RT proteins cannot be identified using this approach. In general, the approach based on identifying RT homologs is more sensitive than the one based on detecting LTRs (Xu et al., 2018). Nevertheless, old ERVs are often more difficult to detect than young ERVs.

Classification and taxonomy of retroviruses
Based on the latest release of International Committee on the Taxonomy of Viruses (ICTV), the Retroviridae family is classified into two subfamilies, namely Orthoretrovirinae and Spumaretrovirinae. The Orthoretrovirinae subfamily includes six genera (Alpharetrovirus, Betaretrovirus, Gammaretrovirus, Deltaretrovirus, Epsilonretrovirus, and Lentivirus), and the Spumaretrovirinae subfamily (also known as foamy viruses) includes five genera (Bovispumavirus, Equispumavirus, Felispumavirus, Prosimiispumavirus, Simiispumavirus) (Hendrickson et al., 2018) (Fig. 4A). However, the ICTV classification only accounts for exogenous retroviruses. ERVs have been traditionally classified into three classes based on their similarity and phylogenetic relationship with exogenous retroviruses. Class I ERVs cluster with Gammaretrovirus and Epsilonretrovirus, Class II ERVs cluster with Alpharetrovirus, Betaretrovirus, Deltaretrovirus, and Lentivirus, and Class III ERVs cluster with foamy viruses and ERV-L (Gifford and Tristem, 2003;Gifford et al., 2018) (Fig. 4A). For historical reasons, human ERVs (HERVs) had been grouped into families primarily based on the specificity of tRNA primer binding sites (PBS) (Bannert and Kurth, 2006;Jern and Coffin, 2008). However, this classification system had problems: ERVs that belong to the same phylogenetic group might have different types of PBS, while unrelated ERVs might share the same type of PBS, and mutations and deletions were present in PBS regions (Bannert and Kurth, 2006;Jern and Coffin, 2008). Therefore, the classification of exogenous and endogenous retroviruses has not been well incorporated. Most of ERVs cannot be assigned into certain genus. It is also challenging to  classify and name ERVs, when ERVs arise from recent endogenization events are nested with the diversity of exogenous retroviruses (Blomberg et al., 2009).
It is logic to include both exogenous and endogenous retroviruses in a taxonomy of retroviruses (Blomberg et al., 2009;Gifford et al., 2018). Based on phylogenetic analyses of Pol proteins from exogenous and endogenous retroviruses sampled across vertebrates, retroviruses have been classified into five major clades, designated Jin, Mu, Shui, Huo, and Tu (Xu et al., 2018). Clades Jin and Mu include viruses related to gammaretroviruses and epsilonretroviruses, respectively, and they include class I ERVs. Clade Shui includes viruses related to alpha-, beta-, delta-retroviruses, lentiviruses, and class II ERVs. Clade Huo is related to snakehead retrovirus, ERV-L, and class III ERVs. Clade Tu is related to foamy viruses (Fig. 4A). However, RH and IN proteins were acquired multiple times in retroviruses, making phylogenetic analyses based on Pol proteins should be taken with cautions (Wang and Han, 2021a). On the other hand, based on phylogenetic analyses of the hallmark RT proteins, the positions of some retroviruses are unstable, and the support values for some groups are low. The recent discovery of lokiretroviruses makes the situation more complicated (Wang and Han, 2021a). Lokiretroviruses are sister to all the known retroviruses, and might represent a novel subfamily of retroviruses. We anticipate a classification system that fully integrates endogenous and exogenous retroviruses and follows phylogeny will be developed in the future.

Dating ERV invasion
The time for ERV invasion into host genomes can help establish time scales for retrovirus evolution. Several methods have been developed to estimate the integration time of ERVs. (i) The most commonly used approach is based on the divergence of 5 0 -and 3 0 -LTRs of an ERV. When integrated into host genomes, 5 0 -and 3 0 -LTRs of an ERV are identical at sequence level (Coffin et al., 1997). Then, they accumulate mutations independently as distinct genomic loci. Therefore, the sequence  (Xu et al., 2018;Zheng et al., 2021;Hayward et al., 2015;Han and Worobey, 2012a, 2012b, 2012c, 2015Chen et al. 2019Chen et al. , 2021aChen et al. , 2021bWang and Han, 2021a). divergence between 5 0 -and 3 0 -LTRs are proportional to the integration time (Johnson and Coffin, 1999). The integration time of an ERV, T ¼ D/2μ, where D indicates the genetic distance between 5 0 -and 3 0 -LTR, and μ indicates host neutral evolutionary rate (Fig. 5A). However, the approach comes with at least two caveats: 5 0 -and 3 0 -LTRs may evolve at a rate different from host neutral evolutionary rate (Martins and Villesen, 2011); gene conversion frequently occurs in LTR regions, which decreases their genetic distance (Johnson and Coffin, 1999). Moreover, given ERVs frequently accumulate large-scale deletion and become fragmented, not all the ERVs are flanked by two LTRs and thus their ages cannot be estimated using the approach based on the divergence between 5 0 -and 3 0 -LTRs. (ii) The age of an ERV can also be estimated based on its orthologous insertions within closely related species. The possibility that two ERVs integrated independently at an orthologous position of closely related species is extremely low (Stoye, 2001). Therefore, ERV insertion should occur in the last common ancestor of species who share orthologous ERVs, and the divergence time of their last common ancestor provides a minimum time for the ERV insertion. This approach might sometimes underestimate the invasion time, if an orthologous insertion is lost in certain outgroup species (Fig. 5B). (iii) ERV insertions might undergo segmental duplication. ERVs loci generated by segmental duplication can be identified by ERVs flanking with highly similar host sequences. Calculating the divergence time of segmental duplicates provides a minimum estimate of the invasion time (Katzourakis et al., 2007) (Fig. 5C). (iv) Rarely, the invasion time of an ERV can be estimated based on the occurrence of stop codons within it. Under neutral evolution, stop codons are generated at a rate of~1/2310 per codon per million years in primate coding sequences (Belyi et al., 2010). However, it is unclear whether ERVs evolve neutrally after integration, and the reliability of this method should be taken with caution (Fig. 5D).

Origin, evolution, and diversity of retroviruses
Systemic analyses of ERVs in vertebrates provide important insights into the origin, evolution, and diversity of retroviruses. Large-scale phylogenomic analyses show that ERVs are present in all the genomes of jawed vertebrates studied to date (Hayward et al. 2013(Hayward et al. , 2015Xu et al., 2018;Zheng et al., 2021). Within jawless vertebrates, ERVs (later classified to endogenous lokiretroviruses) were also identified in sea lamprey (Petromyzon marinus) and Arctic lamprey (Lethenteron camtschaticum) (Hayward et al., 2015;Xu et al., 2018;Wang and Han, 2021a). ERVs have not been identified in any species outside vertebrates (Xu et al., 2018;Wang and Han, 2021a). The distribution of retroviruses suggests that retroviruses may have originated during the early evolution of vertebrates >450 million years ago (Xu et al., 2018). However, the distribution might also be explained by a later origin and spread across vertebrate through frequent host switching. Foamy viruses, a subfamily of retroviruses, have been found to mainly co-diverge with their jawed vertebrate hosts for >450 million years (Han and Worobey, 2012a;Aiewsakun and Katzourakis, 2017;Xu et al., 2018;Wei et al., 2019;Chen et al., 2021a). Moreover, retroviruses of tetrapods (amphibians, reptiles, birds, and mammals) nest within the diversity of fish retroviruses, and fish retrovirus occupy basal positions of all major retrovirus clades. These findings reveal an ancient aquatic origin of retroviruses >450 million years ago (Xu et al., 2018).
The discoveries of ERVs also help establish or refine the time scales of retrovirus evolution. Molecular clock analyses of exogenous lentiviruses and delta-retroviruses indicate their evolutionary time scales in the unit of hundred years, and both have been thought to be young viral groups (Van Dooren et al., 2001;Wertheim and Worobey, 2009). Endogenous lentiviruses have been identified in lagomorphs, primates, carnivores, and Dermoptera (Katzourakis et al., 2007;Keckesova et al., 2009;Gifford et al., 2008;Gilbert et al., 2009;Worobey, 2012b, 2015;Cui and Holmes, 2012). These endogenous lentiviruses invaded host genomes from 1.9 to 3.8 million years ago (Mya) (in Malagasy lemurs) to 21-40 Mya (in Colugo) Worobey, 2012b, 2015;Cui and Holmes, 2012;Hron et al., 2016). Endogenous delta-retroviruses have been identified in bats, carnivores, cetaceans, and insectivores (Farka sov a et al., 2017;Hron et al., 2018;Zheng et al., 2021). A delta-retrovirus was estimated to integrate into carnivore genomes >24.6 Mya (Hron et al., 2018). These findings strongly support that both lentiviruses and delta-retroviruses are not young viral groups, and they have been circulating in mammals for millions of years.

Biological function of endogenous retroviruses
Upon infecting host germ cells, ERVs are inherited in the Mendelian way. ERVs have various effects on their hosts: deleterious, neutral, and beneficial (Fig. 2). If ERVs are inserted into host genes, they are likely to interrupt the normal open reading frames and function of the host genes, which are generally detrimental. Indeed, ERVs are mainly found outside of the host genes (Brady et al., 2009;Emerman and Malik, 2010). ERVs have been proposed to be responsible for a wide range of diseases, such as cancer, autoimmune diseases, and neurodegenerative diseases, but these hypotheses remain controversial and unproven (Voisset et al., 2008;Stoye, 2012;Bhardwaj and Coffin, 2014;Küry et al., 2018). On the other hand, the transcripts or proteins of HERVs might represent biomedicine markers for the diagnosis of those diseases (Song et al., 2021). The deleterious ERV insertions will be removed from the host population by purifying selection. Most of the remaining ERVs are neutral and will become degraded due to the accumulation of disruptive mutations, and insertions or deletions. Beneficial ERVs will be recruited by hosts for host biological functions, a process known as co-option or exaptation Han, 2020, 2021b) (Fig. 2).
Broadly, the functionality of ERVs can be divided into three categories, namely physiological, immunological, and regulatory functions (Frank and Feschotte, 2017;Feschotte, 2008;Chuong et al. 2016Chuong et al. , 2017. (i) Physiological roles. Some ERV Env proteins, termed as syncytins, were repurposed to function in placenta formation and development (Mi et al., 2000;Blaise et al., 2003). Synctins originated tens of times in mammals and a viviparous lizard, which might be pivotal in the emergence of mammals (Dupressoir et al., 2005;Lavialle et al., 2013;Cornelis et al., 2017). Moreover, Synctins might also contribute to myoblast fusion and muscle sexual dimorphism in mice (Redelsperger et al., 2016). (ii) Immunological roles. Some ERV-derived proteins can protect the hosts from exogenous retrovirus infection. Expression of the ERV env gene is known to inhibit exogenous retrovirus infection by saturating receptors recognized by exogenous retroviruses (Aswad and Katzourakis, 2012;Kozak et al., 1984;McDougall et al., 1994;Palmarini et al., 2004;Tikhonenko and Lomovskaya, 1990). Hosts also employ the ERV-derived proteins to block the replication of exogenous retroviruses. Some endogenous jaagsiekte sheep retroviruses (JSRVs) express defective Gag proteins, which prevent virion release and interfere with exogenous JSRV (Arnaud et al., 2007;Palmarini et al., 2004;Mura et al., 2004). Fv1 gene was derived from the gag gene of a murine ERV-L that was inserted into the genomes of species in Muroidae~45 million years ago, and serves as the restriction factor to protect the hosts against infections of various retroviruses (Boso et al., 2018;Young et al., 2018;Best et al., 1996;Yap et al., 2014;Yan et al., 2009). Moreover, the level of ERV expression controls host's respond to the microbiota in homeostatic or inflammatory manner (Lima-Junior et al., 2021). (iii) Regulatory roles. Given that ERVs contain specific regulatory sequences modulating their RNA synthesis and processing, they may provide new regulatory elements and thus influence the expression pattern of nearby host genes (Feschotte, 2008;Stoye, 2012). ERVs have been found to shape the evolution of transcriptional networks underlying immune response and placental function (Finnegan, 2012;Chuong et al. 2016Chuong et al. , 2017Pasquesi et al., 2020;Sun et al., 2021). Moreover, ERV-derived long noncoding RNA, designated lnc-EPAV, serves as a positive regulator of innate immune responses (Zhou et al., 2019). HERVH might act as an enhancer associated lncRNA that interacts with mediator kinase proteins and pluripotency-associated transcription factors to regulate the transcription of pluripotency-associated genes (Sexton et al., 2022). Besides, homologous recombination of ERVs at different chromosome regions can mediates chromosome rearrangement, which might play a major role in vertebrate genome evolution (Hughes, 2001).
Recently, a comprehensive analysis of vertebrate genomes identified 177 independent retroviral gene (gag and env) co-option events in vertebrates, suggesting ERV genes had been frequently co-opted during the evolution of vertebrates (Wang and Han, 2020). The evolution of retroviral gene co-option followed a birth and death model, in which retrovirus genes were frequently co-opted, and co-opted retroviruses were frequently lost Han, 2020, 2021b). These co-opted retroviral genes were subject to different selection pressure, implying potentially diverse cellular functionalities, but the specific cellular functionality remained to be studied (Wang and Han, 2020).

Outstanding questions
The discoveries of ERVs have greatly refined our understanding of the origin, deep evolution, and diversity of retroviruses. ERVs also shape the evolution of the genome complexity of their hosts. However, many outstanding questions remain: (i) Retroviruses have been thought to originate during the early evolution of vertebrates >450 Mya. However, it remains unclear how retroviruses originated (Hayward, 2017). (ii) Multiple classification systems of retroviruses have been developed, but all these systems have their own shortcomings. A classification system that fully integrates endogenous and exogenous retroviruses and follows phylogeny is eagerly needed. (iii) Retrovirus genes or regulatory sequences have been frequently co-opted in vertebrates, but their functionalities need to be experimentally characterized. (iv) The relationship between ERVs and human diseases should be fully evaluated. (v) ERVs are ubiquitously distributed in the vertebrate genomes, but only tens of exogenous retroviruses have been isolated. More efforts should be paid to isolate exogenous retroviruses in wild animals and assess their potential pathogenicity, which are crucial for predicting and controlling zoonotic disease outbreaks in the future. (vi) ERVs are highly abundant in vertebrates, but some vertebrate genomes harbor more ERVs than others (Hayward et al. 2013(Hayward et al. , 2015Xu et al., 2018;Zheng et al., 2021). It remains obscure about the factors that shape the landscape of ERVs within a certain genome.

Conflict of interest
The authors declare that they have no conflict of interest.

Acknowledgement
This work was supported by National Natural Science Foundation of China (31922001).


s2fieldsofstudy: Biology Medicine

text:
A pan-cancer analysis of the oncogenic role of ribonucleotide reductase subunit M2 in human tumors

Background Recent studies have identified ribonucleotide reductase subunit M2 (RRM2) as a putative promoter of tumors. However, no systematic analysis of its carcinogenicity has been conducted. Methods The potential functions of RRM2 in various tumor types were investigated using data from the Genotype-Tissue Expression (GTEx), the Clinical Proteomic Tumor Analysis Consortium (CPTAC), the Cancer Genome Atlas (TCGA), the Human Protein Atlas (HPA), cBioPortal, GEPIA, String, and Gene Set Enrichment Analysis (GSEA). We analyzed the difference in mRNA and protein expression, pathological stage, survival, mutation, tumor microenvironment (TME), and immune cell infiltration in relation to RRM2. Meanwhile, using TCGA and the Tumor Immune Estimation Resource 2 (TIMER 2), the associations between RRM2 expression, immune infiltration, and immune-related genes were assessed. Additionally, CCK-8, Edu and RT-PCR assays were used to validate that RRM2 acts as an oncogene in liver cancer cells and its association with HBx. A cohort of liver hepatocellular carcinoma (LIHC) patients (n=154) from Huashan Hospital was analyzed for the expression of RRM2 and the association between RRM2 and immune infiltration. Results Using the GTEx and TCGA databases, we discovered that 28 tumors expressed RRM2 at significantly higher levels than the corresponding normal tissues. Increased RRM2 expression may be predictive of a poor overall survival (OS) in patients with seven different cancers. GO, KEGG, and GSEA analyses revealed that the biological process of RRM2 was associated with the regulation of carcinogenic processes and immune pathways in a variety of tumor types. The expression of RRM2 was highly correlated with maker genes involved in immune activation and immunosuppression, immune checkpoints, DNA mismatch repair system (MMR), and the infiltration levels of Tregs and macrophages (TAMs), suggesting that the carcinogenic effect of RRM2 may be achieved by regulating immune related genes. Moreover, as demonstrated by CCK-8 and Edu assays, RRM2 was an oncogene in liver cancer cells. We confirmed for the first time that RRM2 was significantly upregulated by HBx, suggesting that RRM2 may be a key regulator of LIHC induced by HBV. IHC analysis validated the upregulated expression of RRM2 protein and its correlation with immune infiltration makers in a LIHC patient cohort. Conclusion RRM2 may be a valuable molecular biomarker for predicting prognosis and immunotherapeutic efficacy in pan-cancer, particularly in LIHC.

INTRODUCTION
The tumor immune microenvironment (TIME) is among the leading causes of malignant progression, involving a complex network of interactions between tumor cells and the immune system (Nelson et al., 2021;Cioni et al., 2020;Da Silva et al., 2021). Cancer cells can develop sophisticated strategies to evade the immune system's attack, most notably through the use of immune checkpoint-related genes (Zeng et al., 2021). Immune checkpoint-related genes may be involved in immunosuppressive mechanisms, posing a significant barrier to successful immunotherapy (Mo et al., 2020;Cao et al., 2020). Immune checkpoint blockade therapy, which revolutionized cancer treatment, has now become the gold standard for cancer immunotherapy (Abril-Rodriguez & Ribas, 2017). It may cause an immune system blockage, resulting in a long-term anti-cancer response (Palmieri & Carlino, 2018;Deng et al., 2021;Sheng et al., 2020;Xu et al., 2020). However, the response rate of immunotherapy patients with cancer varies, and only a small percentage of patients respond to immunotherapy, limiting its practical application (La-Beck et al., 2015). This could be due to a lack of biomarkers to direct personalized immunity toward a specific target (Hu et al., 2022). As a result, it is critical to identify potential therapeutic biomarkers for immune checkpoint inhibitors, which will aid in predicting treatment effectiveness by clearly separating responders from non-responders (Sharma et al., 2017;Friedman, Proverbs-Singh & Postow, 2016). Given that immunotherapy efficacy is related to TIME, we sought to identify oncogenes that play a critical role in TIME (Nishino et al., 2017).
RRM2, the small subunit of ribonucleotide reductase, has been identified as a tumor promoter (Zhan et al., 2021). Previous studies demonstrated a strong association between RRM2 expression level and cancer risk. RRM2 overexpression has been reported in patients with LIHC and may function as an endogenous ferroptosis controller (Mah et al., 2015). RRM2 overexpression in multiple myeloma patients' bone marrow mononuclear cells is highly associated with international staging system (ISS) staging, bone destruction, and extramedullary infiltration (Liu et al., 2020). RRM2 overexpression has also been implicated in oral squamous cell carcinoma, where it was associated with an advanced pathological grade and a poor prognosis for patients . RRM2 overexpression may be predictive of endometriosis malignancy (Yang et al., 2021). Despite these trends, further investigation is required to elucidate the additional roles of RRM2 in tumor progression and development, specifically its molecular contributions to tumor immunomodulation.
Given the variable and disease-specific nature of the molecular role, it is essential to investigate the profile of an oncogene in cases of generalized carcinoma (Majumder et al., 2021). So far, no systematic analysis of the carcinogenicity of RRM2 has been conducted. In this pan-cancer study, multiple databases were used to investigate the expression patterns of RRM2 and to visualize its prognostic value landscape. We identified the genes and signaling pathways regulated by RRM2 during cancer development. Then, we explored the connection between RRM2 expression and TIME, particularly immune checkpoints and tumor-infiltrating immune cells (TIICs).
Notably, the hepatitis B pathway was identified as a KEGG significant enrichment pathway of RRM2 binding and interacting genes, indicating that RRM2 plays an important role in the occurrence and progression of cancer, particularly LIHC. Therefore, HBVrelated liver cancer was selected to verify the biological function of RRM2. Combined, these results suggest that RRM2 may serve as a predictor of both cancer clinical outcomes and immunotherapy among a variety of cancers. and it may act as a carcinogen by modulating the tumor immune microenvironment, which justifies targeting RRM2 as a novel therapeutic approach.

Gene and protein expression analysis of RRM2 in pan-cancer
The UCSC Xena platform was used to retrieve RNA-seq data and clinical information for the TCGA transcriptome (Goldman et al., 2020). Using log2 (TPM), the expression data were transformed to eliminate duplicate and missing values (Vivian et al., 2017). In the Single Gene Analysis module of GEPIA2, the correlation between differential RRM2 mRNA expression analysis and pathological stages analysis of TCGA tumors was determined (Tang et al., 2017). The violin plots were created using log2 (TPM (transcripts per million) + 1) transformed expression data (Huo et al., 2021). The UALCAN portal was used to analyze data from CPTAC dataset. The RRM2 total protein or phosphoprotein level was compared (Chen et al., 2019;Qiu et al., 2021). The RRM2 protein expression was validated using IHC analysis from HPA database, the results were compared to TCGA RRM2 gene expression data (Uhlén et al., 2015).

Clinical value analysis of RRM2 in pan-cancer
To begin, Cox regression analyses were performed on TCGA datasets to determine the effect of RRM2 on pan-cancer patient survival with the R packages ''survival'',''survminer'' ''rms'' ''Hmisc'' and ''ggplot2''. We assessed OS, progression-free survival (PFI), and disease-specific survival (DSS). We visualized the univariate Cox regression analysis using the forestplot package and created a nomogram based on the optimal multivariate Cox regression analysis. GEPIA2 was utilized to obtain OS and DFS significant level map data as well as survival plots for RRM2 across tumors (Tang et al., 2019). To determine the hypothesis, the log-rank test was utilized. Then, the R packages ''pROC'' and ''ggplot2'' were employed to generate the receiver operating characteristic (ROC) curve using the RRM2 expression profile data from the TCGA databases. The area under the curve (AUC) was computed to quantify the ability to discriminate between the tumor and normal groups (Chen et al., 2021;Xie et al., 2021a;Xie et al., 2021b;Qiu et al., 2022).

Genetic alteration analysis of RRM2 in pan-cancer
We collected data on RRM2 alteration frequency, mutation type, mutated site information, CNA (Copy number alteration), DNA methylation alterations, and 3D structure of the protein throughout TCGA malignancies using the cBioPortal for Cancer Genomics platform (Cerami et al., 2012;Gao et al., 2013).

RRM2-related gene enrichment analysis
For the subsequent analysis of the protein-protein interaction (PPI) network, we employed the Search Tool for the Retrieval of Interacting Genes (STRING) (Szklarczyk et al., 2019). Using Cytoscape (version 3.8.2), the PPI network was constructed, visualized, and analyzed (Otasek et al., 2019). Using GEPIA2, the top 100 RRM2-associated target genes were identified, and Pearson correlation coefficients were calculated between RRM2 and selected genes. Meanwhile, using RNA-seq data from the TCGA database, the top ten correlation results in GEPIA2 correlation analysis were verified using the adjusted Spearman's rank correlation test.
A heatmap representation of the selected genes' expression profiles. Then, we used the ''clusterProfiler'' and ''org.Hs.eg.db'' and ''ggplot2'' R packages to merge and filter the top 100 RRM2-correlated genes data in order to conduct GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (Yu et al., 2012;Cui et al., 2020). In addition, GSEA analysis was utilized to explore the possible biological processes involving RRM2 in pan-cancer using ''clusterprofiler'' R package. The sets of pathway genes were curated using data from the Molecular Signatures Database (MSigDB). Enrichment significance was defined as gene sets with P adjust < 0.05 and FDR < 0.25 (Subramanian et al., 2015).

Association analysis of RRM2 expression with tumor immune microenvironment
To investigate the function of RRM2 in TIME, we determined the relationship between RRM2 and marker genes for immune activating genes, immunosuppressive genes, MMRrelated genes, and immune checkpoint genes using the R packages ''stat'' and ''ggplot2'' as well as Spearman's correlation analysis. TIMER is a tool for analyzing the relationship between RRM2 expression and immune infiltrates in all TCGA tumors (Li et al., 2020;Li et al., 2017).

Clinical samples
From 154 patients undergoing hepatic resection at Huashan Hospital, LIHC tumor tissues and paracancerous tissues were obtained. Then, all of the samples were examined using histology. Following that, histology was used to verify all of the samples. The Huashan Hospital Human Ethics Committee approved the experimental protocols, and all patients sign an informed consent form (HIRB: 2021-783).

Cell viability and cell growth assay
The CCK-8 and the 5-ethynyl-20-deoxyuridine (Edu) assay were used to determine cell viability and growth. According to manufacturer's protocol, HepG2 and Huh-7 cells were seeded at a density of 2000 cells/well in 96-well plates and treated with negative control and si-RRM2. Cell viability of transfected cells at 48 h was detected using CCK-8 assay (Dojindo, Rockville, MD, USA). Furthermore, the cells were incubated with Edu (Ribobio, Guangzhao, China). Data were obtained from three separate experiments.

Statistical analysis
R language version 3.6.2 was utilized for statistical analysis. This study displayed all data as means ± standard deviation. Student's t test and GraphPad Prism 6 were applied to determine the significance of the differences.

RRM2 is upregulated in pan-cancer
Using GEPIA2, we analyzed the mRNA expression status of RRM2 across numerous cancer types. RRM2 expression was found to be elevated in the majority of human tumors, according to the findings (Fig. S1). RNA-seq data from the TCGA and GTex databases were used to further confirm the expression of RRM2 in pan-cancer. As depicted in Fig. 1A, RRM2 expression significantly overexpressed among 28 types of cancer. In contrast, it showed little expression in acute myeloid leukemia (LAML). In paired tumor and normal tissues from the TCGA pan-cancer database, RRM2 was highly expressed in bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), thyroid carcinoma (THCA), uterine corpus (UCEC), esophageal carcinoma (ESCA), head and neck cancer (HNSC), kidney renal papillary cell carcinoma (KIRP), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), LIHC, lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC) and prostate adenocarcinoma (PRAD) (Fig. 1B).
We investigated the relationship between RRM2 expression and pathological stages of cancer. Remarkably, we discovered clear associations in the majority of instances, including adrenocortical carcinoma (ACC), BRCA, COAD, kidney chromophobe (KICH), skin cutaneous melanoma (SKCM), KIRP, LUAD, kidney renal clear cell carcinoma (KIRC), LUSC, ovarian serous (OV), LIHC and THCA (Fig. 1C, all P < 0.05). In addition to evaluating RRM2 at the transcriptional level, we evaluated it at the protein level using the National Cancer Institute's CPTAC dataset, which contains large-scale proteome data. The total protein expression of RRM2 was strongly elevated in BRCA, KIRC, COAD, HNSC, LIHC, LUAD, OV, pancreatic adenocarcinoma (PAAD), and UCEC patients ( Fig. 2A, P < 0.05, P < 0.001). Meanwhile, we examined IHC information from the HPA database to validate RRM2 protein expression in CPTAC dataset and compared the results to TCGA data on RRM2 gene expression. The outcomes of the analyses of data from the HAP and TCGA database were congruent. Normal liver, thyroid gland, stomach, and colon were negative for RRM2 immunohistochemical staining, whereas tumor tissues were moderately stained (Fig. 2B). Regarding the protein phosphoprotein of RRM2, we analyzed three types of pan-cancer tumors in greater detail using the CPTAC dataset. Figure 2C summarizes RRM2 phosphorylation sites and significant differences. The S20 locus of RRM2 is more highly phosphorylated in BRCA than in normal tissues, but less so in LUAD.

Multifaceted clinical value analysis of RRM2 in Cancers
We downloaded the TCGA mRNA sequencing and clinical data of 33 cancer types from the UCSC Xena platform and used the Cox proportional hazards model to examine the prognostic assessment value of RRM2 in pan-cancer (Table S1). Figure S2 shows the forest plots indicating RRM2 as a high-risk factor for KIRP, KIRC, ACC, LIHC, LUAD, PAAD, brain lower grade glioma (LGG) and PRAD of OS prognosis assessment in pan-cancer, and RRM2 expression was clearly correlated with worse DSS in LIHC, ACC, BLCA, PAAD, PRAD, KIRP, LGG, LUAD, pheochromocytoma and paraganglioma (PCPG) and KIRC. A nomogram integrating RRM2 expression and pathological stage, an independent clinical risk factor, was developed (Fig. 3). The C-index in ACC, KIRC, and KICH indicated the prediction results match the observation results perfectly, as shown in Table S2. Using the Mantel-Cox test and the GEPIA2 database, we analyzed the survival contribution of RRM2 in various cancer types; the survival maps supplemented by OS and disease-free survival (DFS) curves are depicted in Fig. S3. High RRM2 transcription levels were associated with a poor OS prognosis for cancers, including cervical squamous cell carcinoma and endocervical adenocarcinoma (

Enrichment analysis of RRM2-related partners
Using the STRING tool, we obtained experimentally detected RRM2-binding proteins, and the PPI network of proteins was shown in Fig. 6A. The GEPIA2 tool was then used to identify the top 100 genes that correlate with RRM2 expression in combining TCGA pan-cancer expression. As illustrated in Fig. 6B R = 0.72) genes (all P < 0.001). The majority of tumors exhibited a positive relationship between RRM2 and the nine genes listed above, as indicated by the results of a heat map (Fig. 6C). As shown in Fig. 6D, GO analysis determined that RRM2 was primarily involved in cell cycle checkpoint, ATPase activity, and condensed chromosome. Several pathways, including ''Hepatitis B'', ''p53 signaling pathway'', and ''cell cycle'' were identified as the KEGG pathways most significantly enriched, suggesting that RRM2 played a crucial part in the development and progression of cancers, particularly HBV-related LIHC (Fig. 6E). Meanwhile, GSEA was used to identify functional RRM2 enrichment (Fig. S4). High expression of RRM2 is implicated in numerous carcinogenic and immune-related processes, including PLK1 pathways, Forkhead box M1 (FOXM1) pathways, P53 signaling pathways, immunoregulatory interactions between a lymphoid and a nonlymphoid cell, and DNA repair and methylation.
As for immune infiltration analysis, RRM2 expression level was linked with the estimated infiltration value of Tregs for PRAD, BRCA, LIHC, and THCA (Fig. 8A). In BLCA, COAD, ESCA, KIRP, PAAD, and STAD, RRM2 expression correlated positively with M1-like macrophages and negatively with M2-like macrophages (Fig. 8B). RRM2 expression was significantly correlated with the immune-associated cells T cells in four cancer types, cytotoxic cells in four cancer types, B cells in two cancer types, CD8+T cells in three cancer types, Treg cells in four cancer types, dendritic cells in one cancer type, macrophages in three cancer types, NK cells in six cancer types, and neutrophils in three cancer types. In particular, we discovered a significant correlation between the RRM2 expression level and seven types of immune-associated cells in THYM (Fig. S5).

The validation of RRM2 expression and its oncogenicity
IHC was used to confirm the expression of the RRM2 protein in LIHC patient cohort from Huashan Hospital. We verified the upregulated protein expression of RRM2 in tumor tissue (n = 154) compared to adjacent normal tissue by IHC staining (Figs. 9A and 9B). To demonstrate the oncogene function of RRM2, we conducted experiments with liver cancer cells. We successfully suppressed RRM2 expression in HepG2 and Huh-7 cells using three RRM2 siRNA (Figs. 9C and 9D). Based on the efficiency of knockout, RRM2 siRNA-1 was chosen for subsequent experiments. CCK-8 (Fig. 9E) and Edu assays (Figs. 9F-9H) showed that RRM2 knockdown significantly inhibited cell proliferation in HepG2 and Huh-7, confirming our conclusion. In addition, HBx plasmid was transfected into LO2 cells to explore the regulatory relationship of HBx on RRM2 (Fig. 9I). We found that RRM2 mRNA expression were significantly upregulated in LO2-HBx cells (Fig. 9J).

Validating the associations of RRM2 with immunological makers of Tregs in LIHC patient cohort
Tumors tended to have a high abundance of Tregs cell infiltration, which was associated with protumor immunity. IHC was used to examine the expression of transforming growth factor β1 (TGF-β1) and activator of transcription (STAT)-5B protein from patients tissues.
We demonstrated that both TGF-β1 and STAT5B expression was higher in the RRM2-high group than that in the RRM2-low group (Figs. 10A and 10B). TGF-β1 was significantly more strongly positive in LIHC tissues than in adjacent paracancerous tissues. STAT5B expression followed a similar pattern (Fig. 10C). A significant correlation between STAT5B and RRM2 protein levels was observed in LIHC and paracanerous tissues, as well as mixed tissues from LIHC patients (Fig. 10D). A positive correlation between TGF-β1 and RRM2 was also observed in mixed tissues from LIHC patients (Fig. 10E). All these findings suggested that RRM2 may play a significant role in immune infiltration in LIHC, at least in part.

DISCUSSION
A pan-cancer analysis with a comprehensive understanding of human cancer was typically used to identify recurring characteristics and heterogeneity during vital biological processes that contribute to a dysregulated tumor microenvironment, which is essential for identifying specific targets or characteristics for more precise and tailored treatments (Ye et al., 2019;Barger et al., 2019). This is the first pan-cancer analysis focusing on the significance of RRM2, to our knowledge. Pan-cancer analysis is of clinical value for identifying the similarities, and differences of RRM2 in various tumors. Specifically, the different results of survival analysis and the association with immune-related genes could serve as a guide and theoretical foundation for future clinical treatment, particularly tumor immunotherapy. We analyzed the expression profiles of RRM2 and its association with prognosis in various cancers. Our data revealed that RRM2 expression was significantly The correlation between STAT5B and RRM2 positive area in paracancerous and liver tumor tissues, respectively, and the total correlation in liver and paracancerous tissues. (E) The correlation between TGF-β1 and RRM2 positive area in paracancerous and liver tumor tissues, respectively, and the total correlation in liver and paracancerous tissues. The differences between two groups were estimated by the Student's t test.The Pearson correlation coefficient was used to evaluate the correlation. All values are the mean ± SD. *P < 0 .05, **P < 0.01, ***P < 0 .001, ****P < 0 .0001. Full-size DOI: 10.7717/peerj.14432/ fig-10 up-regulated compared to corresponding noncancerous tissues in a variety of cancers, regardless of gene or protein expression levels, indicating the extensive oncogenic properties of RRM2 in cancers and promising prospects for cancer research. Since phosphorylationdephosphorylation cascade is recognized as a crucial step in oncogenesis. Consequently, we then analyzed the phosphorylation level of RRM2. In comparison to normal tissues, the RRM2 S20 phosphorylation level was significantly altered in BRCA and LUAD. The clinical significance of these posttranslational modification sites remains unknown. Cox proportional hazards model and Kaplan-Meier analysis were used to comprehensively evaluate the association between RRM2 expression and survival time. The results of univariate and multivariate Cox regression analyses suggested that increased RRM2 expression may be associated with shorter OS, DSS, and PFI among KIRP, ACC, KIRC, PAAD, LGG, LUAD, LIHC and PRAD. Kaplan-Meier survival curves of prognostic value analysis revealed that high RRM2 expression predicted poorer OS in ACC, KIRC, KIRP, LGG, LIHC, LUAD, and PAAD, indicating that RRM2 is a risk factor in these cancers. Notably, these results identified RRM2 as a carcinogenic marker for ACC, LGG, KIRP, LIHC, KIRC, PAAD, and LUAD prognosis, regardless of the prognostic algorithm. In addition, ROC curve analysis revealed that RRM2 performed well in the diagnosis of as many as sixteen different types of tumors, especially in CHOL (AUC = 1, 95% CI [1-1]) and UCS (AUC = 1, 95% CI [1-1]). These findings suggest that high RRM2 expression is primarily responsible for carcinogenesis in the majority of tumor types, adding fresh insight to the primary pathogenesis of various tumors and providing an innovative clinical biomarker that can predict the overall survival of cancer patients.
Cancers in humans result from the accumulation of genetic mutations. Thus, we investigated RRM2 genetic mutations in human tumor samples. We discovered that mutation and amplification were the most prevalent RRM2 alterations in pan-cancer, which may trigger the onset and development of cancer. In addition, high RRM2 DNA methylation could serve as a prognostic biomarker in multiple cancers. We used pathway enrichment analyses on known RRM2-interacting proteins and RRM2-correlated genes to uncover the molecular mechanism of RRM2 in tumorigenesis and development. Results indicate that RRM2 was significantly overrepresented in immune-related pathways, specially, immunoregulatory interactions between lymphoid and non-lymphoid cell pathways are a notable immune-related pathway that modify the response of lymphoid cells to self, tumor, and pathogen antigens (Cemerski & Shaw, 2006;Nedvetzki et al., 2007). These findings indicate that RRM2 regulates the tumor immune microenvironment and ligand-receptor interactions between lymphoid and malignant tumor cells.
Next, we investigated the effect of RRM2 on TIME and evaluated whether RRM2 could be used as a biomarker of tumor cell response to immunotherapy and influence clinical outcomes via TIME parameters. Our initial investigation revealed that aberrantly high RRM2 expression was negatively related to immune activating maker genes in the majority of cancers, whereas it was positively correlated with immunosuppressive genes, indicating that RRM2 regulated the tumor immune microenvironment and immuno-oncological interactions. RRM2 expression was highly associated with immune checkpoint genes and MMR-related genes, indicating that RRM2 may regulate the tumor immune microenvironment to promote tumor growth. As principal components of TIME, tumor-infiltrating immune cells frequently participated in tumor behaviors such as Initiation, progression, and metastasis of cancer (Kim et al., 1991). Tumor-infiltrating immune cells are dysfunctional, unable to control tumor growth, and may even foster its advancement, resulting in immune evasion. The link between RRM2 and immune cell infiltration was investigated in this study, we observed a positive correlation between Treg infiltration levels and RRM2 in BRCA, KIRC, LIHC, PCPG, PRAD, and THCA, indicating that RRM2 contributes to Tregs exerting an inhibitory role in aiding malignant tumor cells in evading attack by cytotoxic CD8+ T cells (Shu et al., 2019;Liu, Workman & Vignali, 2016). Meanwhile, we assessed the connection between RRM2 expression and macrophage infiltration (Togashi, Shitara & Nishikawa, 2019). Infiltration of tumorassociated macrophages (TAMs) was significantly correlated with RRM2 expression in the majority of tumor types. In addition, RRM2 may directly or indirectly contribute to macrophage polarization and tumor immune evasion in BLCA, COAD, ESCA, KIRP, PAAD, and STAD (Armitage et al., 2021;. RRM2 could be an option for reducing the number of TAMs. Our finding that RRM2 played a crucial role in cancer immunology emphasized the need for future research involving larger patient cohorts to determine the clinical utility of immune checkpoint inhibitor biomarkers. It is worth noting that KEGG pathway analysis revealed that RRM2 and HBV have a strong relationship. HBV, the leading cause of acute and chronic viral hepatitis, may contribute to the development of LIHC. In addition, the liver, as an important immune organ, provides a large and complex immune microenvironment for the development of LIHC, particularly viral hepatitis involving immune disorders, which plays a significant role in promoting the development of hepatocellular carcinoma (Xie et al., 2022;Vivier et al., 2008). Therefore, we hypothesize that RRM2 is a crucial intermediate link in the development of HBV-related LIHC. On the basis of this hypothesis, we selected LIHC as a representative tumor for future trials to validate the carcinogenicity of RRM2 and its association with HBV as well as immune infiltration. IHC analysis validated the upregulated expression of RRM2 protein and a strong positive association of RRM2 with immunological makers (TGF-β1 and STAT5B) of Tregs in LIHC patients cohort, which suggested that RRM2 may be a novel immune cell infiltration regulator that can control the recruitment and activation of immune cells in LIHC. HBV X protein (HBx), functioning as a viral transactivation factor, can alter multiple genes and result in the development of LIHC (Gao et al., 2021;Niu et al., 2017;Xue et al., 2020;Zhang et al., 2021). We investigated the relationship between RRM2 and HBV by analyzing the level of RRM2 expression after HBx overexpression. Notable is the fact that we confirmed for the first time that RRM2 was significantly upregulated by HBx, suggesting that RRM2 may serve as a key regulator of LIHC induced by HBV. IHC analysis also validated the upregulated expression of RRM2 protein and its correlation with immune infiltration makers in a LIHC patient cohort. Besides, we discovered that inhibiting RRM2 significantly reduced the proliferation rate of HepG2 and Huh-7 cells, which demonstrated RRM2 was an oncogene in LIHC.
Above all, RRM2 may be a valuable molecular biomarker for predicting prognosis and immunotherapeutic efficacy in pan-cancer, particularly in LIHC. We believe these findings may serve as a foundation for future functional experiments and may have eventual clinical implications. Our research has a number of limitations. Importantly, additional in vivo experiments and clinical trials should be conducted to validate RRM2's role as an immune-related biomarker.

CONCLUSIONS
Our pan-cancer analysis provides a comprehensive summary of the oncogenic roles of RRM2 in a variety of human cancers, particularly its role in LIHC and the association between RRM2 and a poor prognosis. We have shown for the first time that HBx significantly upregulates RRM2, indicating that RRM2 may be a key regulatory gene of HBV-induced LIHC. In addition, high RRM2 expression may contribute to the tumor immune microenvironment, suggesting that the carcinogenic effect of RRM2 may be realized through the regulation of immune-related genes. Consequently, RRM2 may be a useful biomarker for the prognosis and prediction of immunotherapeutic efficacy.

TGF-β1
Transforming growth factor β1 STAT Activator of transcription


Agarose-Degrading Characteristics of a Deep-Sea Bacterium Vibrio Natriegens WPAGA4 and Its Cold-Adapted GH50 Agarase Aga3420

Up until now, the characterizations of GH50 agarases from Vibrio species have rarely been reported compared to GH16 agarases. In this study, a deep-sea strain, WPAGA4, was isolated and identified as Vibrio natriegens due to the maximum similarity of its 16S rRNA gene sequence, the values of its average nucleotide identity, and through digital DNA–DNA hybridization. Two circular chromosomes in V. natriegens WPAGA4 were assembled. A total of 4561 coding genes, 37 rRNA, 131 tRNA, and 59 other non-coding RNA genes were predicted in the genome of V. natriegens WPAGA4. An agarase gene belonging to the GH50 family was annotated in the genome sequence and expressed in E. coli cells. The optimum temperature and pH of the recombinant Aga3420 (rAga3420) were 40 °C and 7.0, respectively. Neoagarobiose (NA2) was the only product during the degradation process of agarose by rAga3420. rAga3420 had a favorable stability following incubation at 10–30 °C for 50 min. The Km, Vmax, and kcat values of rAga3420 were 2.8 mg/mL, 78.1 U/mg, and 376.9 s−1, respectively. rAga3420 displayed cold-adapted properties as 59.7% and 41.2% of the relative activity remained at 10 3 °C and 0 °C, respectively. This property ensured V. natriegens WPAGA4 could degrade and metabolize the agarose in cold deep-sea environments and enables rAga3420 to be an appropriate industrial enzyme for NA2 production, with industrial potential in medical and cosmetic fields.

Introduction
Agar is the main component in the cell walls of red algae [1]. As the skeleton structure of agar [2], agarose possesses alternately repeated units of β-D-galactose and 3,6-anhydro-L-galactose [3]. Agarose can be degraded into oligosaccharides and monosaccharides with the potential for anti-inflammatory, anti-cancer, anti-caries, whitening, and moisturizing bio-activities in the industries of medicine, health products, cosmetics, and energy [4]. In addition, the polysaccharides (including cellulose, chitin, and agarose) are the carbon sources in environments, and the degraders of these polysaccharides are also participants in the carbon cycle [5,6]. Therefore, the degradation of agarose is of great value in the studies of industry production [5] and the carbon cycle [7,8].
Microbial degradation is the ideal method for the industrial degradation of agarose due to its high efficiency, environmental protection, energy conservation, and stable reaction product [9]. Agarases, including αand β-agarase [1], are the molecular fundamentals of agarose-degrading microbes. α-agarase, belonging to the GH96 family, cuts the α-1,3-glycosidic bond and degrades agarose into agaro-oligosaccharides (AOS) with 3,6-anhydro-L-galactose at the reducing ends of the product. β-agarase, belonging to the GH16, GH39, GH50, GH42, GH86, and GH118 [10] families [1,11,12], cuts the . Figure 1. The genome analysis of Vibrio natriegens WPAGA4. The circular maps of the two chromosomes of V. natriegens WPAGA4 are shown in (A) and (B). The outer layer of the map is the size of the genome; the second and third circles are the genes on the positive and negative strands, respectively. The colors stand for the COG function classifications; the repetition sequences are shown in the fourth circle; the fifth circle is tRNA (blue) and rRNA (purple); the sixth circle represents the GC content. The light-yellow part indicates that the GC content of the region is higher than the average GC content of the genome, and the blue part indicates that the GC content of the region is lower than the average GC content of the genome; the innermost circle stands for the GCskew. Dark gray represents the region with G content greater than C, and red represents the region with C content greater than G. The CAZyme gene annotation results are shown in (C).

The Determination and Optimization of the Agarose-Degrading Activity in the Broth
The optimum culture temperature for the agarose-degrading activity in the culture broth of V. natriegens WPAGA4 was 31 °C (Figure 2A). Additional agar improved the agarose-degrading activity to a certain extent ( Figure 2C). Furthermore, additional amylum and lactose improved the agarose-degrading activity, whereas glucose, sucrose, maltose, and glycerin could inhibit the activity in WPAGA4 broth ( Figure 2C). The additional nitrogen sources presented in this work could not increase the agarose-degrading activity in the broth. Notably, NH4NO3 and urea obviously inhibited the agarose-degrading activity in the broth of WPAGA4 ( Figure 2D). The maximum agarose-degrading activity of V. natriegens WPAGA4 in the fermentation broth was 19.49 U/mL. . The outer layer of the map is the size of the genome; the second and third circles are the genes on the positive and negative strands, respectively. The colors stand for the COG function classifications; the repetition sequences are shown in the fourth circle; the fifth circle is tRNA (blue) and rRNA (purple); the sixth circle represents the GC content. The light-yellow part indicates that the GC content of the region is higher than the average GC content of the genome, and the blue part indicates that the GC content of the region is lower than the average GC content of the genome; the innermost circle stands for the GC-skew. Dark gray represents the region with G content greater than C, and red represents the region with C content greater than G. The CAZyme gene annotation results are shown in (C).

The Determination and Optimization of the Agarose-Degrading Activity in the Broth
The optimum culture temperature for the agarose-degrading activity in the culture broth of V. natriegens WPAGA4 was 31 • C ( Figure 2A). Additional agar improved the agarose-degrading activity to a certain extent ( Figure 2C). Furthermore, additional amylum and lactose improved the agarose-degrading activity, whereas glucose, sucrose, maltose, and glycerin could inhibit the activity in WPAGA4 broth ( Figure 2C). The additional nitrogen sources presented in this work could not increase the agarose-degrading activity in the broth. Notably, NH4NO3 and urea obviously inhibited the agarose-degrading activity in the broth of WPAGA4 ( Figure 2D). The maximum agarose-degrading activity of V. natriegens WPAGA4 in the fermentation broth was 19.49 U/mL.

Analyses of the Genome Sequence and Putative Agarase Gene aga3420 of Strain WPAGA4
Two circular chromosomes in V. natriegens WPAGA4 containing 3,241,580 bp and 1,994,952 bp were assembled ( Figure 1A,B). A total of 4561 coding genes, 37 rRNA, 131 tRNA, and 59 other non-coding RNA genes were predicted in the genome of V. natriegens WPAGA4, and 5 CRISPRs, 12 gene islands, and 1 prophage were annotated in the genome sequence.

Expression and Characterization of the Recombinant Agarases
Among the four recombinant agarases, only the recombinant Aga3420 (rAga3420) showed obvious agarase activity under the conditions described in Section 4.6. Therefore, rAga3420 was selected for further analysis. The phylogenetic analysis showed that the amino acid sequence of Aga3420 had the closest genetic distance to Agarivorans gilvus agarase (accession number: AQT38174.1) (Figure 3).

Analyses of the Genome Sequence and Putative Agarase Gene Aga3420 of Strain WPAGA4
Two circular chromosomes in V. natriegens WPAGA4 containing 3,241,580 bp and 1,994,952 bp were assembled ( Figure 1A,B). A total of 4561 coding genes, 37 rRNA, 131 tRNA, and 59 other non-coding RNA genes were predicted in the genome of V. natriegens WPAGA4, and 5 CRISPRs, 12 gene islands, and 1 prophage were annotated in the genome sequence.

Expression and Characterization of the Recombinant Agarases
Among the four recombinant agarases, only the recombinant Aga3420 (rAga3420) showed obvious agarase activity under the conditions described in Section 4.6. Therefore, rAga3420 was selected for further analysis. The phylogenetic analysis showed that the rAga3420 with a molecular weight of ~100 kDa, which was consistent with its theoretical molecular weight (104.87 kDa), was purified based on its 6xHis tag in this work ( Figure 4A). The optimum temperature of rAga3420 was 40 °C ( Figure 4B). In addition, 59.7% and 41.2% of the relative rAga3420 activity remained at 10 °C and 0 °C, respectively ( Figure 4B). The activity of rAga3420 is rapidly lost at temperatures higher than 40 °C ( Figure 4B). The optimum pH of rAga3420 was 7.0 ( Figure 4C). rAga3420 had good stability after the incubation at 10-30 °C for 50 min ( Figure 4D); however, 22.3% of the relative activity was lost after the incubation at 30 °C for 60 min ( Figure 4D). The characterization comparison of rAga3420 and other GH50 agarases are listed in Table S1. The metal ions in this work could improve the rAga3420 activity ( Figure 4E). Dithiothreitol (DTT) could increase the activity of rAga3420, whereas sodium dodecyl sulfate (SDS) and Zn 2+ inhibited its activity ( Figure 4F). Salinity had an obvious negative effect on rAga3420 activity ( Figure 4G). The values of Km, Vmax, and kcat were 1.5 ± 0.6 mg/mL, 53.8 ± 1.4 U/mg, and 259.2 ± 21.5 s −1 , respectively, according to the Lineweaver-Burk plot ( Figure S2).  rAga3420 with a molecular weight of~100 kDa, which was consistent with its theoretical molecular weight (104.87 kDa), was purified based on its 6xHis tag in this work ( Figure 4A). The optimum temperature of rAga3420 was 40 • C ( Figure 4B). In addition, 59.7% and 41.2% of the relative rAga3420 activity remained at 10 • C and 0 • C, respectively ( Figure 4B). The activity of rAga3420 is rapidly lost at temperatures higher than 40 • C ( Figure 4B). The optimum pH of rAga3420 was 7.0 ( Figure 4C). rAga3420 had good stability after the incubation at 10-30 • C for 50 min ( Figure 4D); however, 22.3% of the relative activity was lost after the incubation at 30 • C for 60 min ( Figure 4D). The characterization comparison of rAga3420 and other GH50 agarases are listed in Table S1. The metal ions in this work could improve the rAga3420 activity ( Figure 4E). Dithiothreitol (DTT) could increase the activity of rAga3420, whereas sodium dodecyl sulfate (SDS) and Zn 2+ inhibited its activity ( Figure 4F). Salinity had an obvious negative effect on rAga3420 activity ( Figure 4G). The values of Km, Vmax, and kcat were 1.5 ± 0.6 mg/mL, 53.8 ± 1.4 U/mg, and 259.2 ± 21.5 s −1 , respectively, according to the Lineweaver-Burk plot ( Figure S2). rAga3420 with a molecular weight of ~100 kDa, which was consistent with its theoretical molecular weight (104.87 kDa), was purified based on its 6xHis tag in this work ( Figure 4A). The optimum temperature of rAga3420 was 40 °C ( Figure 4B). In addition, 59.7% and 41.2% of the relative rAga3420 activity remained at 10 °C and 0 °C, respectively ( Figure 4B). The activity of rAga3420 is rapidly lost at temperatures higher than 40 °C ( Figure 4B). The optimum pH of rAga3420 was 7.0 ( Figure 4C). rAga3420 had good stability after the incubation at 10-30 °C for 50 min ( Figure 4D); however, 22.3% of the relative activity was lost after the incubation at 30 °C for 60 min ( Figure 4D). The characterization comparison of rAga3420 and other GH50 agarases are listed in Table S1. The metal ions in this work could improve the rAga3420 activity ( Figure 4E). Dithiothreitol (DTT) could increase the activity of rAga3420, whereas sodium dodecyl sulfate (SDS) and Zn 2+ inhibited its activity ( Figure 4F). Salinity had an obvious negative effect on rAga3420 activity ( Figure 4G). The values of Km, Vmax, and kcat were 1.5 ± 0.6 mg/mL, 53.8 ± 1.4 U/mg, and 259.2 ± 21.5 s −1 , respectively, according to the Lineweaver-Burk plot ( Figure S2).
solution, IPTG-induced cell lysis solution, and purified rAga3420, respectively. The effects of temperature (B) and pH (C) on the relative activity of rAga3420. The stability (D) of rAga3420. The effects of metal ions (E), chemical agents (F), and salinity (G) on the relative activity of rAga3420. DTT, dithiothreitol; ME, β-mercaptoethanol.

Discussion
Polysaccharide degradation is an important component of the carbon cycle in marine environments [23], including the deep sea [24]. Microorganisms are crucial players in driving polysaccharide degradation, which participates in the carbon cycle of deep-sea environments. The genus Vibiro, such as V. diabolicus [25], V. profundi [26], V. bathopelagicus [27], and other species that remain unidentified [28,29], have been widely isolated in the deep sea. Of them, several strains were reported with polysaccharide-degrading activities. For example, Vibrio sp. JAM-A9m, isolated from the deep sea, was detected with alginate lyase activity [29], and another deep-sea Vibrio strain showed cold-adapted amylase activity [30] in previous works. Therefore, the genus Vibrio, from the deep-sea environment, is a contributor to polysaccharide degradation. This work reported the first Vibrio strain isolated from the deep-sea environment with agarose-degrading activity. Although the agarase activity of V. natriegens WPAGA4 (19.49 U/mL) was moderate among the agarose-degrading strains, such as Agarivorans sp. JAMB-A11 (10.3 U/mL) [31], Vibrio sp. PO-303 (63.6 U/mL) [32], and Pseudomonas sp. SK38 (32.3 U/mL) [33], the current study indicated that this strain was potentially involved in the polysaccharide degradation and carbon cycle processes inside the deep-sea environment due to its agar-degrading ability.
As far as we know, the agarose concentration in deep-sea environments has not been reported, and the reason that an agarose-degrading strain exists in the deep sea remains

Discussion
Polysaccharide degradation is an important component of the carbon cycle in marine environments [23], including the deep sea [24]. Microorganisms are crucial players in driving polysaccharide degradation, which participates in the carbon cycle of deep-sea environments. The genus Vibiro, such as V. diabolicus [25], V. profundi [26], V. bathopelagicus [27], and other species that remain unidentified [28,29], have been widely isolated in the deep sea. Of them, several strains were reported with polysaccharide-degrading activities. For example, Vibrio sp. JAM-A9m, isolated from the deep sea, was detected with alginate lyase activity [29], and another deep-sea Vibrio strain showed cold-adapted amylase activity [30] in previous works. Therefore, the genus Vibrio, from the deep-sea environment, is a contributor to polysaccharide degradation. This work reported the first Vibrio strain isolated from the deep-sea environment with agarose-degrading activity. Although the agarase activity of V. natriegens WPAGA4 (19.49 U/mL) was moderate among the agarose-degrading strains, such as Agarivorans sp. JAMB-A11 (10.3 U/mL) [31], Vibrio sp. PO-303 (63.6 U/mL) [32], and Pseudomonas sp. SK38 (32.3 U/mL) [33], the current study indicated that this strain was potentially involved in the polysaccharide degradation and carbon cycle processes inside the deep-sea environment due to its agar-degrading ability.
As far as we know, the agarose concentration in deep-sea environments has not been reported, and the reason that an agarose-degrading strain exists in the deep sea remains unexplored until now. However, the WPAGA4 strain is not the only strain isolated from the deep-sea environment with agar-degrading activity. In former studies, deep-sea strains, including Flammeovirga sp. OC4 [13,16], Shewanella sp. WPAGA9 [16], Flammeovirga pacifica presence of these strains indicates that agar degradation might potentially be a part of the carbon cycle process within the deep-sea environment.
Although V. natriegens isolated from other environments have been used for agarase production in several studies, the agarose-degrading activity and the agarase of V. natriegens have been studied less [36,37]. The current work characterized a GH50 agarase, Aga3420, which supported the agarose-degrading activity of V. natriegens WPAGA4. Similar to other GH50 agarases, rAga3420 was an exo-type that degraded agarose into NA2 [1,38], which is an essential step for microorganisms to metabolize agarose as energy material [39]. As far as we know, rAga3420 is the first GH50 agarase from the deep sea with characterized enzyme properties, including optimum temperature, pH, and thermostability. Compared with the GH50 agarases from other environments, rAga3420 possessed outstanding coldadapted properties (Table S1). For example, two GH50 agarases, AgaA50 and AgaC50, in Microbulbifer elongatus PORT2 isolated from Indonesian coastal seawater lost~80% of its relative activity at 30 • C [40], and GH50 agarase Aga575 from the surface of fresh porphyra only retained <40% of its relative activity at 25 • C [15]. However, rAga3420 could maintain more than 40% and 50% of the relative activities at 0 • C and 10 • C, respectively, thereby ensuring that V. natriegens WPAGA4 could use agarose and agar as energy materials at low temperatures and could thus drive the polysaccharide cycle in the harsh, cold environment of the deep sea [41].
Recently, scientific interest in cold-adapted enzymes has been increasing [42]. Coldadapted enzymes can work at low and moderate temperatures (<40 • C), thereby saving energy expenditure and equipment optimization costs in the heating process [43]. In addition, the by-products of the reaction can be reduced at low temperatures [44]. The cold-adapted rAga3420 in this work could serve as an energy-and cost-saving bio-source for industrial agarose degradation and oligosaccharide production [12]. Meanwhile, low temperatures could also prevent caramelization, a potential by-product of high-temperature polysaccharide degradation [45]. Therefore, the agarase rAga3420, with its low temperature activity, is not only a potential driver of the polysaccharide cycle, but it is also an ideal bio-tool for oligosaccharide production in industrial applications [12].
Although the NAOS, with different degrees of polymerizations, such as NA2, NA4, NA6, and NA8, are reported as active substances [46], NA2 has its advantages in several bioactivities. For example, the antioxidant activity of NA2 is better than that of NA4 and NA6 as it can rapidly remove 2,2-diphenyl-1-picrylhydrazyl [47]. Therefore, the preparation and production of NA2 is important for the development of NAOS industries. However, the endo-type agarases that produce NA2 degrade agarose into mixture products [12], thereby increasing the cost of purification. The exo-type GH50 agarase rAga3420 only produces NA2 during the degradation process. The single product of rAga3420 lowered the cost and reduced process complexity for the purification and production of NA2. In addition, the high optimum temperatures of agarases [13,48,49] increase the energy input during NAOS production. The cold-adapted properties of rAga3420 showed an optimum temperature of 40 • C and retained more than 70% of the relative activity at 20 • C and 30 • C, indicating that rAga3420 could produce NA2 at room temperature. These properties reduced the energy costs of NA2 production. These characterizations make rAga3420 an ideal bio-tool for agarose degradation and NA2 production.
A total of four putative agarase genes were expressed in this work; however, only rAga3420 was detected to have the agarase activity under the conditions provided in this work. We attributed this to the false-positive prediction, minimal expression, protein misfold, and incongruous reaction conditions and substrates in the former work [50]. Nevertheless, it still reminds us that there are several unknown agarase resources in this strain and in deep-sea environments, which remain to be explored in future works.
In conclusion, a V. natriegens WPAGA4 strain with agarose-degrading activity was isolated from deep-sea sediments in the current work. A GH50 agarase gene, aga3420, was sequenced in the genome of V. natriegens WPAGA4 and was expressed in E. coli cells. The exo-type rAga3420 showed a cold-adapted property that maintained~40% of the relative activity at 0 • C. This property ensured that V. natriegens WPAGA4 could degrade and metabolize agarose in a cold, deep-sea environment, and enable rAga3420 to be an appropriate industrial enzyme for NA2 production with further industrial potential in medical and cosmetic fields.

Isolation, Purification, and Identification
The deep-sea sediment sample was collected from the Western Pacific at site CM3MC04-1 (157 • 24 31 E, 19 • 30 30 N; 5378 m deep) in July 2007 using a core sampler. According to the practical guide [51],~2 cm sediments were removed from the outer layer of the sample with a sterilized knife. Then, only the core part was stored in sterile centrifuge tubes which were stored in a −80 • C ultra-low temperature freezer. After the samples were transported back to the labs, they were serially diluted using the liquid 2216E medium (5.0 g of peptone, 1.0 g of yeast extracts, 0.01 g of ferric phosphate, and 1 L sea-water; pH 7.2-7.6) on a sterilized bench. The diluent was spread onto the 2216E plates, and the plates were cultured at 28 • C. The single colony, with agar collapse, was picked and streaked twice on the 2216E plates for further purification.
A colony of the WPAGA4 strain was purified using the above method. The genomic DNA of WPAGA4 was extracted using an EasyPure Bacteria Genomic DNA Kit (Transgen, Beijing, China) according to the instructions provided. Its purity and integrity were tested with Nanopore 2000 (Thermo, Waltham, MA, USA) and 1% agarose gel electrophoresis, respectively. The high-quality genomic DNA was used as the template for the 16S rRNA gene amplification with the primer pairs of 27F (5 -AGAGTTTGATCCTGGCTCAG-3 ) and 1492R (5 -TACGACTTAACCCCAATCGC-3 ) according to the PCR conditions stated in previous work [52]. The PCR product was sequenced by Sangon Biotech (Shanghai, China) Co., Ltd. The sequence of the 16S rRNA gene was aligned to the databases of NCBI (https://www.ncbi.nlm.nih.gov/, accessed on 21 July 2022) and Ezbiocloud (https: //www.ezbiocloud.net/, accessed on 21 July 2022).

The Determination and Optimization of the Agarase Activity in the Fermentation Broth
To optimize the agarose-degrading activity in the fermentation broth of WPAGA4, the strain cells were cultured in 50 mL liquid 2216E medium at 25 • C, 28 • C, 31 • C, 34 • C, and 37 • C for 24 h with an inoculation size of 1% and a rotation speed of 200 rpm for the 250 mL shake flask. Then, the culture broth was centrifugated at 6000 rpm for 30 min at 4 • C, and the supernatant was collected for the determination of agarase activity. The methods for the measurement and analysis were reported in our previous work [40]. To measure the substrate's inducibility, the cells were cultured at 31 • C for 24 h in the blank 2216E medium and the 2216E medium with an additional 0.2% agarose (Westbio, Spanish). Then, the supernatant was collected, and the enzymic activity was determined as mentioned above.
To determine the effects of the carbon and nitrogen sources on the supernatant activity of WPAGA4, 1% (w/v) of NH 4 NO 3 , KNO 3 , NH 4 Cl, (NH 4 ) 2 SO 4 , urea, corn steep liquor, soybean cake, glucose, sucrose, maltose, glycerin, amylum, and lactose were added to different samples of 2216E media used in the strain culture. The culture was performed at 31 • C for 24 h, and the supernatant was collected for agarase activity determination as mentioned above. All determinations were repeated three times.

The Analysis of Agarase Gene Aga3420 in the Genome of WPAGA4
An agarase gene, aga3420, was annotated in the genome of WPAGA4. The nucleoti-de sequence of aga3420 was aligned using BLAST software against the NCBI nr data-base for the similarity analysis. The analysis of the glycoside hydrolase family of the putative agarases was performed by the dbCAN2 meta server (https://bcb.unl.edu/dbCAN2/index. php, accessed on 5 January 2022) [65]. Phylogenetic analysis of the amino acid sequence of aga3420 was conducted using Mega 5.0 software with the neighborhood-joining method.

The Cloning and Expression of the Agarase Genes
Primer pairs in Table S4 were used for the amplification of the putative agarase genes aga3420 using a C1000 Thermal Cycler (Bio-Rad, California, CA, USA). The PCR program was as follows: 3 min at 95 • C; 25 cycles of 95 • C for 15 s; 55 • C for 15 s; 72 • C for 1 min/kb; and an extension at 72 • C for 6 min. The PCR products were further purified using a Universal DNA Purification Kit (Tiangen Biotech, Beijing, China). The purified PCR products were ligated with the pEAZY ® -Blunt E2 expression vector (TransGen, Beijing, China). The recombinant vector was transformed into the competent cells of E. coli BL21(DE3) pLysS.
The E. coli BL21(DE3) pLysS cells harboring the recombinant vector were cultured at 37 • C in LB broth (5.0 g yeast extracts, 10.0 g peptone, 10.0 g NaCl, and 1 L deionized water) with additional 100 µg/mL ampicillin. The culture broth was induced by 0.1 mM IPTG at 16 • C for 18 h until the OD600 value reached 0.6-0.8. The cells were collected by centrifugation at 10,000 rpm for 10 min at 4 • C, and the cells were fractured using a sonicator (Scientz, Ningbo, China). The Ni-nitrilotriacetic acid (NTA) Sefinose™ Resin (Tiangen Biotech, China) was used to further purify rAga3420 with 6 × His-tagged from the crude enzyme solution. The molecular weight and purity of rAgaW1540 were subsequently estimated using SDS-PAGE and Coomassie brilliant blue staining.

Activity Determination and Characterization of rAga3420
The protein concentration was determined using the Bradford method according to the former study [66]. The agarase activity (U) was defined as the amount of the enzyme that released 1 µmol of reducing sugars per minute under certain conditions. The measurement method for the agarase activity was conducted according to the previous study [16]. Briefly, 10 µL of the enzyme solution (20.7 µg/mL) was mixed with 240 µL agarose (0.4%; w/v) dissolved in different pH buffers, including sodium citrate-dibasic sodium phosphate (50 mM; pH 3.0-8.0), Tris-HCl (50 mM; pH 8.0-9.0), and glycine-NaOH (50 mM; pH 9.0-10.6), and the reaction mixture was incubated at temperatures of 0 • C, 10 • C, 20 • C, 30 • C, 40 • C, 50 • C, 60 • C, 70 • C, and 80 • C for 20 min. Then, 750 µL 3,5-dinitrosalicylic acid (DNS) solution was added to the reaction mixture, and the reducing sugar was measured [67]. To determine the thermostability of rAga3420, the enzyme was pre-incubated at 10 • C, 20 • C, and 30 • C for 60 min, and the relative activity of rAga3420 was determined every 10 min. The effects of metal ions and chemical agents on rAga3420 were measured by determining the enzymatic activity in the solution containing 1 mM of Mg 2+ , Fe 2+ , Fe 3+ , Ca 2+ , Mn 2+ , Cu 2+ , Na + , K + , Zn 2+ , EDTA, SDS, DTT, and β-mercaptoethanol. The impact of salinity was determined by measuring the relative activity with 0.5-3.5 mmol/L of NaCl. The kinetic parameters of rAga3420 were measured by determining the activity within the agarose at final concentrations of 0.1-0.6 mg/mL, and the Km and Vmax values were calculated by the Lineweaver-Burk plot. All determinations were repeated three times. The relative activity was defined as the percentage of activity determined with respect to the maximum agarase activity.

The Determination of Degradation Products
The degradation products of agarose by rAga3420 were identified by liquid chromatography mass spectrometry (LC-MS; Agilent XEVO G2-XS QTOF, Palo Alto, PA, USA). According to the former work [68], the degradation solution was filtered with a filter membrane of 0.22 µm pore size and then mixed with a triple volume of absolute ethanol for removing the protein and undegraded polysaccharides. Subsequently, a C18 column (2.1 × 100 mm, 1.8 µm; Waters, Milford, MA, USA) was used for the separation of the oligosaccharides by gradient elution with 0.1% formic acid-water and formic acidacetonitrile at a flow rate of 0.5 mL/min. The injection volume was 1 µL. The values of m/z from 50 to 3000 were analyzed by a quadrupole mass analyzer with positive-ion ESI mode. The interface voltage, detector voltage, curved desolvation line temperature, and heat-block temperature were set as 4.5 kV, 1.65 kV, 200 • C, and 200 • C, respectively.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/md20110692/s1, Figure S1: The dbCAN annotation result of aga3420, Figure S2: Lineweaver-Burk plot of the kinetic parameter determination of rAga3420. Table S1: Comparison of the properties of rAga3420 and other GH50 agarases, Table S2: The values of ANI and DDH of Vibrio natriegens WPAGA4 with other Vibrio species, Table S3: The dbCAN annotation result of the 4 putative agarase genes, Table S4: The primer pairs for the putative agarase gene amplification.

Data Availability Statement:
The genome sequence of WPAGA4 was deposited into the GenBank database under the accession numbers CP094880 and CP094881.