Role of Epitranscriptomic and Epigenetic Modifications during the Lytic and Latent Phases of Herpesvirus Infections

Herpesviruses are double-stranded DNA viruses occurring at a high prevalence in the human population and are responsible for a wide array of clinical manifestations and diseases, from mild to severe. These viruses are classified in three subfamilies (Alpha-, Beta- and Gammaherpesvirinae), with eight members currently known to infect humans. Importantly, all herpesviruses can establish lifelong latent infections with symptomatic or asymptomatic lytic reactivations. Accumulating evidence suggest that chemical modifications of viral RNA and DNA during the lytic and latent phases of the infections caused by these viruses, are likely to play relevant roles in key aspects of the life cycle of these viruses by modulating and regulating their replication, establishment of latency and evasion of the host antiviral response. Here, we review and discuss current evidence regarding epitranscriptomic and epigenetic modifications of herpesviruses and how these can influence their life cycles. While epitranscriptomic modifications such as m6A are the most studied to date and relate to positive effects over the replication of herpesviruses, epigenetic modifications of the viral genome are generally associated with defense mechanisms of the host cells to suppress viral gene transcription. However, herpesviruses can modulate these modifications to their own benefit to persist in the host, undergo latency and sporadically reactivate.

Introduction
Herpesviruses are double-stranded DNA viruses with genomes that range between 152-172 kbp in length. These viruses are members of the Herpesviridae family, which encompasses at least 100 different viruses [1]. From these, eight herpesviruses have been shown to affect humans, which are distributed within three subfamilies: Alpha-, Betaand Gammaherpesvirinae [2]. While herpes simplex viruses type 1 and 2 (HSV-1 and HSV-2, respectively) and the varicella-zoster virus (VZV) are classified within the Alphaherpesvirinae subfamily, cytomegalovirus (CMV) and the human herpesviruses 6 and 7 (HHV-6, HHV-7) belong to the Betaherpesvirinae subfamily, and Epstein Barr virus (EBV) and the Kaposi's sarcoma-associated herpesvirus (KSHV) are members of the Gammaherpesvirinae subfamily [2]. All these viruses can infect humans asymptomatically or symptomatically with a variety of symptoms and clinical manifestations, leading in some circumstances to mild up to severe diseases. Other commonly studied herpesviruses that affect animals are pseudorabies virus (PRV) and Marek's disease virus (MDV), both belonging to the Alphaherpesvirinae subfamily [3,4].

Herpesviruses: Epidemiology and Illnesses
Human viruses of the Alphaherpesvirinae subfamily, such as HSV-1 and HSV-2, are highly prevalent worldwide with a prevalence estimated at 66.6% and 13.2%, respectively [5]. Data gathered in Asia indicate that these viruses have a frequency of 50% in children and 76.5% in adults [6], while in other countries lower values have been reported, such as 47.1% for HSV-1 and 12.1% for HSV-2 in the USA [7]. Diseases generated by HSV-1 and HSV-2 are usually mucocutaneous or affect the central nervous system [8]. The main transmission route of HSVs is through oral, oro-sexual, or sexual contact [8]. During early stages of infection, viral glycoproteins interact with their ligands on the cell surface for viral binding. HSV virions are overall composed by a lipid bilayer envelope decorated with numerous proteins and glycoproteins on the outer surface. Beneath the viral envelope, there are approximately 20 viral proteins, in a zone called the tegument, which harbors several viral proteins that play essential roles over the host antiviral response early on after infection [9]. Below the tegument is the viral capsid which harbors the viral genome in a lineal form [9]. A similar structure is found in other herpesviruses, although with particularities for each of them. At least five glycoproteins are critical for HSV-1 entry and four for HSV-2 [10]. Glycoprotein B (gB) is key for virus attachment to the heparan sulfates on the host cell for both HSVs, while gC is also important for HSV-1, but not for HSV-2 [11]. gD is also present on the virion surface, and binds to any of the receptors mentioned above, namely nectin-1, nectin-2, the herpes virus entry mediator (HVEM), or 3-OS heparan sulfates [10,12]. Upon binding to one of its ligands gD activates two other essential glycoproteins on the virion surface, gH and gL, which act as a complex (gH/gL) that in turn activates the fusogenic properties of gB [10]. Once these interactions are established, HSV entry into the cytoplasm is mediated by membrane fusion, endocytosis or phagocytosis [13,14]. The viral capsid then enters the cell and is transported through microtubules to the nuclear pore, where the genome is injected into the nucleus and viral gene transcription initiates [15]. At this stage, productive viral infection begins with the transcription of a particular set of genes named immediate early genes (IE or alpha) [16]. A pivotal transcription factor called ICP4 (infected cell protein 4) is transcribed and translated early on (referred to as an immediate early gene; IE or alpha gene), which allows the beginning of a downstream sequential gene expression cascade that promotes the transcription and translation of early genes (E or beta genes), such as ICP8 (infected cell protein 8) or TK (thymidine kinase), and finally late genes (L or gamma genes), such as VP16 (viral protein 16) or gD (glycoprotein D) [16]. The latter are sometimes further classified as early late (gamma-1) and late (gamma-2) genes [17]. The classification of IE, E and L genes is based on the timing of their transcription during infection and dependency on different viral transcription factors [14,18].
Varicella-zoster virus (VZV), another alphaherpesvirus, causes varicella (also known as chickenpox) and zoster (also known as shingles) [19]. Current evidence suggests that primary VZV infection begins with the replication of the virus in respiratory epithelial cells [20]. Similar to other alphaherpesviruses, VZV binds to heparan sulfate receptors on the cell surface for its entry into the cell. Moreover, VZV can use the insulin-degrading enzyme (IDE) as the cellular receptor through an interaction that is mediated by the viral glycoprotein E (gE) in the MeWo cell line, which is a model of fibroblast. Once the virus is inside the cell, the expression of IE, E and L genes is initiated [21]. Infection is followed by a widely distributed vesicular rash that is a typical trait of varicella after the incubation period (approximately 10 to 21 days) [19]. This pattern likely reflects viral spread from the tonsils and other local lymphoid tissues, from where infected T cells can transport the virus via the bloodstream to the skin [22]. VZV is found in a worldwide geographic distribution, but annual epidemics are more prevalent in temperate climates, occurring most often during late winter and spring [23]. The prevalence of VZV infection in the population is elevated and, in the USA, has been reported to be over 93% in the 6-to 19-year old group and 98% in adults aged 20-49 years [24,25]. Cases of herpes zoster provide a source of VZV transmission to susceptible close contacts, causing varicella; the virus then spreads rapidly to other susceptible individuals, in part because, in contrast to other herpesviruses, VZV is transmissible by the respiratory route [19]. Varicella attack rates among susceptible household contacts exposed to VZV are approximately 90%; more limited exposures, such as those occurring in school classrooms, result in transmission rates of about 10 to 35% [26].
Two other alphaherpesviruses of interest are pseudorabies virus (PrV) and Marek's disease virus (MDV), which can infect animals. PrV is a neurotropic virus that causes Aujeskzy's disease in swine, although a broad spectrum of mammals can be infected with this virus, such as rodents, cats, dogs and cattle [27]. The seroprevalence of PrV in pigs in China has been reported to be approximately 35% [28]. Similar to HSV-1, neuron infection with PrV occurs through trigeminal neuron infection involving the interactions of the viral glycoprotein gD with the host nectin-1 receptor, which activates cytoskeleton remodeling, allowing viral particle entry into the cell [29]. MDV infection on other hand causes one of the most prevalent cancers described in the animal kingdom and mainly poultry [30]. The prevalence of infection can reach up to 49.5% in chicken, as reported in Iraq [31]. MDV infection is mediated by the interaction of the viral glycoprotein gC with glycosaminoglycans (GAGs) on cell surfaces and this virus mainly infects and replicates in immune cells, such as macrophages, natural killer (NK cells), as well as B cells, and establishes latency in T cells [30].
Betaherpesviruses such as CMV have a high prevalence in the infant population, with approximately 59% of children older than six years old having been exposed to the virus [32]. Strong differences in the epidemiology of CMV amongst the adult population in industrialized and emerging countries exist, where the latter have the highest frequency of infection reaching 100%, while industrialized countries display a range of infection between 60-70% [33]. CMV transmission occurs by direct contact with infected body fluids and there is evidence regarding vertical and breast-milk transmission [33]. The main cells infected by CMV are epithelial cells, fibroblasts and myeloid cells [34]. The mechanism of entry described for CMV involves molecular interactions between viral glycoprotein complexes, such as gH/gL and gH/gL/gp42, with integrins on the surface of epithelial cells and major histocompatibility complex class II (MHC-II) molecules in B cells [35]. Symptoms associated with CMV infections vary significantly, from asymptomatic to severe illness, with mononucleosis being the most common clinical presentation [36].
Another betaherpesvirus that can infect humans is HHV-6 (also known as HHV-6A). This virus infects children at a high rate and 90% of the human population is likely to being infected before the age of three [37]. The general route of infection is by contact with contaminated saliva, after which the virus enters mainly into CD4 + T cells through the CD46 receptor [38]. Symptoms after HHV-6 infection depend on the immunological state of the host. Infection in immunosuppressed individuals is related with the development of ataxia, hypersomnia, dementia, encephalitis, and myocarditis, among others [37]. The human herpesvirus 7 (HHV-7, also known as HHV-6B), is another betaherpesvirus able to infect humans and is highly prevalent worldwide [39]. In the United States of America, the prevalence of HHV-7 is above 85% [40]. However, prevalence rates can significantly vary between ethnicities [40]. This virus enters to the target cells through the CD134 receptor which is expressed only in activated CD4 + T cells [41]. Usually, HHV-7-related symptoms after infection only occur in immunosuppressed individuals and, in these subjects, they are somewhat similar to those reported for HHV-6 [39].
Epstein Barr virus and the Kaposi's sarcoma-associated herpesvirus (KSHV, or HHV-8) are gammaherpesviruses that can infect humans [42]. EBV infection in the population is widespread, over 90%, with asymptomatic infections mainly occurring during childhood [43]. EBV mainly infects B cells through the interaction between the viral glycoprotein-350 (gp-350) and the B cell receptor CD21 [44]. If infection occurs during adulthood, it can result in infectious mononucleosis, a self-limiting disease characterized by the inflammation of lymph nodes in the neck region [43]. Transmission of EBV is mainly through oral secretions [42], although infections after blood transfusions and organ transplants also occur [42]. EBV is known to be an oncogenic virus, as its infection is associated mainly with Hodgkin's lymphoma [45]. Additionally, EBV reactivation in immunosuppressed individuals may cause lymphoproliferative diseases [45]. On the other hand, recent evidence links EBV with the development of multiple sclerosis [46].
As its name suggests, KSHV is responsible for the development of all Kaposi's sarcomas [47]. Similar to other types of herpesviruses, KSHV interacts with heparan sulfates on the cell surface to enter the cell through the use of its viral glycoprotein gB [48]. In children, cases of hemophagocytic lymphohistiocytosis are described after infection with this virus [49]. The epidemiology of KSHV significantly depends on the geographical region [50]. While the prevalence of this virus in North America, Europe and Asia is lower than 5%, in Eastern Europe, the Middle East and the Caribbean the seroprevalence is above 50% [50]. Noteworthily, the highest seroprevalence of KSHV is found in Africa and regions of the Amazonas in Brazil, with values reaching 60% [50]. While the transmission route for this virus is not totally clear, its high prevalence in men that have sex with men supports the hypothesis of sexual transmission [50]. However, there is a high level of seroprevalence in children in Africa ,also suggesting salivary transmission [50]. Vertical and blood transfusion transmission are also well described [51]. Importantly, this virus is able to infect numerous cell types, such as dendritic cells (DCs), monocytes, B lymphocytes and oral epithelial cells [48,51].
Overall, herpesviruses are present in humans at a high prevalence, with the exception of HHV-8, which has a more limited penetration in the population. Importantly, these viruses elicit lifelong infections through the establishment of persistent infections and can undergo latency. Lifelong infections and frequent reactivations for some of these viruses evidences the existence of immune evasion mechanisms that allow them to persist in the host. Although there are numerous antivirals available that limit severe diseases caused by these viruses, they do not clear the virus from the host and thus, reactivations may occur during the entire life of an infected individual. Furthermore, some clinical manifestations elicited by these viruses, such as skin lesions produced by HSV-1 and HSV-2, are only mildly affected by such drugs, reducing in only a few days what are usually week-long lesions [52]. Therefore, a better understanding of the molecular determinants that favor herpesvirus infections and their persistence in the host is needed for developing better therapeutic approaches. In this review we revise and discuss common epitranscriptomic and epigenetic modifications reported in virus-derived nucleic acids and how these are modulated by host and viral determinants.

Epitranscriptomic Modifications and Herpesviruses N 6 -methyl-6-adenosine (m 6 A):
The N 6 -methyladenosine (m 6 A) modification of RNA molecules consists of the addition of a methyl group at position N 6 of the adenosine nucleotide and is the most prevalent internal modification occurring in mRNA [53]. The reaction is catalyzed by enzymes called "writers" and takes place co-transcriptionally in the nucleus [54]. Deposition of m 6 A onto mRNA molecules is catalyzed by a methyltransferase complex formed by METTL3 and METTL14, which uses S-adenosylmethionine (SAM) as a substrate [55]. Additionally, METTL3/METTL14 complex is associated with the protein WTAP, which is important for substrate targeting [55]. On the other hand, a set of "reader" proteins recognize m 6 A-modified mRNA molecules to exert the function of the modification [56]. The best characterized m 6 A reader proteins belong to the YTH-domain family (YTH) and are associated with the regulation of different steps of mRNA metabolism ranging from alternative splicing, nuclear export, translation efficiency, to degradation and subcellular localization [57][58][59][60]. Lastly, other enzymes that interact with RNAs having m 6 A modifications are the "erasers" such as alkb homolog 5 RNA demethylase (ALKBH5) and alpha-ketoglutarate-dependent dioxygenase (FTO), which can reverse this modification through demethylation reactions [54].
Importantly, the m 6 A modification in RNA has been widely described in the context of viral infections amid both RNA and DNA viruses, including SARS-CoV-2, HIV-1, Zika, dengue, hepatitis C virus (HCV) and HSV-1, among others [61][62][63][64][65][66]. The role of the m 6 A RNA modification has been related to different processes occurring during the viral replication cycle, including viral gene expression and assembly of new virions, but also immune evasion [67].

Pseudouridine:
Pseudouridine (Ψ) is the most abundant internal RNA modification found in transfer RNA (tRNA) and ribosomal RNA (rRNA) [68]. This modification is irreversible and occurs through two different pathways, one that is RNA-dependent and another that is RNA-independent. RNA-dependent mechanisms involve a family of ribonucleoproteins called H/ACA small nuclear ribonucleoproteins (sRNPs) that use a guide RNA (gRNA) complementary to the sequences that will be modified. This gRNA hybridizes with target sequences and adopts a secondary structure so that the sRNPs recognize this structure and can catalyze the isomerization of uracil into pseudouridine [68]. RNA-independent mechanisms involve pseudouridine synthase enzymes (PUS) [68], which can recognize secondary structures within RNA molecules that contain uridine targets [68]. These enzymes are known to modify tRNA (PUS 2/3/4/6/9), small nucleolar RNA (snRNA) (PUS 1/7), rRNA (PUS 5/7) and mRNA (PUS 1/2/3/4/6/7/9). Interestingly, reader proteins related to this modification have not been identified to date [69].
Importantly, pseudouridylation is known to occur in both tRNA and rRNA, yet its role in mRNA is unclear to date [70]. Current evidence suggests that pre-mRNA pseudouridylation can modulate alternative splicing through diverse mechanisms, such as the modification of the affinity of RNA-binding proteins to RNA, or by modulating secondary structures within pre-mRNAs [70]. Pseudouridylation might stabilize RNA, as it can modify its structure by eliciting base-pairings that promote a rigid backbone [71]. Additionally, there is evidence suggesting that pseudouridylation confers RNA resistance to degradation, as there are studies reporting that pseudouridine-modified RNAs are more resistant to snake venom and spleen phosphodiesterases, as compared to RNAs with non-modified uridine [72]. Another process in which pseudouridylation is involved is during translation, as pseudouridylation in mRNAs can be recognized by the MetRS (which is an tRNA synthetase of methionine specifically), for its translation with regulatory purposes [73]. This may reduce the binding of antisense RNA molecules improving translation efficiency [73]. An interesting finding regarding pseudouridylation is its role in immune escape. Recent evidence suggests that a complete replacement of uridine with pseudouridylation dramatically decreases the levels of IFN-β mRNA due to RNA containing pseudouridylation, failing to activate conformational changes in retinoic acidinducible gene I protein (RIG-I), which is directly related to the activation of the antiviral interferon pathway [74][75][76][77]. While RNA molecules carrying uridine modifications can activate some Toll-like receptors, such as TLR3, to switch on host immune responses [76], pseudouridylation is not recognized by such TLRs due to steric incompatibility [76].
Importantly, pseudouridylation modifications have been described in viral RNAs, such as HIV-1 and SARS-CoV-2 [78,79]. For example, in SARS-CoV-2 pseudouridylation was detected in mRNA, yet its role is unclear and the question of whether pseudouridylation may be associated with immune escape remains open [79]. Taken together, more studies are needed to better understand the effect and role of pseudouridylation in the context of viral infections. Interestingly, pseudouridylation modifications have been used in the development of mRNA-based vaccines due to its characteristics that are related to a decreased activation of anti-RNA innate immune system components [72] and increased production of the antigen [72].

Methyl-5-cytosine (m 5 C):
The 5-methylcytosine modification in RNA molecules consists in the addition of a methyl group to the carbon in position 5 of cytosine, which is catalyzed by methyltransferase enzymes. To date, NOL1/NOP2/SUN domain (NSUN) proteins (NSUN 1-7) have been shown to introduce such a methyl group in mRNAs [80]. Interestingly, m 5 C has been identified in mRNA molecules isolated from diverse types of viruses, such as SARS-CoV-2 and HIV-1, as well as in non-coding RNA of the herpesvirus EBV [81][82][83]. These modifications have been reported to be related with both the function and stability of mRNA molecules and have been found in at least 1.4% of transcripts derived from retroviruses. Additionally, occurrence is ten-times higher in these mRNAs than in cellular mRNAs during infection [82]. Noteworthily, m 5 C has been associated with an increase in the translation of HIV-1 mRNAs [84].

-O-Methylation (-O-me):
The 2 -O-methylation modification in RNA molecules consists in the addition of a methyl group in the 2 -OH of the ribose [85]. This modification has been identified in rRNA and tRNA molecules and evidence suggests that it can also be found in mRNA transcripts [85]. Importantly, all nucleotides, either canonical or non-canonical, can be 2 -O-methylated by sets of enzymes called MTases [85].
In the cellular context, 2 -O-Methylation is mainly associated with RNA stability and translational modulation [85]. A role for this modification in viral infections has been described for the West Nile virus, in which the 2 -O-methylation of the 5 cap enhances its virulence in macrophages, because this modification allows the viral evasion of intrinsic cellular defense mechanisms, such as interferon-induced proteins with tetratricopeptide repeats (IFIT), which are associated with antiviral activity. These proteins inhibit translation by interactions between their tertiary structures and partner proteins used by viruses to translate their proteins, such elongation factors eIF2, eIF3 and eIF4 [86,87]. Cap 2 -O-methylation has been described during SARS-CoV-2 infection, with this modification being performed by the nonstructural proteins NSP10/14 and NSP10/16 heterodimers. Interestingly, this allows the virus to evade host nucleases and protects mRNA from degradation, providing the virus additional immune escape capabilities [88,89]. Moreover, internal 2 -O-methylation was shown to avoid immune sensing of the HIV-1 RNA in dendritic cells and macrophages [90].
Overall, RNA modifications in the context of viral infections are diverse, with multiple different nucleotides being affected and different effects over the stability and functions of RNA molecules (Table 1). For instance, some of these modifications can enhance viral mRNA stability and prevent degradation or enhance the efficiency of mRNA translation. Moreover, these modifications can be helpful for the virus for evading the host antiviral immune response by affecting the recognition of viral RNAs by Toll-like receptors and RIG-I-like receptors (RLRs), therefore dampening the activation of the cellular antiviral signaling pathways. Whether these modifications occur or not in the context of viral infections with herpesviruses and their roles are discussed below.

Episodic Splicing Modifications
Constitutive splicing is an RNA maturation step that consists of the removal of introns and the ligation of exons within pre-mRNA molecules. This process is mediated by small nuclear ribonucleoproteins (i.e., U1, U2, U3, U4, U5 and U6), which form a complex known as the spliceosome [91]. Importantly, the cleavage of introns and the ligation of exons can vary depending on the RNA molecule and eventually generate different types of mRNA molecules starting from the same transcript in a process called alternative splicing [91].

Alternative Splicing Events
Alternative splicing can be generated by the inclusion or exclusion of introns and exons in the final mRNA transcript [92]. Among alternative splicing events, exon skipping is the most common, which consists of the inclusion or exclusion of an exon into the mature mRNA molecule. While the exact mechanism of exon skipping is somewhat unclear, there is evidence that relates this process with nonsense mutations, although additional investigations are needed [92]. Another alternative splicing process consists of mutually exclusive exons in which the splicing of exons is performed in a coordinated manner, where one exon is retained and another is spliced-out. Finally, intron retention occurs when a complete and unspliced intron is retained within the mature mRNA [91]. Other methods generating alternative mRNA transcripts consist of the use of alternative transcription start sites (TSS), as displacing the transcription initiate site can lead to the exclusion of exons [92].

Splicing Modifications and Herpesvirus Infection
For HSV-1, very few genes have been described to undergo splicing. However, there is evidence for a role of the viral protein ICP27 in regulating alternative splicing during infection. Indeed, ICP27 has been reported to inhibit alternative splicing during lytic infection of human epithelial kidney cells (HEK-293 cells) [93]. Interestingly, the inhibition of alternative splicing has been described to be mediated through the interaction of ICP27 with snRNPs involved in this process, which were redistributed [91]. Conversely, when ICP27 is absent during the latent phase of HSV-1, viral transcripts undergo splicing to restrict viral antigen expression and evade the host immune response [93]. Alternative splicing has also been described for MDV [94]. In this case, the ICP27 protein of MDV can interact with splicing regulator proteins (SR) in lymphoid cells (MSB-1 cells, which is the lymphoblastoid cell line) [91]. However, as opposed to what has been described with HSV-1, the role of MDV ICP27 in alternative splicing of viral transcripts has not been reported, and thus it remains to be clarified if this protein has any impact over the virus' lytic or latent phases. Regarding PRV, alternative splicing events have been described in mRNA transcripts of the epo and ul21 genes in epithelial cells (i.e., PK15 cell line) [95,96]. However, similar to MDV, the role of alternative splicing events during the latent or lytic phases of this virus remains an important question to address experimentally.
For the betaherpesvirus CMV, an interesting role for alternative splicing has been described, because alternative splicing in the major immediate early gene locus (MIE), particularly in exon 4, generates an abortive infection in human foreskin fibroblast cells. Therefore, alternative splicing of MIE during the lytic phase could be detrimental for viral replication [97].
Regarding the gammaherpesvirus EBV, there is evidence indicating that EBNA-2 and EBNA-LP viral proteins can modulate alternative splicing and regulate the expression of genes related with the survival and proliferation of infected B cells, as B cell lymphoma protein (Bcl-L) and Numb endocytic adaptor protein (NUMB); both are regulated by alternative splicing [98]. For KSHV, the alternative splicing of viral transcripts has been described for ORF57, which promotes the expression of lytic genes and can regulate the splicing of the biscistronic ORF70/K3 RNA, which is associated with viral DNA replication in the lytic phase of this virus and immune evasion [99].
Taken together, alternative splicing of viral RNAs in herpesvirus has been reported and relates to immune escape and optimal viral lytic phases, although further studies are needed to fully grasp all the viral and host factors involved and how they are modulated during both the lytic and latent phases of viral infection.

Epitranscriptomic Modifications in Alphaherpesvirus Infections
Overall, the mechanisms by which m 6 A contributes to viral infection with HSV-1 are poorly understood and evidence is contrasting [65]. At present, there are studies that report that m 6 A modifications do exist in HSV-1 transcripts in host cells, such as human epithelial cells (HeLa), and this modification is critical for viral replication [65]. Interestingly, infection by HSV-1 increases the expression of m 6 A writers METTL3 and METTL14 in infected cells [65]. Moreover, silencing of these proteins results in a significant reduction in viral replication [65]. Consistently, an enhancement in viral replication and replication was observed when the eraser enzyme ALKBH5 was knocked-down [65]. However, recent evidence suggests that the alphaherpesvirus protein US3 of the pseudorabies virus (PRV), which has kinase activity, can inactivate the m 6 A writer complex METTL3/METTL4/WTAP through phosphorylation in a fibroblast cell line (ST), producing a decrease in m 6 A levels and evidencing that m 6 A seems not to be important for PRV replication [100].Therefore, effects of m 6 A in viral replication may be specific for the cellular type infected and the alphaherpesvirus evaluated. Noteworthy, it has been observed that upon HSV-1 infection of fibroblasts, the m 6 A writer METTL3 is re-distributed from the nucleus to the cytoplasm. Importantly, these changes in the distribution of this protein were mediated by the HSV-1 protein ICP27, particularly by the carboxy-terminus domain of this protein which has both structural and sequence homology with regulators of RNA processing and export [101].
Regarding m 6 A reader proteins, these have also been found to be important in the replication of HSV-1 by unknown mechanisms. Given the generally described effects of m 6 A over the increase in stability of mRNAs and translation efficiency, HSV-1 may recruit writer proteins to its modified mRNAs to improve their stability and translation efficiency, which would likely produce enhanced viral replication [65,101].
For other alphaherpesviruses, such as MDV, infection has been reported to modulate the m 6 A landscape in host cells, mainly by increasing the frequency of the m 6 A modification in long non-coding RNAs (lncRNAs) in a chicken fibroblast cell line [102]. However, additional evidence is needed to understand how these changes could affect the replication of this virus.
Noteworthily, the effects of m 5 C, Ψ and 2 -O-methylations in the context of alphaherpesvirus infections have not been reported to date ( Figure 1). Nevertheless, during recent years new technologies, such as direct RNA sequencing, have been emerging, which may help to elucidate any effects these modifications may have in the replicative cycle of alphaherpesviruses relevant to human health, such as HSV-1, HSV-2 and VZV [103].
Given these methodological advances, we foresee in the near future the identification and characterization of pseudouridylations and 2-O-methylations within alphaherpesvirus infections.

Epitranscriptomic Modifications in Betaherpesvirus Infections
CMV infections have been associated with an increase in reader proteins, such as YTHDF2 in NK cells, and inhibiting this protein in vivo reduced the host's capacity to control this virus [104]. Interestingly, CMV infection of NK cells affects proteins that read or erase the m 6 A modification in the host cell. For example, an increase of the YTHDF2 reader protein was observed 1.5 days post infection, which was accompanied by a decrease in YTHDF3. Increased YTHDF2 expression was observed in response to viral infection and was associated with the activation of the transcription-factor signal transducer activation of transcription 5 (STAT-5), which activates host antiviral pathways mediated by IFN-γ [104]. However, this response may vary depending on the cell type, as an increase in the expression of the writer proteins METTL3, METTL14, together with the YTHDF2 reader, was observed after CMV infection in fibroblasts. Noteworthily, these enzymes were found to be necessary for virus propagation [105]. Indeed, when the writer protein METTL3 was inhibited in fibroblasts infected with CMV, an enhanced type-I interferon response was observed, suggesting that CMV may use this modification as an evasion mechanism [105].
Given these methodological advances, we foresee in the near future the identification and characterization of pseudouridylations and 2-O-methylations within alphaherpesvirus infections.

Epitranscriptomic Modifications in Betaherpesvirus Infections
CMV infections have been associated with an increase in reader proteins, such as YTHDF2 in NK cells, and inhibiting this protein in vivo reduced the host's capacity to control this virus [104]. Interestingly, CMV infection of NK cells affects proteins that read or erase the m 6 A modification in the host cell. For example, an increase of the YTHDF2 reader protein was observed 1.5 days post infection, which was accompanied by a decrease in YTHDF3. Increased YTHDF2 expression was observed in response to viral infection and was associated with the activation of the transcription-factor signal transducer activation of transcription 5 (STAT-5), which activates host antiviral pathways mediated by IFN-γ [104]. However, this response may vary depending on the cell type, as an increase in the expression of the writer proteins METTL3, METTL14, together with the YTHDF2 reader, was observed after CMV infection in fibroblasts. Noteworthily, these enzymes were found to be necessary for virus propagation [105]. Indeed, when the writer protein METTL3 was inhibited in fibroblasts infected with CMV, an enhanced type-I interferon response was observed, suggesting that CMV may use this modification as an evasion mechanism [105].
Findings from another study propose that CMV uses the host m 6 A RNA modification machinery to modify the expression of host proteins during infection. Zu and colleagues. observed that CMV infection in endothelial cells elicits an aberrant and increased expression of writer enzymes, such as METTL3, and YTHDF2 reader proteins. This effect pro- Findings from another study propose that CMV uses the host m 6 A RNA modification machinery to modify the expression of host proteins during infection. Zu and colleagues. observed that CMV infection in endothelial cells elicits an aberrant and increased expression of writer enzymes, such as METTL3, and YTHDF2 reader proteins. This effect produced a modification of mRNA encoding host proteins, such as MCU (mitochondrial calcium uniporter), with METTL3 modifying the mRNA of this host factor in the 3 UTR region, which in turn was read by YTHDF2, resulting in an aberrant expression of MCU that facilitates apoptotic processes within the cell [106]. Similar to alphaherpesviruses, further evidence regarding the role of other RNA modifications in betaherpesvirus infections require studies yet to be carried out and reported (Figure 1).

Epitranscriptomic Modifications in Gammaherpesvirus Infections
The effects of the m 6 A modification have also been studied in EBV. Interestingly, the EBV early protein BZLF1 was shown to bind the promoter of the writer enzyme METTL3, negatively modulating its expression [107]. Moreover, Zheng et al. further described that during the pre-latency phase of EBV infection, the expression of the viral protein EBNA-2, which is involved in transcriptional regulation of viral transcripts and is associated with the immortalization of EBV-infected cells [108], decreased upon knock-down of METTL3 in B cell lymphoma cells (BJAB cells) [109]. An important observation reported in this work was that EBV infection can induce the hypomethylation of the mRNA of TLR9, affecting its stability, which may affect the initiation of an innate immune response, altogether favoring the virus' capacity to achieve long-term latency in the host [109]. Another study found that the deletion of the writer protein METTL3 in B cells negatively affected the expression of viral lytic proteins and decreased the yield of progeny virions, further suggesting that METTL3 may have a role in the replicative cycle of EBV [110].
Another modification that has been studied in the context of EBV infection is m 5 C. Some studies suggest that the EBV viral non-coding RNA of the EBER1 gene is modified with m 5 C [81], and that the writer enzyme that performs this modification is NSUN2, which is associated with a decrease in EBER1 levels [81,110]. However, the non-coding RNA EBER can bind to Toll-like receptors and activate the interferon pathway, hence the inhibition of EBER1 levels in the host cell, due to m 5 C modifications maybe acting as an immune escape mechanism evolved by this virus [81].
Another gammaherpesvirus undergoing the m 6 A modification in its transcripts is the Kaposi's sarcoma-associated herpesvirus (KHSV), whose viral mRNA is extensively methylated with m 6 A during the lytic and latent phases, with a marked increase in mRNA displaying the m 6 A modification when lytic replication is stimulated [111]. Importantly, enzymes involved in this modification are critical for KSHV infection, with evidence suggesting that the reader protein YTHDF2 suppresses KSHV infection. This suppressive effect has been reported to be mediated by the deadenylation and degradation of viral transcripts [88]. However, findings in this area are contrasting and effects of writer and reader proteins apparently depend on the cell type used [112]. For example, Hesser and colleagues. showed that YTHDF2 and METTL3 are necessary for virion production in Kaposi's sarcoma-derived cells (iSLK.BAC16 cells). Yet, the effects were different in other cells obtained from KS, such as iSLK.219 cells, in which YTHDF2 and METTL3 were associated with the accumulation of the viral lytic transactivator ORF50, and not in virion production [112]. Moreover, the non-YTH reader protein staphylococcal nuclease and tudor domain containing 1 (SND1) has been reported to be involved in the replication of KSHV and its disruption found to interfere with the replication cycle of KSHV in a cell line of B-cells (TREx BCBL1-Rta cells) due to the protective role of SND1 in ORF50 mRNA, stabilizing and protecting against degradation [113]. ORF50 is the first lytic protein produced by this virus and is necessary to switch between lytic and latent phases [113].
Overall, the effects of m 6 A during KSHV indicate that this modification is necessary for viral replication, and that virus infection per se can modify the landscape of viral and host m 6 A RNA modifications. However, similar to the infection with CMV, the reader protein YTHDF2 has also been involved with an antiviral activity; in this case specifically mediating viral transcript degradation. On the other hand, YTHDF3 would seem to have a positive role in viral replication through a mechanism yet to be determined. The influence of other RNA modifications during both the latent or lytic phases of this virus remains to be identified (Figure 1).

Histone Modifications
Some viruses containing double-stranded DNA as genomes can be packaged with host histones into the virion. Indeed, this has been described for HSV-1 and some adenoviruses, for example [114,115]. Histone modifications are catalyzed by enzymes that add or delete chemical groups (mainly acetyl and methyl groups) in lysine or arginine residues. For instance, enzymes adding an acetyl group are called histone acetyltransferases (HAT), while the opposite reaction is carried out by histone deacetyltransferases (HDACs) ( Table 2) [116]. Importantly, acetyl groups modify the electrostatic charges of histones affecting molecular interactions between the histones and DNA, altering the opening or closure of chromatin structures [116]. Differences in the quantity and types of acetyl groups added to histones are associated with the activation or repression of transcriptional activity [116]. For example, acetylation of lysine 27 in histone 3 (H3K27ac) is associated with active promoters and enhancers because this modification favors open chromatin, and therefore transcriptions factors can bind to promoters to start transcription [117]. While H4K20ac is associated with gene repression in human cells because it promotes chromatin closure [118].
Additionally, histones can also be methylated in a reaction catalyzed by histones methyltransferases (HMT), and different modifications have been described, such as the trimethylation of histone 3, lysine 27 (H3K27me3), which is associated with transcriptional repression, while methylation in H3K4me is a modification related to the transcriptional activity of genes [119]. Table 2. Viral DNA modifications by writer and eraser enzymes.

Writers Erasers Function
Viruses in Which Effects Have Been Described
Diverse modifications are involved in gene regulation.

5-methylcytosine (5mC)
Epigenetic modifications are responsible of phenotypic or gene expression changes that cannot be explained by alterations in DNA sequences [120]. Although these modifications are widely known to occur in the genome of host cells, DNA methylation, particularly 5-methylcytosine (5mC), has also been described in the context of infections with DNA viruses (Table 2) [121]. This modification consists in the addition of a methyl group at the carbon in position 5 of cytosines, in a reaction that is catalyzed by DNA methyl transferase enzymes (DNMTs) DNMT1 and DNMT3 [120]. This modification mainly occurs at CpG dinucleotides, also called "CpG islands" [120]. Importantly, methylation is a reversible modification and can be edited by demethylase enzymes called ten-eleven translocation methylcytosine dioxygenase (TET), which oxidize 5 methylcytosine to intermediary compounds that are finally converted into cytosine [126]. Noteworthily, 5mC can be recognized by methyl-binding proteins (MBPs), a type of reader protein, which can act as a "scaffold" protein between the modified DNA and repressor complexes, such as Rest Corepressor (CoREST), with gene silencing as a result [127].

Epigenetic Modifications in Alphaherpesvirus Infections
Although the incoming DNA genome of HSV-1 lacks 5mC, it is methylated by DNMT's of the host upon infection [122]. Methylation of the HSV-1 genome by DNMT3A, which decreases at later stages of the viral replication cycle, was found to be important for viral replication [122]. Indeed, knocking down of DNMT3A induced a dramatic decrease in viral titers in human fibroblast foreskin-infected cells, highlighting the critical relevance for an active viral replication [122]. Interestingly, during the latency phase of this virus in rat Schwann cells (RSC), CpG methylations were not identified, relating this modification with gene repression, and suggesting that infected cells might use this process to silence viral gene expression [121]. Therefore, methylation of viral DNA occurs during the lytic phase and subsequently seems to be deleted during the establishment of latency [121,128]. Thus, an interesting field for further investigations is how HSV-1 regulates the deletion of methylations during the switch of lytic and latent phases and how viral proteins can regulate the host machinery of DNA methylation during latency.
Another epigenetic modification studied in the context of HSV-1 infections is posttranslational modifications in histones. Evidence suggests that the alphaherpesvirus protein US3 from HSV-1, VZV and PrV can hyperphosphorylate histone deacetylase 2 (HDAC2) and that this effect can reduce viral genome silencing during latency phases to allow efficient viral replication upon the lytic phase [129]. Although the HSV-1 genome does not encode histone proteins, 1-h post infection the viral genome has been found to be associated with H3 histones and can actually be packaged into nucleosomes in neurons (SH-SY5Y neuroblastoma cells) [130]. These nucleosomes with H3 are similarly distributed amongst immediate early, early and late HSV-1 genes [130]. Interestingly, an increase in histone modifications that promote DNA availability, such as H3K27ac, were observed to be increased during HSV-1 replication, while repressive modifications, such as H3K27me3, were found to be decreased in THP-1 cells, a monocyte cell line [131]. Therefore, it is possible that the expression of viral genes is regulated by the accessibility of its chromatin [131]. Indeed, the permissive histone modification H3K27ac was found to be critical for viral replication, as its pharmacological inhibition associated with a significant repression of HSV-1 replication [128].
Other studies have found that during early stages of infection with HSV-1, the viral transcriptional co-activator Host cell factor 1 (HCF-1) interacts with different proteins and protein complexes that remodel chromatin, such as methyltransferases and demethylases that specifically modify histones [132,133]. Indeed, during infection with HSV-1 or VZV, HCF-1 was found to interact with remodeling proteins, such as lysine-specific demethylase 1A (LSD1) and mixed-lineage leukemia protein 1 (MLL1), to actively elicit modifications that promote immediate early gene transcription. Furthermore, an interaction between HCF-1 and the chaperone protein anti silencing factor-1 (ASF-1), which manages the flux of histones, was found to be critical for replication of this virus [132].
During latent phases of infection, lytic genes of HSV-1 are condensed into heterochromatin, repressing viral gene expression ( Figure 2). This repression has been reported to be associated with the recruitment of a repressor protein complex known as Polycomb repressive complex 1 (PRC1), by proto-oncogene ring finger protein (Bmi1), which negatively impacts the binding of positive regulators of chromatin and promotes modifications such as acetylation of histones by a mechanism unclear to date [134]. Hence, histone modifications associated to viral DNA can modulate the expression of viral transcripts, and during the latent phases these viruses are likely to have a critical role in inhibiting the transcription of viral RNAs [134]. This approach would allow the virus to inhibit host cell antiviral countermeasures, because the genome of HSV-1 is packaged with histones similar to the host DNA [134].
Similar to human herpesvirus, animal alphaherpesviruses, such as MDV, regulate the latent and lytic phases through modifications of the chromatin. Indeed, MDV represses its gene expression in T cells through repressive histone modifications, such as H3K27me3 and H3K9me3 nearby lytic genes in the viral genome [135].
Epigenetic modifications have also been described in other types of alphaherpesviruses, such as VZV. Interestingly, despite the fact that VZV has 61 CpG dinucleotides in its genome, these are not located in promoter regions and therefore, their methylation may not be related to transcriptional suppression. However, more evidence is needed to confirm this hypothesis as 5mC in the body of genes also impact gene expression. Noteworthily, it has been found that during latency after VZV infection, only some viral proteins are expressed, namely the proteins ORF21, ORF29, ORF62 and ORF63 [136]. Importantly, the regulated expression of these ORFs during latency is necessary, and regulation is mediated by epigenetic markers, such as H3K9ac, which are maintained during latency [136,137].
Overall, epigenetic modifications have been reported to regulate both the lytic and latent phases of alphaherpesvirus infections (Figure 2). For example, during latency, alphaherpesviruses can condensate their genomes into heterochromatin and only genes related with latency, such as LAT transcripts in HSV-1, and VLT in VZV, are accessible for transcription. Condensed chromatin during latency also allows virus evasion of the immune system, because double-stranded packaged genomes with histones are somewhat similar to host DNA, and because this DNA is episomal and is not producing mRNAs translated into proteins, antiviral responses may not be activated (Figure 2) [138]. lated with latency, such as LAT transcripts in HSV-1, and VLT in VZV, are accessible for transcription. Condensed chromatin during latency also allows virus evasion of the immune system, because double-stranded packaged genomes with histones are somewhat similar to host DNA, and because this DNA is episomal and is not producing mRNAs translated into proteins, antiviral responses may not be activated (Figure 2) [138]. The methylation of the DNA is associated with methyl binding proteins (MBP), a type of reader protein, acting as a "scaffold" between 5mC and repressor of complexes such as CoREST. Next, the eraser enzymes (TET) can erase the methylation and form intermediates that are ultimately converted to cytosine such as 5'-hydroxy-methyl-cytosine (5hmC). (b) The histones can also suffer epigenetic modifications by enzymes such as HAT, HDAC and HMT, that favor the latency phase, depending on the modification. Acetylation of lysine 27 in histone 3 (H3K27ac), is associated with active promoters and enhancers. H4K20ac is associated with gene repression. Methylation of lysine 27 in histone 3 is associated with transcriptional repression. H3K4me is related to transcriptional activity of genes.

Epigenetic Modifications in Betaherpesvirus Infections
DNA methylation has been reported in the genome of CMV, specifically within the promoters of immediate early genes being associated with gene repression [123]. Importantly, the pharmacological inhibition of DNA methylation by compounds such as azacytidine have been shown to improve viral replication [123]. This led to the notion that the methylation of viral DNA may be a host defense mechanism for suppressing the The methylation of the DNA is associated with methyl binding proteins (MBP), a type of reader protein, acting as a "scaffold" between 5mC and repressor of complexes such as CoREST. Next, the eraser enzymes (TET) can erase the methylation and form intermediates that are ultimately converted to cytosine such as 5'-hydroxy-methyl-cytosine (5hmC). (b) The histones can also suffer epigenetic modifications by enzymes such as HAT, HDAC and HMT, that favor the latency phase, depending on the modification. Acetylation of lysine 27 in histone 3 (H3K27ac), is associated with active promoters and enhancers. H4K20ac is associated with gene repression. Methylation of lysine 27 in histone 3 is associated with transcriptional repression. H3K4me is related to transcriptional activity of genes.

Epigenetic Modifications in Betaherpesvirus Infections
DNA methylation has been reported in the genome of CMV, specifically within the promoters of immediate early genes being associated with gene repression [123]. Importantly, the pharmacological inhibition of DNA methylation by compounds such as azacytidine have been shown to improve viral replication [123]. This led to the notion that the methylation of viral DNA may be a host defense mechanism for suppressing the transcription of viral genes in infected cells. Interestingly, there is evidence indicating that CMV infection can overall modify the methylation landscape in the infected cell, with a global hypomethylation observed within fibroblasts after infection [123]. Moreover, CMV infection was also found to change the distribution of DNMT1 and DNMT3b from the nucleus to the cytoplasm, with a "DNMT-free environment" in the nucleus being beneficial for virus replication, by evading host mechanisms intended at suppressing viral gene transcription through epigenetic modifications by these enzymes [123].
Other studies reported that CMV replication requires the viral genome to be available for transcription, needing to switch from heterochromatin to euchromatin during productive infection [139]. Accordingly, histone modifications have been described, wherein H3K4me2 has been related to replicating or post-replicative CMV genomes. Other posttranslational histone modifications, such as H3K9me2 and H3K9/14ac, have been associated with the viral genome being structured as heterochromatin (Figure 2) [139]. This packaging of the viral genome in heterochromatin and as an episomal form is associated with latency because condensed chromatin is not associated to the production of transcripts nor translation of viral proteins ( Figure 2). As mentioned, this would allow the virus to not be recognized by immune system components, avoiding the activation of the host antiviral responses during latency [139].

Epigenetic Modifications in Gammaherpesvirus Infections
Similar to the other herpesviruses discussed in this review, the EBV genome has also been found to be methylated in CpG islands, with this modification being associated with transcriptional repression and alternative splicing of viral transcripts [140]. For example, when the cp promoter (first promoter in the EBV genome) is methylated at its upstream sites, a smaller transcript is generated (EBNA 1) and no other EBNA genes are expressed. This gene is related with viral latency in B cells, and the inclusion or exclusion of different EBNA genes is conditioned if latency is established in vivo or in vitro [140,141].
Interestingly, an epigenetic mechanism that controls the transition of the EBV virus from a lytic to a latent phase in B cells has been described. During early stages of infection, the transcription factor BZLF1 of this virus is able to bind to 5mC in DNA domains known as meZRE for gene transcription regulation [142]. On the contrary, CpG methylation of these DNA domains was associated with the activation of the transcription of viral lytic genes. Furthermore, genes that did not have CpG methylated dinucleotides were not regulated by BZLF1 [142]. Therefore, DNA methylation of EBV genes may not always be synonymous with gene repression. Importantly, the viral tegument protein (BNRF1) has been reported to prevent histone modifications during early phases of infection and its regulation is indirectly mediated by interactions with the protein Death-associated protein (DAXX), which is associated with chromatin remodelers such as ATRX in an epithelial cell line Hep-2 [124]. The displacement of the interaction between DAXX and ATRX would prevent nucleosome formation with the viral genome, specifically at latent genes [124].
Interestingly, inhibition of the expression of BZLF1 has been reported to be needed to maintain EBV latency. In this case, it has been suggested that modifications over histones that limit access to the DNA are critical, particularly with the H3K27me3 and H4K20me3 modifications being involved due to BZLF1 silencing. Also similar to other herpesviruses, EBV uses nucleosome packaging of the viral DNA to maintain latency ( Figure 2) [142].
Finally, methylation of the genome of KSHV has been reported, although its methylation pattern did not change up to 6 days post infection in SLK cells. Like other herpesviruses, the genome of KSHV encodes protein factors that can interact with DNMTs enzymes, such as the viral protein LANA, which is involved in the induction and maintenance of latency [143]. This protein also can interact with methyl-binding proteins to suppress genes related to the lytic phase. As KSHV can interact with enzymes involved in epigenetic modifications, it may also modify the methylation pattern of the host cell through the interaction between viral proteins and the cellular methylation machinery, such as DNMT3a and methyl-binding proteins, such as MBP and MeCP. Noteworthily, LANA may use this machinery to suppress the expression of lytic genes [125].
Histone modifications patterns in KSHV infection have also been explored during the latent phase of this virus and its reactivation, wherein repression of lytic genes is reportedly mediated by enzymes of the Polycomb complex, specifically the methylase EZH2, which generates repressive marks in lytic genes to maintain its suppression for maintaining latency in the TRExBCBL1-RTA cell line, which is a cellular line of lymphoma [144]. During reactivation onto the lytic phase, the ORF59 protein has been found to interact with histone demethylases UTX and JMJD3, forming a protein complex with PAN RNA motifs as intermediaries [144]. This interaction allowed the loss of repression modifications in histones, such as H3K27me3, leading to increased expression of viral genes [143,145].

Concluding Remarks
To date there is no effective treatment for clearing herpesvirus infections or effectively blocking their reactivations in the host. In this review, we revised and discussed how different post-transcriptional RNA modifications are associated with increased or decreased viral fitness, which hints to potential future targets for the development of new drugs that could inhibit RNA-modifying enzymes. However, the development of drugs that can modify RNA or DNA within the host may be double-edged swords, because modifying the epigenome landscape could result in aberrant gene expression patterns in host cells that could result in the malignant transformation of cells.
While the m 6 A is so far the most studied post-transcriptional RNA modification in mRNA and viral RNA, with its presence being described in RNA molecules during alphaherpesvirus, betaherpesvirus and gammaherpesvirus infections, its effect on viral replication would benefit from further studies. Furthermore, a role for this modification will depend on the presence of reader proteins that bind to m 6 A, converting this modification into functional outputs. On the other hand, those reader proteins identified so far in the context of these viral infections, such as YTHDF2, could be modulated positively to have antiviral activity against some herpesviruses, such as beta-and gammaherpesviruses. Moreover, readers such as YTHDF3, which are associated with improved viral replication, could be targeted to dampen viral replication to block the replication cycle of these viruses. Nevertheless, in either case, it is important to consider the effects that such interventions may have over normal cellular functions and thus could be restricted to local applications, particularly over viral infections that act in tissues such as the skin. Still, whether viral proteins have the ability to bind m 6 A-modified RNA is completely unknown. Regarding writer enzymes, such as METTL3, for all viruses discussed in this review these proteins seem to be necessary for viral fitness and this expression is translated in impairments in viral fitness.
The development of high throughput technologies and new antibodies will hopefully allow for advancing this field, to better identify those modifications occurring during infection with herpesviruses and better understand their replication processes. Importantly, while studying RNA modifications is relevant during the infective replication cycle of these viruses, little is known on the role of RNA modifications in the establishment and maintenance of latency, as some viral proteins are not expressed during this process [140]. Naturally, in the context of latency, epigenetic modifications are more likely to play critical roles during this phase, and yet this area has been relatively poorly explored. Interestingly, epigenetic marks such as H3K36me3 were shown to drive m 6 A deposition co-transcriptionally [146]. Whether this crosstalk between epigenetic and epitranscriptomic marks exists during viral infections has never been explored and deserves further exploration.
Although to our knowledge herpesviruses do not encode nucleic-acid-modifying enzymes in their genomes able to write or erase chemical modifications, they encode proteins that can interact with the epigenetic machinery of the host. These interactions allow the regulation of gene expression at different stages of the life cycle of these viruses and, importantly, regulate latency in host cells. Altogether, histone modifications are related with gene regulation and the establishment of latency in host cells. Funding: This work was supported by ANID-Millennium Science Initiative Program-ICN09_016: Millennium Institute on Immunology and Immunotherapy (ICN09_016/ICN 2021_045; former P09/016-F) and FONDECYT grants #1190864 and #1190156 from ANID. This work was also supported by the Regional Government of Antofagasta through the Innovation Fund for Competitiveness FIC-R 2017 (BIP Code: 30488811-0). Data Availability Statement: Data sharing does not apply to this article as no new data were created or analyzed in this review.


s2fieldsofstudy: Biology Medicine

text:
Genome-scale metabolic modelling of the human gut microbiome reveals changes in the glyoxylate and dicarboxylate metabolism in metabolic disorders

Summary The human gut microbiome has been associated with metabolic disorders including obesity, type 2 diabetes, and atherosclerosis. Understanding the contribution of microbiome metabolic changes is important for elucidating the role of gut bacteria in regulating metabolism. We used available metagenomics data from these metabolic disorders, together with genome-scale metabolic modeling of key bacteria in the individual and community-level to investigate the mechanistic role of the gut microbiome in metabolic diseases. Modeling predicted increased levels of glutamate consumption along with the production of ammonia, arginine, and proline in gut bacteria common across the disorders. Abundance profiles and network-dependent analysis identified the enrichment of tartrate dehydrogenase in the disorders. Moreover, independent plasma metabolite levels showed associations between metabolites including proline and tyrosine and an increased tartrate metabolism in healthy obese individuals. We, therefore, propose that an increased tartrate metabolism could be a significant mediator of the microbiome metabolic changes in metabolic disorders.

INTRODUCTION
In recent years we have seen significant advances in elucidating the importance of the gut microbiome in human health and disease (Gentile and Weir, 2018;Zomorrodi and Maranas, 2012). Microbial communities have an intimate symbiotic relationship with their host, promoting protection against pathogenic microbes, maintenance of homeostasis, and processing of nutrients otherwise indigestible by humans . Shifts and alterations in the microbiome have been linked to different environmental factors including age, geography, and body mass index (BMI) (Le Chatelier et al., 2013b;Yatsunenko et al., 2012). In addition, diet has a significant impact on the microbiome (David et al., 2014). Various diseases have been associated with altered microbiome composition, and many studies have shown these communities can contribute to metabolic inflammation and metabolic disorders (Fan and Pedersen, 2021;Johnson et al., 2017;Proffitt et al., 2020;Tilg et al., 2020). There are several studies linking major metabolic disorders (type 2 diabetes (T2D), obesity, and atherosclerosis cardiovascular disease (ACVD)) with changes in gut microbiota composition (Karlsson et al., 2012(Karlsson et al., , 2013aLe Chatelier et al., 2013b). These dysbiotic changes resulting from a loss in stability of the gut microbiome drive changes in the gut ecosystem that result in a reduction of microbial diversity (Litvak et al., 2018). In turn, this results in shifts in key metabolite production. For example, short-chain fatty acids (SCFAs) such as acetate regulate immune cell production and help maintain intestinal homeostasis (Yang et al., 2020) and have been linked with promoting obesity and insulin resistance (Perry et al., 2016). Likewise, amino acid (AA) metabolism in the gut microbiome has been observed to have a large impact on health and disease (Cook and Sellin, 1998;Oluwagbemigun et al., 2020). In particular, plasma glutamate levels positively correlate with increased BMI and with fasting triglycerides, both of which can lead to insulin resistance (Palomo-Buitrago et al., 2019).
Despite a large number of gut metagenomic studies in metabolic diseases, the contribution of specific microbes to host metabolism during metabolic disease, systematic analysis, and modeling has not been mechanistically studied. In particular, little is known about the specific impact of microbe-microbe metabolic interactions. Recent advances in the metabolic modeling and reconstruction of genome-scale metabolic models (GEMs) have enabled the study of species-specific metabolisms and metabolic interactions within microbial communities (Gu et al., 2019;Magnusdottir et al., 2017;Shoaie et al., 2015). These computational models provide a platform that can be used to explore the genotype-phenotype relationships thereby allowing us to predict different phenotypic possibilities for microbes under different sets of constraints (Karlsson et al., 2011). This facilitates the study of interactions between species, and the quantification of growth, consumption, and production of metabolites. Using constraintbased modeling, diverse phenotypes under different conditions, such as nutrient availability can be simulated (Bordbar et al., 2014).
In this study, we focused on gut microbial metabolic modeling of three main metabolic disorders, ACVD (Jie et al., 2017;Karlsson et al., 2012), obesity (Le Chatelier et al., 2013b;Qin et al., 2010), and T2D (Karlsson et al., 2013a;Qin et al., 2012a). We used publicly available gut metagenomics datasets from these cohorts along with a previously generated list of significantly associated species in each study from our recent human gut microbiome atlas; Database: www.microbiomeatlas.com (Shoaie et al., 2021). These datasets were previously analyzed based on the recent gut integrated gene catalog (Wen et al., 2017), and the metagenome species profiled using our recently updated 1,989 Metagenomic Species Pan-genomes (MSPs) (Plaza Oñ ate et al., 2018). We applied GEMs for these species to investigate how gut microbiome metabolism varies between these three metabolic disorders. Individual and community-level modeling were performed together with reaction abundance to pinpoint specific bacterial metabolites and reactions associated with each metabolic disorder. Our results provide insights into the mechanistic role of gut microbiome in metabolic diseases and how this role is potentially similar between these three pathologies. Our models predicted the results of previous studies such as an increase in acetate (Schwiertz et al., 2010) and depletion of butyrate levels (Arora and Bä ckhed, 2016). In addition, our modeling indicated a disparity in the production of amino acids including glutamate, proline, tyrosine, and valine in the gut microbiome between the healthy and disease cohorts. Performing reaction abundance analysis between cohorts demonstrated multiple reactions involving tartrate metabolism which were enriched in all three disorders. We performed personalized community modeling analysis for two of the cohorts, leading to further clarification of bacterial metabolism differences between healthy and diseased subjects. These findings were further evaluated by investigating the link between the tartrate metabolism and circulating amino acid levels using the host blood metabolites thus showing the potential impact of the tartrate on the host metabolism.

RESULTS
To investigate the individual metabolic roles of the microbiome in metabolic disorders of the host, we used gut metagenomics data from six previous studies on obesity, T2D, and ACVD (Jie et al., 2017;Karlsson et al., 2012Karlsson et al., , 2013bLe Chatelier et al., 2013b;Qin et al., 2010Qin et al., , 2012a. This amounted to 1,443 subjects in total; 278 obese patients and 263 matched controls; 271 patients with T2D and 231 matched controls; and 219 patients with ACVD and 181 matched control samples (Tables 1 and 2). Previously, we have analyzed these data with an updated gut gene catalog and MSP profile (Shoaie et al., 2021) while here we performed the analysis to choose the significant MSPs and their corresponded GEMs. We then applied individual-level constraint-based modeling of these MSPs to understand the contribution of each species to the overall metabolic changes. In addition, we performed personalized community modeling together with personalized reaction abundance analysis. We also investigated the plasma metabolite association from a separate cohort of Swedish patients to evaluate the modeling predictions. iScience Article Identifying common and unique metagenome species across the three metabolic disorders Phyla level analysis of the MSP profiles between patient and control groups in all three disorders showed Bacteroidetes and Firmicutes had consistent significant differences (FDR <0.01) ( Figure S1). Comparing obese samples to their matched controls showed an increase in the abundance of Bacteroidetes, and decrease in the abundance of Melainabacteria, Verrucomicrobia, Tenericutes, Synergistetes, and Firmicutes phyla. In contrast, we found that the abundance of Bacteroidetes was decreased in patients with T2D and ACVD while other phyla, including Actinobacteria and Firmicutes, were increased.
For a deeper insight at the species level, we identified significant MSPs within patients from all diseases (FDR <0.01) (Tables S1, S2, and S3). The highest abundance of enriched MSPs (fold-change) was in ACVD (65%), while there were 44 and 49% disease-enriched MSPs in obesity and T2D, respectively ( Figure S2). We identified 202 statistically different (FDR<0.01) MSPs between patient and control groups in all cohorts; 139 significant MSPs in Obesity, 12 significant MSPs in T2D and 51 significant MSPs in ACVD.
To evaluate the median abundance of the MSPs for the patients and matched controls, MSPs were selected based on having different medians between the two groups (17 MSPs for obesity, 6 for T2D, and 25 for ACVD). Grouping these together and removing any duplicated MSPs yielded 42 unique significant MSPs ( Figure 1) ( Figure S3, Table S4), with 33 MSPs out of the selected 42 from the phylum Firmicutes. The chosen 42 MSPs were not all significant in each of the 3 disorders. To determine if each MSP was enriched in the disease cases, we only considered the disease in which the MSP showed significance. This gave 15 disease-enriched MSPs and 26 control-enriched MSPs. Out of the six MSPs significant in ACVD and T2D, four were depleted in both diseases and one was enriched while Clostridium phoceensis was not consistent across the two diseases, being decreased in patients with ACVD but increased in patients with T2D.
Genome-scale metabolic modeling showed changes in short-chain fatty acids and amino acid metabolism in metabolic diseases In order to investigate the microbe-microbe interactions, we applied 37 functional GEMs for the MSPs associated with pathology. 11 models represented MSPs increased in ACVD, T2D, and obesity while 26 models represented MSPs decreased in all three disorders. A further model represented the inconsistent species C. phoceensis identified above. The selected models contained on average 1,717 reactions, 1,684 metabolites, and 1,022 genes with an average growth rate of 0.89 h À1 ( Figure S4). The majority of the selected GEMs belonged to Lachnospiraceae and Ruminococcaceae families with 19 and 5 models, respectively. Furthermore, we showed that the models are metabolically distinct as determined by their Jaccard index, the average Jaccard index was equal to 0.637 ( Figure 2A, Table S5). The models were assessed by constraining exchange reactions with a high fiber omnivorous diet as an input and setting biomass as the objective function (Bidkhori et al., 2021). We clustered the flux predictions from the models based primarily on the models' enrichment in diseased patients and secondly on which disease the MSP was significantly enriched/depleted in ( Figure 2B). The metabolism of sugars and amino acids was predicted in these models. Using sucrose, starch, fructose, and glucose which are available from the high fiber diet as inputs, models showed they produce indole, ammonia, and amino acids. In the species-level simulations, SCFAs were also produced; 35 models secreted acetate, 10 secreted butyrate, and 12 secreted propionate (Table S6). After comparing the fluxes in individual GEMs, we compared the average flux across those MSPs increased in disorder versus decreased. Average flux based on microbiome compositional changes in diseases showed the flux exchange of key metabolites such as SCFAs and amino acids was altered. Acetate and propionate flux were increased in species with higher abundances in patients from iScience Article all disease cohorts ( Figure 2C). The predicted flux of butyrate showed lower amounts being produced by species increased in disease. There were also notable differences between the branch chain amino acid (BCAA) flux profiles, with average valine production higher than those MSPs increased in disease. Finally, there were increases in the average production of ammonia and proline, as well as increased consumption of glutamate, glycine, alanine, and serine by those models representing species enriched in disease.

Enrichment of tartrate metabolism among the three metabolic disorders
To further understand the functional behavior of the microbiome, we performed reaction abundance analysis for each individual using the MIGRENE Toolbox (Bidkhori et al., 2021) (Tables S7, S8, and S9, STAR Methods). A reaction abundance matrix was constructed for all samples, based on the MSP abundance table and the reactions present within the GEMs. Statistical analysis was then performed based on comparing the diseased cohorts to their matched controls. There were in total 820 significant reactions in obesity, 939 reactions in T2D, and 1,105 reactions in ACVD, compared to their healthy cohorts (FDR <0.01). From this, there were 214 unique reactions that were significant and present in all three diseases. We mapped the 214 significant reactions to KEGG orthologs (KOs). This resulted in 86 significant KOs (FDR <0.01) ( Figure S5, Table S10). However, only five of these reactions from two different KOs were enriched consistently across all three diseases. To define the molecular functionality of these reactions,  For each bacterial species, the red bar shows represent the enrichment in T2D, the blue bar represents the enrichment in obesity and the yellow bar represents the enrichment in ACVD. Fold-change>0 shows the species is increased in disease, fold-change <0 shows the species is decreased in disease. An asterisk (*) next to a bar indicates the bacterial species was statistically significant in that disease when compared to the matched controls. A filled-in black circle next to the species name indicates a genome-scale metabolic model was available for that bacterial species, an empty circle shows there is no model available. Microbial community modeling shows decreased production of acetate and butyrate in patients with T2D Following the analysis of individual GEMS and their potential contribution to metabolism in the gut microbiome, it is important to look at the microbial community as an integrated system. Hence, we developed 147 community models for each individual based on their species abundance for the two Swedish cohorts (93 T2D case, 36 T2D control, 5 ACVD case, and 10 ACVD control), using the MIGRENE Toolbox (Bidkhori et al., 2021). We ran simulations on these community models and implemented an optimization for the total biomass of the community model (Method). This was maximized such that the biomass of individual models within the community was maximized. The community models were constrained with a high fiber omnivorous diet, which only limited the uptake exchange reactions in the models. The average biomass production was 1.64 mmol/gDW/hr for the community models.
The community models for individuals in the Swedish T2D cohort showed a common theme of lower production of key metabolites in the T2D communities. We observed on average the T2D community models predicted lower secretion of acetate and butyrate ( Figure 4A), BCAAs ( Figure 4B) and lysine, tyrosine, and histidine ( Figure 4C) when compared to the average production flux from the control communities.
The community modeling for individuals in ACVD cohort predicted higher production of the butyrate and acetate, and the tyrosine and histidine in the diseased models ( Figures 4D-4F), while the valine and the lysine were predicted with higher production in the control community models.
These findings are relevant as they highlight how the community behaves differently from the individual models. For example, the production of acetate which was shown to be highly secreted by the individual models increased in disease now shows lower production from the T2D community models when compared to the control community models. Also, in the ACVD community models, we observed the BCAA production differs from that predicted in the individual modeling. To evaluate the importance of the microbiome tartrate and glyoxylate cycle in obesity, we investigated the plasma level of metabolites and gut metagenome from publicly available data on the longitudinal Swedish wellness study (Olsson et al., 2022;Tebani et al., 2020). Plasma metabolite profile allows us to link the gut microbial metabolism to the host physiology. This study used 101 healthy subjects with no clinically significant health complications ranging between 50 and 65 years of age. In this study stool and plasma samples have been taken from the patients. Previously, we performed the reaction abundance analysis on the gut microbiome of these individuals (Bidkhori et al., 2021).
We performed association analysis between the reaction abundance profile and metabolomics using multivariate association analysis. This metabolic analysis showed the association between reactions from the glyoxylate and dicarboxylate pathway and plasma metabolites in obese individuals of this cohort compared to non-obese ones. We used plasma metabolite profiles and BMI to find associations between iScience Article the abundance of glyoxylate and dicarboxylate pathways in the gut using a multivariate random-effects model (MVEM). Focusing on reactions relating to glyoxylate metabolism we found associations between 14 reactions and 45 plasma metabolites belonging to carnitines, fatty acids and lipids, amino acids, steroids, and organic acids. (FDR <0.01, Figure S6, Table S11).
Malate dehydrogenase [NAD(P)+] in particular had five positive associations with proline, ketoleucine, histidine, glutamate, and alanine. Other reactions which were positively associated with alanine plasma levels were tartrate decarboxylase and tartrate dehydratase along with four other enzymes ( Figure S7). Other amino acids which showed positive associations reactions were proline which had five positive associations and tyrosine which showed four. Whilst the amino acids histidine, tryptophan and ornithine all showed multiple negative associations with reactions within the glyoxylate metabolism.
To investigate if these microbial reactions were also associated with BMI, we split the Swedish cohort by individual patients' BMI and compared the two subgroups. We ran statistical tests (Wilcoxon rank-sum test) to ascertain if the reactions were correlated with one group or the other. Of the 14 reactions, 7 showed a significant positive association with BMI, including tartrate dehydratase, tartrate decarboxylase, and malate dehydrogenase. 4 reactions (malate dehydrogenase (oxaloacetate-decarboxylating), malate dehydrogenase (oxaloacetate-decarboxylating) (NADP+), malate dehydrogenase [NAD(P)+], tartrate decarboxylase) were significantly increased in abundance in samples from obese participants compared to other non-obese participants (P-value < 0.05 Wilcoxon rank-sum test, Figure 5, Table S12). The increase of these reactions in the obese Swedish subjects validates the discovery of the enrichment of the bacterial glyoxylate and dicarboxylate metabolism in obese cohorts.

DISCUSSION
This study was undertaken to gain deeper understanding of gut microbiome dysbiosis and the role of bacteria in metabolic disorders. Our increasing understanding of the interactions between the gut microbiome and its host during disease is poised to reveal the mechanistic role of the microbiome and its collective metabolism on human health and disease pathophysiology. Metabolic diseases have all previously been linked to changes in the gut microbiome, but to date, they have not been looked at in an integrative fashion. This study demonstrates differences in three metabolic diseases when compared to healthy controls and highlights the link between the diseases.
Metagenomic analysis identified the bacteria present during dysbiosis, including highly prominent species from the Melainabacteria, Verrucomicrobia, Tenericutes, Synergistetes, and Firmicutes phyla. Our results showed Clostridium bolteae, Clostridium symbiosum, Eggerthella lenta, and Escherichia coli were increased in T2D, and these species have been shown to be increased in T2D before (Karlsson et al., 2013a;Qin et al., 2012b). E. lenta and E. coli have also been shown to be increased in ACVD before (Jie et al., 2017). Ruminococcus gnavus, Ruminococcus torques, and Klebsiella pneumoniae were also noted in our results as increased in ACVD and obesity, which have also previously been noted to be increased in ACVD and obesity (Jie et al., 2017;Le Chatelier et al., 2013b). Although Eubacterium eligens, Clostridiales bacterium, Roseburia intestinalis, and Faecalibacterium prausnitzii were all depleted in our results and have been shown to also be decreased in metabolic diseases before (Karlsson et al., 2013a;Qin et al., 2012b).
Further analysis using GEMs predicted phenotypic behavior -in particular, the substrate uptake and secretion profiles for SCFAs and amino acids. A large proportion of the models were from Lachnospiraceae, a genus whose members are known to ferment polysaccharides to SCFAs (Boutard et al., 2014). We further explored the bacterial behavior within an ecosystem by constructing personalized community models. These community models showed notable differences in the production of SCFAs and AAs, in particular BCAAs.
Based on individual model simulations, we predicted an increase in acetate and propionate production from species that were increased in the microbiome of patients of metabolic diseases. Acetate has been closely linked to inflammasome activation before, seeming to have a reducing effect on the inflammation (Xu et al., 2019). Propionate has also shown anti-inflammatory effects, suppressing NF-kB reporter activity (Tedelind et al., 2007). Elevated levels of BCAAs have been observed in patients with obese and T2D (Sikalidis and Maykish, 2020) and they have become biomarkers of insulin resistance helping predict iScience Article the development of T2D (Zhu and Goodarzi, 2020). In addition, individual modeling showed an increase in glutamate consumption in species enriched in disease. Glutamate can induce obesity when administered in rodents, as well as has a positive correlation with obesity in Chinese adults (Lee et al., 2017). This individual modeling, therefore, shows that the metabolome profile of MSPs associated with these metabolic diseases is associated with the inflammation and physiology present in them. Here, we have shown how the gut microbiome, inflammation, and metabolic disease could be linked via these metabolites.
Reaction abundance analysis gave a deeper understanding of the mechanisms behind these results.
Consistently throughout the three explored diseases, we found increased abundance of reactions in glyoxylate and dicarboxylate metabolism, in particular reactions surrounding tartrate. Tartrate is usually found in foods such as grapes and when consumed, it feeds into metabolism via the TCA cycle or is converted to glycerate (Kohn and Jakoby, 1968). Increased plasma levels of glycerate are positively correlated with T2D (Reddivari et al., 2017). Only 20% of tartrate ingested in food is eliminated in urine meaning the remaining ingested tartrate is potentially consumed by the microbes in the gut (Finkle, 1933   iScience Article bulk of it (Chadwick et al., 1978) highlighting the close connection between tartrate and the gut microbiome. Enteric bacteria can use both malate and tartrate (C 4 -dicaroboxylates) for anaerobic growth . Furthermore, tartrate can feed into proline and arginine metabolism, potentially causing increased production of these metabolites. Increased plasma levels of amino acids including proline and arginine have been seen in disease and could potentially lead to exacerbating disease (Proffitt et al., 2020).
Glyoxylate and dicarboxylate metabolism has previously been linked to atherosclerosis (Chen et al., 2018) and to obesity (Beaumont et al., 2016). The glyoxylate cycle can metabolize fatty acids to glucose which can then contribute to insulin resistance (Song, 2000). As tartrate (a metabolite within the glyoxylate and dicarboxylate metabolic pathways) has been associated with glutamate, this metabolite could be the missing link in the causative effect of glutamate on increased BMI. The increased production of acetate by the GEMs seen during metabolic disease could be linked to the glyoxylate and dicarboxylate metabolism as these pathways catabolize acetate and amino acids for energy production in the microbe (Maniscalco et al., 2017). There was also an increase in the arginine and proline metabolism. This increase also has a connection to the increased consumption of glutamate, which is necessary for the synthesis of arginine and proline (Cappelletti et al., 2018).
The consistency in the present results highlighting tartrate metabolism provide strong support for our models correctly predicting the phenotypes observed in the three studied metabolic diseases. The finding that reactions involved with tartrate along with functional findings from the plasma metabolite data in the separate study from Swedish cohorts provides strong support for this hypothesis. These findings validated the predictions of increased production of acetate, proline, and arginine from the metabolic modeling. It is also known that malate dehydrogenase [NAD(P)+] and tartrate dehydratase both produce oxaloacetate while malate dehydrogenase (oxaloacetate decarboxylate) produces pyruvate . Both oxaloacetate and pyruvate are used in the amino acid synthesis where aspartate, alanine, asparagine, methionine, lysine, and threonine are synthesized (Marty-Teysset et al., 1996), further validating our findings.
From these findings, we would hypothesis that an increased abundance of tartrate metabolism is not beneficial for health. This is based on the reactions showing an increase across all three of our diseases, also the fact that the gut microbiome metabolizes this metabolite. Tartrate is fermented in the colon to SCFAs (EFSA Panel on Food Additives and Flavourings et al., 2020), we saw acetate had a higher production from the individual modeling which correlates with tartrate metabolism producing SCFAs. Increased acetate can affect host metabolism by causing increased pancreatic b-cell activity, glucose-stimulated insulin secretion (GSIS), hyperphagia, and obesity (Bose et al., 2019). To take this work furthermore, computational protein or gene expression analysis could ensure the presence of the tartrate pathway within the bacterial species we have seen to be enriched in diseases. Additionally, more metabolomic analysis of stool and plasma from patients with T2D or ACVD could further validate these findings. Taking this idea further microbial culture metabolomics could be analyzed for comparison with the patient metabolomics to confirm the hypothesis from patient metabolomics.
This analysis can be used to unravel the association between impaired gut microbiomes and metabolic disorders. We have identified tartrate metabolism in the microbiome as a significant pathway with high potential to impact on each of these metabolic syndromes. As well as this giving a potential new biomarker, it also suggests potential novel intervention targets for disease. Further research into the metabolism of tartrate in the gut environment is needed to understand the direct impact this compound has on gut microbes and host health. It is clear that the gut microbiome is a key factor in maintaining health and our highlighted novel areas for research give areas to focus on for both disease understanding and biomarker discovery for obesity and its co-morbidities.

Limitations of the study
There are a number of limitations to our study. Firstly, the cohorts chosen to represent the metabolic diseases are from different geographical regions which will have its own impact and implications on the gut microbiome profile. There are also limitations within the modeling of the bacteria. There are several bacterial species in the gut (Almeida et al., 2019), while we only had 1, 333 bacterial species-specific models available from the MIGRENE toolbox. Another limitation within the modeling was the computational power ll OPEN ACCESS iScience 25, 104513, July 15, 2022 iScience Article limits which restricted the number of models within the community modeling to 10 GEMs per community. Clearly, if this could have been increased the community models would have been a better representation of the ecosystem within the gut microbiome.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: iScience Article iScience Article steps: (1) the MSP's FDR score, FDR<0.01, (2) distinctly different median values in the patients and controls, one equal to 0 and one equal to a real positive number. We estimated the MSP enrichment or depletion in diseases by comparing the average abundance in case samples, and the average abundance in control samples per MSP. The output identified each MSP as either enriched in disease or enriched in control.

Longitudinal Swedish cohort samples metabolomics
The initial analysis for the associations between the plasma metabolite levels and reaction abundancies within the gut microbiome was done previously (Bidkhori et al., 2021), using available metabolomics and metagenomics (Olsson et al., 2022;Tebani et al., 2020). Here, the plasma and stool samples from the Swedish cohort were linked using the reaction abundance. The plasma samples gave circulating metabolite profile and the stool samples gave reaction abundance from the gut microbiome. Following this we ran multivariable association analysis on the data to find the corelation between the clinical metadata and microbial omics features. This analysis was done in R using packages taRifx, Maaslin2 and car.