# Parental Selenium Nutrition Affects the One-Carbon Metabolism and the Hepatic DNA Methylation Pattern of Rainbow Trout (Oncorhynchus mykiss) in the Progeny

Pauline Wischhusen, Takaya Saito, Cécile Heraud, Sadasivam Kaushik, Benoit Fauconneau, Philip Antony Jesu Prabhu, Stéphanie Fontagné-Dicharry, Kaja

Skjaerven

► To cite this version:

Pauline Wischhusen, Takaya Saito, Cécile Heraud, Sadasivam Kaushik, Benoit Fauconneau, et al. Parental Selenium Nutrition Affects the One-Carbon Metabolism and the Hepatic DNA Methylation Pattern of Rainbow Trout (Oncorhynchus mykiss) in the Progeny. Life, 2020, 10 (8), 10.3390/life10080121. [hal-02911045]

HAL Id: hal-02911045

https://hal.inrae.fr/hal-02911045

Submitted on 3 Aug 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

life

<figure><img src="image_3.png" /></figure>

Article

# Parental Selenium Nutrition Affects the One-Carbon Metabolism and the Hepatic DNA Methylation Pattern of Rainbow Trout (*Oncorhynchus mykiss*) in the Progeny

Pauline Wischhusen <fnref n="1" />,† ID, Takaya Saito <fnref n="2" />,† ID, Cécile Heraud <fnref n="1" /> ID, Sadasivam J. Kaushik <fnref n="1" /> ID,
Benoit Fauconneau <fnref n="1" /> ID, Philip Antony Jesu Prabhu <fnref n="2" /> ID, Stéphanie Fontagné-Dicharry <fnref n="1" />,*,† ID
and Kaja H. Skjærven <fnref n="2" />,† ID

Received: 11 June 2020; Accepted: 23 July 2020; Published: 25 July 2020

<figure><img src="image_4.png" /></figure>

check for
updates

**Abstract:** Selenium is an essential micronutrient and its metabolism is closely linked to the methionine cycle and transsulfuration pathway. The present study evaluated the effect of two different selenium supplements in the diet of rainbow trout (*Onchorhynchus mykiss*) broodstock on the one-carbon metabolism and the hepatic DNA methylation pattern in the progeny. Offspring of three parental groups of rainbow trout, fed either a control diet (NC, basal Se level: 0.3 mg/kg) or a diet supplemented with sodium selenite (SS, 0.8 mg Se/kg) or hydroxy-selenomethionine (SO, 0.7 mg Se/kg), were collected at swim-up fry stage. Our findings suggest that parental selenium nutrition impacted the methionine cycle with lower free methionine and S-adenosylmethionine (SAM) and higher methionine synthase (mtr) mRNA levels in both selenium-supplemented treatments. DNA methylation profiling by reduced representation bisulfite sequencing (RRBS) identified differentially methylated cytosines (DMCs) in offspring livers. These DMCs were related to 6535 differentially methylated genes in SS:NC, 6890 in SO:NC and 7428 in SO:SS, respectively. Genes with the highest methylation difference relate, among others, to the neuronal or signal transmitting and immune system which represent potential targets for future studies.

**Keywords:** selenium; methionine cycle; transsulfuration; nutritional programming; DNA methylation; rainbow trout

## 1. Introduction

Selenium (Se) is an essential micronutrient in humans and animals, with selenoproteins exerting various metabolic functions [1]. Among vertebrates, fishes have been described to have a well-developed selenoproteome [2], but there is concern within the aquaculture sector with present feed formulations. The ongoing replacement of Se-rich fishmeal with plant protein sources [3,4] is associated with a decrease in dietary Se level provided to farmed fish reared over a long period [5,6]. Many of the characterized selenoproteins are known to influence antioxidant metabolism, but knowledge on the effects of dietary Se on other metabolic pathways is not well characterized [7].

As shown *in vivo*, in the case of a Se deficiency, an increase in glutathione levels possibly relates to a feedback mechanism by changes in redox state [8,9]. The major source for glutathione is cysteine,

<footnote n="1">¹ INRAE, University of Pau and Pays de l'Adour, E2S UPPA, NUMEA, 64310 Saint Pee sur Nivelle, France; pauline.wischhusen@inrae.fr (P.W.); cecile.heraud@inrae.fr (C.H.); sadasivam.kaushik@inrae.fr (S.J.K.); benoit.fauconneau@inrae.fr (B.F.)</footnote>

<footnote n="2">² Institute of Marine Research, P.O. Box 1870, 5020 Bergen, Norway; takaya.saito@hi.no (T.S.); antony.philip@hi.no (P.A.J.P.); kaja.skjaerven@hi.no (K.H.S.)</footnote>

<footnote>* Correspondence: stephanie fontagne-dicharry@inrae.fr</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
which is synthesized from homocysteine via the transsulfuration pathway [10]. On the other hand, Se deficiency can impair the transsulfuration pathway with decreased levels of cysteine, cystathionine and homocysteine [11,12]. Homocysteine is a key metabolite in the methionine cycle, which links the antioxidant system to one-carbon (1C) metabolism [13]. Studies in mice have confirmed that Se affects methionine metabolism with decreased betaine homocysteine methyltransferase activity and S-adenosylhomocysteine (SAH) levels [14–16]. Comparable studies in fish are lacking.

In the methionine cycle, methionine is activated by S-adenosylmethionine synthetase to form S-adenosylmethionine (SAM). In the cell, SAM is the universal donor for methylation reactions forming SAH. DNA methylation is a major regulatory mechanism for epigenetic modifications [17]. DNA methylation at repeated cytosine phosphate guanine (CpG) residues, especially when localized at the promoter region, is considered to influence gene expression [18]. Dietary supplementation of Se has been associated with both hyper- and hypomethylation in mice, but the relationship between Se and epigenetic mechanisms is still not fully understood [19]. The present work therefore aims to study the effect of parental Se nutrition in rainbow trout (*Oncorhynchus mykiss*) on the 1C metabolism and the hepatic DNA methylation pattern of the progeny.

The period of embryonic development is extremely sensitive to environmental-induced epigenetic modifications. For example, the allocation of maternal gene products and nutrients to the yolk has been associated with regulation of key embryonic developmental processes and persisting changes in the phenotype of the progeny [20,21]. In zebrafish (*Danio rerio*), dietary inclusion of methyl group donors did not lead to changes in hatching rate or survival, but mRNA sequencing of the embryos revealed “hidden” effects of parental nutrition [22] which led to phenotypic changes at later life stages [23]. In rainbow trout, the maternal Se nutrition during oogenesis increased not only the number of spawning females, but also the Se levels in the oocytes, especially when provided in the form of organic Se [6]. Changing Se levels in the progeny during embryonic development were also associated with modifications in the oxidative status.

In fish diets, Se supplementation becomes increasingly important to make up for the low Se levels detected in diets based on plant protein sources [24]. With regard to Se supplements, in addition to the widespread use of sodium selenite in terrestrial livestock nutrition, selenomethionine is the naturally dominant dietary seleno-compound, known to be a highly bioavailable form of Se also in mineral premixes [25]. These seleno-compounds, however, might exert different impacts on the 1C metabolism as they are metabolized through different routes [26]. Seleno amino acids are metabolized interchangeably with their sulfur analogues making selenomethionine to follow the methionine cycle, while inorganic Se compounds such as sodium selenite can be directly reduced to selenide to be incorporated into selenoproteins as selenocysteine [24].

In this context, the present study aims to make a comparison between the use of sodium selenite and hydroxy-selenomethionine (OH-SeMet), a pure form of the hydroxy-analogue of selenomethionine, as dietary supplements in plant protein-rich feeds for rainbow trout broodstock on the 1C metabolism and the hepatic DNA methylation pattern of the progeny.

## 2. Results

### 2.1. Parental Selenium Affects Transsulfuration Metabolites in Swim-Up Fry

A decrease in cysteine and cysteinyl-glycine was detected in liver of female broodstock only when fed sodium selenite (SS) compared to the non-supplemented control (NC). In addition, homocysteine levels were higher in fish fed OH-SeMet (SO) compared to the two other groups. No effect of the dietary Se on hepatic aminothiol concentrations was detected in males, which had generally lower hepatic aminothiol levels compared to females (Table 1).
<table><caption>**Table 1.** Aminothiol concentrations (µg/g sample) measured in liver and oocytes of rainbow trout (*Oncorhynchus mykiss*) broodstock fed diets containing different levels and source of Se.</caption><thead><tr><th></th><th></th><th>Homocysteine</th><th>Cysteine</th><th>Cysteinyl-Glycine</th><th>Glutathione</th><th>γ-Glutamyl-Cysteine</th></tr></thead><tbody><tr><td rowspan="3">Oocyte</td><td>NC</td><td>0.3 ± 0.0</td><td>6.7 ± 1.0</td><td>3.7 ± 0.4</td><td>16 ± 1</td><td>1.0 ± 0.1</td></tr><tr><td>SS</td><td>0.3 ± 0.0</td><td>6.7 ± 0.6</td><td>4.1 ± 0.5</td><td>17 ± 1</td><td>1.0 ± 0.1</td></tr><tr><td>SO</td><td>0.2 ± 0.0</td><td>8.1 ± 1.1</td><td>3.1 ± 0.3</td><td>13 ± 1</td><td>1.0 ± 0.1</td></tr><tr><td></td><td><i>p-value</i></td><td>0.21</td><td>0.48</td><td>0.35</td><td>0.07</td><td>0.48</td></tr><tr><td rowspan="3">Female liver</td><td>NC</td><td>1.4 ± 0.1 <sup>b</sup></td><td>33 ± 4 <sup>a</sup></td><td>53 ± 4 <sup>a</sup></td><td>551 ± 32</td><td>30 ± 4</td></tr><tr><td>SS</td><td>1.1 ± 0.1 <sup>b</sup></td><td>17 ± 2 <sup>b</sup></td><td>37 ± 3 <sup>b</sup></td><td>530 ± 32</td><td>23 ± 2</td></tr><tr><td>SO</td><td>3.2 ± 0.3 <sup>a</sup></td><td>38 ± 3 <sup>a</sup></td><td>48 ± 3 <sup>ab</sup></td><td>526 ± 47</td><td>28 ± 3</td></tr><tr><td></td><td><i>p-value</i></td><td><strong>&lt;0.01</strong></td><td><strong>&lt;0.01</strong></td><td><strong>0.01</strong></td><td>0.88</td><td>0.29</td></tr><tr><td rowspan="3">Male liver</td><td>NC</td><td>0.7 ± 0.1</td><td>21 ± 3</td><td>26 ± 3</td><td>511 ± 64</td><td>21 ± 2</td></tr><tr><td>SS</td><td>1.0 ± 0.2</td><td>24 ± 5</td><td>23 ± 3</td><td>300 ± 66</td><td>10 ± 2</td></tr><tr><td>SO</td><td>1.1 ± 0.3</td><td>33 ± 8</td><td>27 ± 4</td><td>465 ± 30</td><td>16 ± 5</td></tr><tr><td></td><td><i>p-value</i></td><td>0.49</td><td>0.33</td><td>0.80</td><td>0.06</td><td>0.48</td></tr><tr><td rowspan="2">Average</td><td>Female</td><td>1.9 ± 0.2 <sup>a</sup></td><td>29 ± 2</td><td>46 ± 2 <sup>a</sup></td><td>536 ± 22 <sup>a</sup></td><td>27 ± 2 <sup>a</sup></td></tr><tr><td>Male</td><td>0.9 ± 0.1 <sup>b</sup></td><td>26 ± 3</td><td>25 ± 2 <sup>b</sup></td><td>434 ± 39 <sup>b</sup></td><td>13 ± 2 <sup>b</sup></td></tr><tr><td></td><td><i>p-value</i></td><td><strong>&lt;0.01</strong></td><td>0.43</td><td><strong>&lt;0.01</strong></td><td><strong>0.02</strong></td><td><strong>&lt;0.01</strong></td></tr></tbody></table>

Values are the mean ± SEM (n = 8 in female tissue and n = 5 in males). <sup>a,b</sup> Within-rows values not sharing a common superscript letter are significantly different (p < 0.05) according to one-way ANOVA followed by Tukey's HSD.

A parental effect of Se in swim-up fry could be detected for cysteine as well as cysteinyl-glycine, which were both significantly lower in fry originating from parents fed Se-supplemented diets compared to the NC group (Table 2). This was accompanied by a decrease in pyridoxamine levels, but other B vitamins (folate and vitamin B12) including the pyridoxamine derivate pyridoxal were not significantly affected. Cystathionine, glutathione and γ-glutamyl-cysteine levels were not significantly different between the Se treatments. Similarly, parental Se treatment had no significant effect on the homocysteine level detected in swim-up fry.

<table><caption>**Table 2.** Free amino acid, aminothiol, and B vitamin composition of swim-up fry from broodstock fed the different diets.</caption><thead><tr><th>Dietary Group</th><th>NC</th><th>SS</th><th>SO</th><th><i>p</i>-Value</th></tr></thead><tbody><tr><td>Essential amino acids <sup>1</sup></td><td>1972 ± 79 <sup>a</sup></td><td>1737 ± 54 <sup>b</sup></td><td>1400 ± 65 <sup>c</sup></td><td><strong>&lt;0.01</strong></td></tr><tr><td>Non-essential amino acids <sup>1</sup></td><td>2415 ± 50 <sup>a</sup></td><td>2351 ± 41 <sup>a</sup></td><td>2109 ± 60 <sup>b</sup></td><td><strong>&lt;0.01</strong></td></tr><tr><td>Methionine <sup>1</sup></td><td>99 ± 5 <sup>a</sup></td><td>83 ± 3 <sup>b</sup></td><td>56 ± 4 <sup>c</sup></td><td><strong>&lt;0.01</strong></td></tr><tr><td>Homocysteine <sup>1</sup></td><td>1.2 ± 0.1</td><td>1.1 ± 0.1</td><td>1.2 ± 0.1</td><td>0.86</td></tr><tr><td>Cystathionine <sup>1</sup></td><td>9 ± 1</td><td>6 ± 1</td><td>7 ± 1</td><td>0.21</td></tr><tr><td>Cysteine <sup>1</sup></td><td>21 ± 1 <sup>a</sup></td><td>17 ± 1 <sup>b</sup></td><td>17 ± 0 <sup>b</sup></td><td><strong>0.01</strong></td></tr><tr><td>Cysteinyl-glycine <sup>1</sup></td><td>28 ± 1 <sup>a</sup></td><td>24 ± 1 <sup>b</sup></td><td>24 ± 1 <sup>b</sup></td><td><strong>0.02</strong></td></tr><tr><td>Glutathione <sup>1</sup></td><td>179 ± 7</td><td>159 ± 7</td><td>169 ± 13</td><td>0.35</td></tr><tr><td>γ-Glutamyl-cysteine <sup>1</sup></td><td>18 ± 1</td><td>17 ± 1</td><td>16 ± 1</td><td>0.25</td></tr><tr><td>Taurine <sup>1</sup></td><td>688 ± 17</td><td>751 ± 18</td><td>724 ± 16</td><td>0.05</td></tr><tr><td>Pyridoxamine <sup>2</sup></td><td>0.24 ± 0.01 <sup>a</sup></td><td>0.21 ± 0.02 <sup>b</sup></td><td>0.18 ± 0.01 <sup>b</sup></td><td><strong>0.01</strong></td></tr><tr><td>Pyridoxal <sup>2</sup></td><td>1.82 ± 0.08</td><td>1.65 ± 0.06</td><td>1.85 ± 0.10</td><td>0.17</td></tr><tr><td>Folate <sup>2</sup></td><td>0.36 ± 0.03</td><td>0.37 ± 0.02</td><td>0.28 ± 0.02</td><td>0.05</td></tr><tr><td>Cobalamine <sup>2</sup></td><td>0.04 ± 0.00</td><td>0.04 ± 0.00</td><td>0.03 ± 0.00</td><td>0.51</td></tr></tbody></table>

## 2.2. Parental Selenium Nutrition Affects the Methionine Metabolism in Swim-Up Fry

In the whole body of swim-up fry, the methionine concentration was significantly decreased, when parents received Se-supplemented diets compared to fry from the NC group, with the lowest concentration observed in fries from SO treatment (Table 1). The decreased methionine levels were accompanied by a general decrease in both essential and non-essential amino acids. PCA analysis of free amino acids and N-metabolites in fry revealed a strong clustering of the data according to the three parental groups dominated by essential amino acids, with the main contributing variables being lysine, isoleucine, valine, leucine, methionine, threonine and histidine besides glutamine, ammonium chloride

<footnote n="1">¹ (µg/g sample); ² (µg/mg sample). Values are the mean ± SEM (n = 8). <sup>a,b,c</sup> Within-rows values not sharing a common superscript letter are significantly different (p < 0.05) according to one-way ANOVA followed by Tukey's HSD.</footnote>
and glycine (Figure 1). The only amino acid that was significantly higher in Se-supplemented treatments compared to the control was asparagine, with a 18 ± 4 µg/mg sample in NC vs. a 34 ± 3 µg/mg sample in SS and a 44 ± 5 µg/mg sample in SO.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** PCA biplot of free amino acids and related compounds measured in whole-body swim-up fry. Arrows represent the 10 most contributing variables to the model. Ellipses represent the 95% confidence intervals around a center of eight pooled samples per dietary treatment.</figcaption></figure>

In broodstock liver tissue, a reduction in the SAM/SAH ratio was detected for the Se-supplemented groups (Figure 2A). SAM levels in males and females were, however, strongly affected by inorganic Se, without a significant difference between NC and SO in males that showed higher SAM as well as SAH levels compared to females. In oocytes, no significant difference in SAM or SAH levels and the SAM/SAH ratio was detected between treatments.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** (A) SAM, SAH and the SAM/SAH ratio in whole-body swim-up fry; (B) SAM, SAH and the SAM/SAH ratio in broodstock tissue. Bars are the mean ± SEM (n = 8 in swim-up fry and female tissues and n = 5 in males). Means not sharing a common superscript letter are significantly different (*p* < 0.05) according to one-way ANOVA followed by Tukey's HSD.</figcaption></figure>

In the whole body of swim-up fry, the SAM/SAH ratio was low in both Se-supplemented treatments compared to the control with the lowest SAM/SAH ratio detected in the SO group (Figure 2B). The decrease in the SAM/SAH ratio can be related to the comparatively lower SAM levels observed in this group.

## 2.3. Parental Selenium Affects mRNA Levels of Genes Related to the One-Carbon Metabolism in Swim-Up Fry

Gene expression levels in female liver tissue were not significantly different between groups. Methionine synthase (mtr) expression was higher in male liver tissue when the fish were fed Se-supplemented diets, with the highest expression in the SS treatment (Figure 3A). In addition,
in male liver tissue, the expression of adenosylmethionine decarboxylase 1a (*amd1a*) was higher in SO compared to NC and that of glycine N-methyltransferase (*gnmt*) in SS compared to the other two dietary treatments. Except for adenosylmethionine decarboxylase 1b (*amd1b*) and *gnmt*, the expression of the 1C metabolism-related genes analyzed was higher in the females than in the males.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** (A) Relative mRNA levels in whole-body swim-up fry from rainbow trout subjected to different Se treatments; (B) relative mRNA levels in parental liver tissue from rainbow trout subjected to different Se treatments. Data are normalized to β-actin and expressed as fold changes compared with the control group NC. In A, values are expressed relative to NC males. Bars are the mean ± SEM (A: n = 8; B n = 8 in female liver and n = 5 in male liver). Means not sharing a common superscript letter are significantly different (*p* < 0.05) according to one-way ANOVA followed by Tukey's HSD.</figcaption></figure>

Parental feeding of both SS and SO increased *mtr* gene expression in the swim-up fry compared to NC feeding (Figure 3B). In addition, the expression of *amd1b* was higher in SS compared to NC and that of adenosylhomocysteinase (*sahh*) was higher in SO compared to the two other groups.

## 2.4. Parental Selenium Resulted in a Weak Group-Wise DNA Methylation Clustering

Reduced representation bisulfite sequencing (RRBS) data were first processed and aligned to the rainbow trout genome (Table A1). For downstream analysis, only uniquely mapped reads (47.6 ± 1.4%) were used. Of the 12 samples sequenced (4 per dietary group), none appeared as an outlier. In the search for the global methylation pattern with all the mapped CpG sites, t-SNE (t-distribution stochastic neighbor embedding) was used. The individual methylation pattern was stronger compared to group-wise global patterns, with only a weak group-wise clustering identified when using the 95th
percentile of the CpC variance (Figure 4). The stronger individual variation compared to group-wise clustering was confirmed using other methods including PCA, hierarchical clustering and correlation analysis (Figure A1).

<figure><img src="image_8.png" /><figcaption>**Figure 4.** t-SNE analysis with the CpG sites using either the > 50% quantile or > 95% quantile.</figcaption></figure>

2.5. *Data Alignment Gives a Balanced Hepatic Methylation Pattern between Groups*

The regional annotation showed that most of the mapped CpG originated with 51.8% from the gene bodies compared to the whole rainbow trout genome, where it accounts for 22.9% (Figure 5A). With 37.6%, most of the mapped CpG in the gene body were coming from the intron region. Further, promoters were more targeted by RRBS compared to the whole rainbow trout genome, with an increase from 4.7% to 7.2%.

The total number of differentially methylated cytosine (DMC) was comparable for the groups SS:NC (10904), SO:NC (11806) and SO:SS (13179) and, within each group, the number of hyper- and hypomethylated CpG sites was balanced even when divided into different sub-regions, exon and intron for the body and P250 for the proximal promoter, P1K for the promoter and P6K for the distal promoter region as well as flanks (Figure 5B).
Regional distributions of CpGs

<figure><img src="image_9.png" /></figure>

(A)

Methylation differences by region (%)

<figure><img src="image_10.png" /></figure>

SS:NC (total: 5124 5780)

<figure><img src="image_11.png" /></figure>

SO:NC (total: 5518 6288)

<figure><img src="image_12.png" /></figure>

(B)

**Figure 5.** (A) Regional distributions of mapped/original CpG. (B) Regional distributions of methylation differences of the differentially methylated cytosines (DMC) with a 20% threshold. Three violin plots show the density of overall methylation differences for SS:NC, SO:NC, and SO:SS, with scattered dots indicating DMC.

2.6. Parental Selenium Affects the DNA Methylation Pattern in Several Metabolic Pathways

Comparing NC with the SS treatment showed a total of 6535 differentially methylated genes (DMGs) from which 1142 DMGs had DMCs located in the promoter region (Figure 6). Similarly, in total, 6890 DMGs were detected between NC and the treatment receiving SO with 1250 DMGs showing DMCs located in the promoter region. The highest number of DMGs was detected between the two different Se-supplemented treatments with a total of 7428 genes of which 1340 DMGs had DMCs located in the promoter region. A synergetic effect of Se in NC vs. SS and NC vs. SO was detected on 3663 genes, whereas SS vs. SO displays a specific effect of the Se source on 2387 genes.
<figure><img src="image_13.png" /><figcaption>**Figure 6.** Venn diagrams summarizing the analysis of genes with different methylation patterns in NC vs. SS, NC vs. SO and SS vs. SO. (A) Genes are included with at least one DMC in gene body or promoter region; (B) Genes are included with at least one DMC in the promoter region.</figcaption></figure>

In all datasets, multiple KEGG pathways were significantly enriched—22 in SS:NC, 18 in SO:NC and 20 in SO:SS (Figure A2). These KEGG pathways relate to diverse biological mechanisms, mainly cellular metabolism and environmental information processing, but also the organismal system and cellular processing.

Among the five genes with the highest number of DMCs listed for each of the sub-regions (exon, intron, proximal promoter, promoter and distal promoter) in Tables A2–A4, three genes were common. The limbic system-associated membrane, transcript variant X5 protein (Isamp) belongs to the immunoglobulin super-family, known to be expressed and excreted in the developing forebrain showed 15 DMCs in SS:NC, 24 in SO:NC and 22 in SO:SS, all located in the intron region. The methylation pattern revealed both hyper- and hypomethylated CpG sites. The DMCs of the other two genes were located in the promoter region. The radical S-adenosyl methionine domain containing protein 2-like (*viperin*) is a cytoplasmic antiviral protein that is induced by interferons. *Viperin* had five DMCs in SS:NC and three DMCs in SO:NC and SO:SS, respectively. In the inorganic Se treatment, the CpG sites were hypomethylated, but they were hypermethylated in the organic Se treatment. The third gene was gamma-aminobutyric acid receptor subunit rho-2 (*gabrr2*), which is an inhibitory neurotransmitter in the vertebrate brain. The DMCs of *gabrr2* were located in the distal promoter region and similar to *viperin* the gene was hypomethylated in SS:NC and hypermethylated in SO:NC.

## 2.7. Parental Selenium also Affects Methylation in Genes Related to the 1C Metabolism

Several genes related to the methionine cycle and transsulfuration pathway were identified, but mostly they contained only single DMC sites (Table 3). An effect on the genes that provide selenocysteine for the selenoprotein synthesis was only detected in the organic Se treatment. Nevertheless, selenoprotein I and selenoprotein U had DMCs in both Se-supplemented treatments.
<table><caption>**Table 3.** Differentially methylated genes (DMGs) related to sulfur and selenium metabolism.</caption><thead><tr><th colspan="2">DMGs</th><th colspan="3">Hyper-/Hypomethylated DMC</th></tr><tr><th>Gene ID</th><th>Gene name</th><th>SS:NC</th><th>SO:NC</th><th>SO:SS</th></tr></thead><tbody><tr><td colspan="5"><strong>Methionine Cycle</strong></td></tr><tr><td>110530927</td><td>S-adenosylmethionine synthase</td><td>0/1</td><td></td><td></td></tr><tr><td>110537066</td><td>S-adenosylmethionine synthase-like</td><td></td><td>0/1</td><td>0/2</td></tr><tr><td>110502651</td><td>S-adenosylmethionine synthase-like</td><td>1/0 <sup>P</sup></td><td>0/2</td><td>1/0</td></tr><tr><td>110529528</td><td>S-adenosylmethionine decarboxylase proenzyme-like</td><td></td><td>0/1</td><td>0/1</td></tr><tr><td>110538418</td><td>DNA (cytosine-5)-methyltransferase 1-like, transcript variant X1</td><td></td><td>0/1</td><td></td></tr><tr><td>110505844</td><td>DNA (cytosine-5)-methyltransferase 3A-like</td><td></td><td>1/0</td><td></td></tr><tr><td>110532515</td><td>DNA (cytosine-5)-methyltransferase 3A-like, transcript variant X2</td><td></td><td></td><td>1/0</td></tr><tr><td>110497603</td><td>DNA (cytosine-5)-methyltransferase 3A-like, transcript variant X6</td><td>1/0</td><td></td><td></td></tr><tr><td>110492301</td><td>DNA (cytosine-5)-methyltransferase 3B-like, transcript variant X1</td><td></td><td></td><td>1/0</td></tr><tr><td>110500231</td><td>Putative adenosylhomocysteinase 3</td><td></td><td>0/1 <sup>P</sup></td><td></td></tr><tr><td>110494352</td><td>S-adenosylhomocysteine hydrolase-like protein 1 transcript variant X1</td><td></td><td>1/0</td><td></td></tr><tr><td>110490243</td><td>Adenosylhomocysteinase 3-like</td><td>2/1</td><td>4/2</td><td>0/1</td></tr><tr><td>110522167</td><td>Putative adenosylhomocysteinase 3, transcript variant X2</td><td></td><td></td><td>0/1</td></tr><tr><td colspan="5"><strong>Glutathione Metabolism</strong></td></tr><tr><td>110521555</td><td>Glutamate—cysteine ligase regulatory subunit-like</td><td>0/1</td><td></td><td></td></tr><tr><td>110502703</td><td>Glutamate—cysteine ligase catalytic subunit-like, transcript variant X2</td><td>0/1</td><td></td><td></td></tr><tr><td>110532297</td><td>Glutathione synthetase</td><td></td><td>0/1 <sup>P</sup></td><td></td></tr><tr><td>110522620</td><td>Glutathione-specific gamma-glutamylcyclotransferase 1-like</td><td></td><td>2/0</td><td></td></tr><tr><td>110537206</td><td>Gamma-glutamyltransferase 5-like, transcript variant X1</td><td>0/1 <sup>P</sup></td><td>1/0 <sup>P</sup></td><td>2/0 <sup>P</sup></td></tr><tr><td>100305229</td><td>Glutathione S-transferase kappa 1, transcript variant X1</td><td></td><td></td><td>1/0 <sup>P</sup></td></tr><tr><td>110492369</td><td>Glucose-6-phosphate 1-dehydrogenase-like, transcript variant X2</td><td></td><td>2/0</td><td>0/1</td></tr><tr><td>100305228</td><td>Peroxiredoxin 6, transcript variant X2</td><td></td><td>1/0</td><td></td></tr><tr><td>110532317</td><td>Spermidine synthase</td><td></td><td>0/1 <sup>P</sup></td><td></td></tr><tr><td>110535309</td><td>5-oxoprolinase (ATP-hydrolyzing) Isocitrate dehydrogenase [NADP] cytoplasmic-like</td><td>1/0 <sup>P</sup></td><td>0/1</td><td>1/2</td></tr><tr><td>110501851</td><td>Isocitrate dehydrogenase [NADP] cytoplasmic-like</td><td></td><td>1/1</td><td>0/1 <sup>P</sup></td></tr><tr><td>110520228</td><td>Isocitrate dehydrogenase [NADP] cytoplasmic-like</td><td></td><td>1/1</td><td></td></tr><tr><td colspan="5"><strong>Selenoprotein Synthesis and Selenoproteins</strong></td></tr><tr><td>110494778</td><td>Methionyl-tRNA synthetase 1</td><td></td><td>0/1</td><td>0/1</td></tr><tr><td>110488988</td><td>Methionyl-tRNA synthetase 2, mitochondrial</td><td>1/0</td><td></td><td></td></tr><tr><td>110500323</td><td>Sep (O-phosphoserine) tRNA:Sec (selenocysteine) tRNA synthase</td><td></td><td></td><td>1/0 <sup>P</sup></td></tr><tr><td>110512109</td><td>tRNA selenocysteine 1-associated protein 1-like</td><td></td><td>0/1</td><td></td></tr><tr><td>110528137</td><td>Eukaryotic elongation factor, selenocysteine-tRNA specific</td><td></td><td>1/0</td><td></td></tr><tr><td>110529243</td><td>Selenoprotein I</td><td>0/3</td><td>5/0</td><td>7/0</td></tr><tr><td>100499413</td><td>Selenoprotein U</td><td>0/1 <sup>P</sup></td><td>1/0 <sup>P</sup></td><td>1/0 <sup>P</sup></td></tr><tr><td>110532070</td><td>Selenoprotein K, transcript variant X1</td><td></td><td>0/1</td><td></td></tr><tr><td>110497881</td><td>Selenoprotein O-like</td><td></td><td></td><td>1/0</td></tr><tr><td>110525853</td><td>Thioredoxin reductase 2</td><td></td><td></td><td>0/1</td></tr></tbody></table>

# 3. Discussion

## 3.1. Parental Selenium Nutrition and the Transsulfuration Pathway in the Progeny

Decreased levels of cysteine in rainbow trout fry originating from the parental group fed Se-supplemented diets are in contrast to reports in adult rats, mice and chicken, where rather an impaired transsulfuration with decreased cystathionine and cysteine levels has been described under

<footnote>P located at the promoter region.</footnote>
dietary Se deficiency [12,14,15]. Decreased levels of cysteine were also detected in liver tissue of female broodstock of the sodium selenite treatment, indicating an effect independent of life stage in rainbow trout (Figure 7). In the present study, swim-up fry of the Se-supplemented treatments had lower pyridoxamine levels compared to the control group, which indicates an increased demand for vitamin B6 by parental Se nutrition. The combined effect of maternal Se and pyridoxal nutrition has been studied in porcine embryos by Dalto et al. [27,28], who described that the co-supplementation increased plasma seleno-dependent glutathione peroxidase levels in the progeny in the long term. In addition, the supply of organic Se and vitamin B6 stimulated the expression of elongation factors, biological processes related to translation and the mitotic cell cycle in five-day-old embryos [27]. Although the sulfur and seleno amino acids follow a similar pathway, the most important reaction for selenocysteine is its reduction to selenide via selenocysteine lyase and further by selenophosphate synthetase that donates Se to the Sec-tRNA for the selenoprotein synthesis [26]. Both these enzymes are vitamin B6 dependent, highlighting the importance of this vitamin in Se metabolism. This reaction is independent of dietary Se form as both inorganic and organic forms undergo the reduction to selenide [29]. The impact of dietary Se form as observed in the present study could be questioned as in an earlier study it was shown that even in the parental sodium selenite treatment more than 94% of the Se in oocytes was either selenocysteine or selenomethionine [6]. Nevertheless, the higher selenomethionine levels corresponding to the higher total Se levels in the organic Se treatment might contribute to changes in the methionine cycle, providing methyl groups and homocysteine/Se-homocysteine for the transsulfuration pathway. Thus, it can be inferred that the higher redox status of Se-homocysteine compared to homocysteine might favor transsulfuration, as its enzymes are readily regulated through the redox status [30].

<figure><img src="image_14.png" /><figcaption>**Figure 7.** Effect of parental Se nutrition on the methionine cycle and transsulfuration pathway in the progeny of rainbow trout. Superscript indicates that the effect was only detected in the respective treatment.</figcaption></figure>

## 3.2. Parental Selenium Nutrition and the Methionine Cycle in the Progeny

In the present study, parental Se had no effect on homocysteine levels in the swim-up fry, while studies in mice indicate an inverse correlation between liver Se and homocysteine levels [31,32]. An increased mRNA level of *mtr* in both Se-supplemented treatments could indicate that in rainbow trout, parental Se favors the re-methylation of homocysteine to methionine in the offspring. An increase in homocysteine can result in the accumulation of SAH, which is a competitive inhibitor of methyl-transferases and therefore associated to global hypomethylation [19]. In mammals, selenomethionine supplementation resulted in decreased hepatic SAH, but selenite on the contrary increased hepatic SAH [16,33]. In the present study, the difference was not significant. Nevertheless,
increased mRNA levels of sahh indicate that SAH might be increasingly metabolized to homocysteine in the organic Se treatment (Figure 7). The decrease in the SAM/SAH ratio in swim-up fry of Se-supplemented treatments due to a decrease in SAM levels is possibly due to the lower methionine levels. In the OH-SeMet treatment, it cannot be excluded that a competition of selenomethionine on active transporters reduced the methionine uptake in the gut [34]. However, considering the small fraction that selenomethionine represents compared to dietary methionine levels in this study, it might be rather indicative of a higher methionine flux in the organic Se treatment. The sequence of Se-compounds in the metabolism creates an additional drain on methyl groups, as selenide and other highly reactive seleno-compounds can be spontaneously methylated [35]. This process is of importance for inorganic Se sources which do not follow the methionine cycle and are directly reduced to selenide [36]. This hypothesis is supported by our data on broodstock liver, where a significant decrease in SAM was only observed in the selenite treatment. This might be an explanation why Speckmann et al. [16] reported a high SAM/SAH ratio in response to selenomethionine supplementation using Se-deficient conditions in mice where no methylation of seleno-compounds for removal could be expected. Several other studies in mice using higher Se levels detected no effect of Se on the SAM/SAH ratio [14,33,37]. A depletion of the methyl donor SAM can result in decreased DNA methyltransferases (DNMT) activity [38]. If in human colon carcinoma cells, administration of selenite inhibits DNMT activity [39], *dnmt1* expression in the present study was not affected by Se, similar to what was reported on hepatic mRNA levels in mice [16].

### 3.3. Parental Selenium Nutrition Affects the DNA Methylation Pattern of Genes Related to Several Metabolic Pathways

Analyzed liver tissue of rainbow trout fry revealed that the DNA methylation patterns of several genes are sensitive to parental Se nutrition. Alterations in DNA methylation by Se have been reported in several murine studies, although with somewhat contradictory results, as Se could be associated with both hyper- and hypomethylation [16,33,37,40]. The methylation of DNA can possibly regulate the spatial-temporal expression pattern of genes driving towards the development of a specific phenotype [18]. Genes directly related to the sulfur and Se metabolism presented methylation differences according to parental Se nutrition. Therefore, epigenetic marks might relate to metabolic differences observed in the present as well as in an earlier study on the expression of genes involved in the glutathione and antioxidant metabolism in rainbow trout fry [6]. Although an expected enrichment of the glutathione pathway could not be detected, genes of the glutathione metabolism including glutathione synthetase and glutathione-s-transferase kappa 1 were detected as DMGs. It has been reported that in cancer cells, selenite supplementation reactivates the transcription of glutathione-s-transferase π, another member of the glutathione-s-transferase family by a hypermethylation of the promoter region [41]. Most studies with cancer cells generally report the methylation of selenoproteins like glutathione peroxidase 1 and 3, methionine sulfoxide reductase B1 and selenium binding protein 1 in the promoter region [42]. In the present study with rainbow trout, parental Se nutrition did not result in changing methylation pattern for these genes, contrary to selenoprotein I, a potential target of parental Se nutrition on the progeny. Selenoprotein I is a protein involved in the formation of the glycerophospholipid phosphatidylethanolamine [43], belonging to the glycophospholipid metabolism KEGG pathway that was enriched in the sodium selenite treatment. A silencing of the gene has been associated with impaired neural development as it is essential in the myelination process [44]. In general, several genes with high changes in DNA methylation like *gabrr2* were related to brain signaling pathways and neurotransmission. Under physiologically relevant conditions, Se nutrition has been associated with a neuroprotective role on γ-aminobutyric acidergic neurons [45]. It remains unclear whether the changes in DNA methylation of neuronal signaling genes as observed here in the hepatic tissue would be similarly detected in other organs such as the brain. DNA methylation works on a time and spatial dimension as genes gain tissue-dependent importance and are also activated and deactivated at different developmental stages [46,47]. High methylation
differences were also identified for genes with a role in immune protection, including the antiviral protein *viperin* that showed several DMC sites in the promoter region. *Viperin* was up-regulated by supra-nutritional Se feeding in rainbow trout as well as Atlantic salmon (*Salmo salar*) in earlier studies [48,49]. This indicates that the impact of Se on the inflammatory response in fish might not be limited to direct feeding effects, but also be exerted through an epigenetic process. In this context, similar to other natural feed additives [50], Se might act as an immunostimulant, improving fish immunity in the long term.

# 4. Materials and Methods

## 4.1. *Experimental Set up*

The experiment was conducted at the INRAE experimental fish farm in Lées-Athas, France. Fish maintenance and experimental procedures were conducted by trained personnel in compliance with the European Directive 2010/63/EU for the protection of animals used for scientific purposes and the French Decree no. 2013-118 for animal experimentation.

Three-year-old rainbow trout (*Oncorhynchus mykiss*) broodstock (initial mean weight: 1.1 ± 0.2 kg in females and 0.9 ± 0.3 kg in males) from the same genetic group produced at the INRAE facilities of Lées-Athas (permit no. A64-104-1) were individually tagged and divided into three groups consisting of 25 females and 15 males. The fish were reared under natural photoperiod, as previously described [6], over six months and fed the respective diets once daily to apparent satiation. At spawning, oocytes from eight females per group were fertilized with pooled sperm received from males of the same dietary treatment collected on the same day. Fertilized eggs from each female were reared separately in small trays until swim-up fry stage supplied with flow-through spring water at 8 ± 1 °C.

## 4.2. *Experimental Diets*

The diets were based on plant ingredients with an 8% fish oil inclusion and designed to differ only in their Se content (Table 4), as previously described [6]. The NC diet at a basal Se level of 0.3 mg/kg was not supplemented with Se. The SS diet was supplemented with sodium selenite to a target level of 0.6 mg/kg (analyzed concentration, 0.8 mg/kg) and the SO diet was supplemented to the same target level of 0.6 mg/kg with OH-SeMet (Selisseo®, Adisseo SAS, Antony, France), resulting in a final Se concentration of 0.7 mg/kg.
<table><caption>**Table 4.** Dietary composition.</caption><thead><tr><th>Diet</th><th>NC</th><th>SS</th><th>SO</th></tr></thead><tbody><tr><td><em>Ingredients</em></td><td></td><td></td><td></td></tr><tr><td>Plant meals <sup>1</sup></td><td>74</td><td>74</td><td>74</td></tr><tr><td>Crystalline amino acids and attractant mixture <sup>2</sup></td><td>3.14</td><td>3.14</td><td>3.14</td></tr><tr><td>Soybean lecithin <sup>3</sup></td><td>2</td><td>2</td><td>2</td></tr><tr><td>Fish oil <sup>3</sup></td><td>8</td><td>8</td><td>8</td></tr><tr><td>Vegetable oils <sup>4</sup></td><td>8</td><td>8</td><td>8</td></tr><tr><td>Astaxanthin (µg/g diet) <sup>5</sup></td><td>40</td><td>40</td><td>40</td></tr><tr><td>Vitamin and mineral mixture without Se <sup>6</sup></td><td>4.82</td><td>4.82</td><td>4.82</td></tr><tr><td>Sodium selenite (µg/g diet) <sup>7</sup></td><td>–</td><td><strong>0.71</strong></td><td>–</td></tr><tr><td>Hydroxy-selenomethione (µg/g diet) <sup>7</sup></td><td>–</td><td>–</td><td><strong>0.75</strong></td></tr><tr><td><em>Analytical composition</em></td><td></td><td></td><td></td></tr><tr><td>Dry matter (DM, %)</td><td>96</td><td>98</td><td>97</td></tr><tr><td>Crude protein (% DM)</td><td>49</td><td>50</td><td>50</td></tr><tr><td>Total lipid (% DM)</td><td>23</td><td>22</td><td>23</td></tr><tr><td>Gross energy (kJ/g DM)</td><td>25</td><td>25</td><td>25</td></tr><tr><td>Ash (% DM)</td><td>6</td><td>6</td><td>6</td></tr><tr><td>Phosphorus (% DM)</td><td>1.2</td><td>1.1</td><td>1.2</td></tr><tr><td>Selenium (mg/kg dry feed) <sup>8</sup></td><td><strong>0.3</strong></td><td><strong>0.8</strong></td><td><strong>0.7</strong></td></tr></tbody></table>

### 4.3. Sampling

The broodstock fish were anaesthetized with benzocaine for stripping and afterwards euthanized by a sharp blow to the head for liver dissection in both males and females. For each individual female, samples of pooled oocytes after stripping (1 g sized samples) and progeny at swim-up fry stage (whole-body fry) killed by an overdose of benzocaine were withdrawn. Moreover, a total of 36 individual swim-up fry livers were randomly dissected on the same day at the Ecology and Fish Population Biology facility in Saint-Pée-sur-Nivelle, France [52], originating from 12 females (n = 4 females per dietary treatment). The three individual livers per female were pooled into a single sample tube for DNA extraction. All collected samples were immediately frozen in liquid nitrogen and stored at -80 °C until further analysis.

### 4.4. Metabolite Analysis

In 0.1 g of pooled whole-body swim-up fry, free amino acids and other N-metabolites were analyzed using the Biochrome Analyzer and post column ninhydrin reaction following deproteinization, as previously described [53]. The aminothiols of the transsulfuration and glutathione pathway including homocysteine, cysteine, γ-glutamyl-cysteine, reduced glutathione and cysteinyl-glycine as well as SAM and SAH were measured by HPLC using one sample extract. First, 0.3 g of broodstock liver tissue, 1 g of pooled oocytes or 1 g whole-body swim-up fry were homogenized with an ultra-turrax in a 20 mM phosphate, 1 mM EDTA (pH = 6.4) buffer. After centrifugation (10,000 g, 15 min, 4 °C), deproteinization of the supernatant was performed using a 10% metaphosphoric acid solution. The protocol for aminothiol analysis was adapted from Toyooka and Imai [54]. Derivatization was performed by adding

<footnote n="1">¹ Plant meals (% diet): 20% wheat gluten (Roquette), 18% corn gluten meal (Inzo), 15% soybean protein concentrate Estril®75 (Sopropêche), 6% soybean meal (Sud-Ouest Aliment), 5% rapeseed meal 00 (Sud-Ouest Aliment), 5% white lupin meal Farilup 500 (Terrena), 3% dehulled pea meal Primatex (Sotexpro), 2% whole wheat (Sud-Ouest Aliment). ² Crystalline amino acids and attractant mixture (% diet): 1.34% L-lysine, 0.3% DL-methionine, 0.5% glucosamine, 0.3% taurine, 0.3% betaine, 0.2% glycine, 0.2% alanine. ³ Soybean lecithin from Louis François and fish oil from Sopropêche. ⁴ Vegetable oils (% diet): 4% rapeseed oil, 2.4% linseed oil, 1.6% palm oil (Daudry). ⁵ Provided as Carophyll®pink (DSM). ⁶ Vitamin and mineral mixture without Se (per kg diet): retinol acetate, 55,000 IU; cholecalciferol, 2,500 IU; DL-α-tocopherol acetate, 50 IU; sodium menadione bisulfate, 10 mg; thiamin-HCl, 1 mg; riboflavin, 4 mg; niacin, 10 mg; D-calcium pantothenate, 20 mg; pyridoxine-HCl, 3 mg; D-biotin, 0.2 mg; folic acid, 1 mg; cyanocobalamin, 10 µg; L-ascorbyl-2-polyphosphate, 50 mg; myo-inositol, 0.3 g; choline, 1 g; CaHPO₄·2H₂O, 33 g; CaCo₃, 2.15 g; Mg(OH)₂, 1.24 g; KCl, 0.9 g; NaCl, 0.4 g; FeSO₄·7H₂O, 0.2 g; ZnSO₄·7H₂O, 40 mg; MnSO₄·H₂O, 30 mg; CuSO₄·5H₂O, 30 mg; NaF, 10 mg; KI, 0.4 mg; CoCl₂·6H₂O, 0.2 mg. All ingredients were diluted with α-cellulose. ⁷ Sodium selenite contained 42% Se (Sigma-Aldrich) and hydroxy-selenomethionine contained 40% Se provided as Selisseo®(Adisseo). ⁸ Total Se was determined using inductively coupled plasma mass spectrometry (ICP MS, Agilent series 7500cx) by Ultra-Trace Analysis Aquitaine (UT2A, Pau, France) according to Vacchina and Dumont [51], with a calculated uncertainty of 15 µg/kg and a limit of quantification of 3 µg/kg.</footnote>
62.5 µL AccQ-Fluor™ borate buffer (Waters, Guyancourt, France), 5 µL 1.55N NaOH and 4.5 mM ABD-F buffer to 25 µL sample aliquot for 20 min at 60 °C. After, the reaction was stopped by addition of 12.5 µL of 1N HCL and cooling at 4 °C for 15 min. Separation was performed using a AccQTag™ column at 40 °C using gradient elution: 0–2 min 97% A, 3% C; 20 min: 96% A, 4% C; 25 min: 20% B, 80% C; 30–35 min: 97% A, 3% C with (A) aqueous solution of AccQTag™ Eluent A; (B) ultra-pure water and (C) methanol. Aminothiols were detected with fluorescence (excitation 385 nm, emission 515 nm). SAM/SAH measurement was adapted from She et al. [55] with separation on a Revolve C18 at 40 °C with the following gradient: 0–10 min 95% A, 5% B; 20 min 30% A, 70% B; 35–45 min 95% A, 5% B with (A) 20mM phosphate buffer with 8mM OSA (pH 2.7, TFA adjusted) and (B) methanol. In 0.1 g of pooled whole-body swim-up fry, pyridoxine, pyridoxal and pyridoxamine were measured by ultra-performance liquid chromatography (UPLC) [56] and vitamin B12 and total folate were analyzed microbiologically using *Lactobacillus delruceckii* ssp. *lactis* and *Lactobacillus rhamnosus*, respectively, as previously described [57].

### 4.5. RNA Extraction and RT-qPCR

The RNA was extracted and analyzed by quantitative RT-qPCR on 0.1 g samples of broodstock liver and a pool of three whole-body swim-up fry, as previously described [6]. The primer sequences are given in Table 5.

<table><caption>**Table 5.** Oligonucleotide primers used to assay mRNA levels by Fluidigm PCR.</caption><thead><tr><td>Gene</td><td>Accession No.</td><td>Forward Primer</td><td>Reverse Primer</td><td>Amplification Size</td></tr></thead><tbody><tr><td>amd1a</td><td>XM_021611778.1</td><td>ccgtaccatcccaaggtttga</td><td>tctgcttgtcgtgtctttgt</td><td>87</td></tr><tr><td>amd1b</td><td>XM_021600287.1</td><td>cagccagattttcccaaacgg</td><td>gcatgctcgttctcccagaa</td><td>108</td></tr><tr><td>bhmt</td><td>FR908041.1</td><td>cagagaagcacggtaactgg</td><td>ttctttgtgctgcatcaggt</td><td>188</td></tr><tr><td> CBS</td><td>NM_001124686.1</td><td>ccacctcaggcaatacaggt</td><td>aacatccaccttctccatgc</td><td>107</td></tr><tr><td> cgl</td><td>EU315111.1</td><td>caccaacccaccatgaaag</td><td>gcgctggaagtaggctgaca</td><td>118</td></tr><tr><td> dnmt1</td><td>XM_021557911.1</td><td>ttgccagaagaggagatgcc</td><td>cccaggtcagcttgcatta</td><td>152</td></tr><tr><td> gnmt</td><td>XM_021585680.1</td><td>ctcaagtcagcgctgaagga</td><td>cactctggtcccctttgaagt</td><td>187</td></tr><tr><td> mtr</td><td>XM_021576690.1</td><td>aatgcaggtctgccaatac</td><td>ctgatgtgtgcaggagtcgt</td><td>137</td></tr><tr><td> sahh</td><td>XM_021609053.1</td><td>atcaaacgggccacagatgt</td><td>tcgtaccttccatggcagc</td><td>167</td></tr><tr><td> β-actin</td><td>AJ438158.1</td><td>gatgggccgaaagacagcta</td><td>tcgtcccgtaggtgacgat</td><td>105</td></tr></tbody></table>

amd1, adenosylmethionine decarboxylase 1; bhmt, betaine-homocysteine S-methyltransferase 1; cbs, cystathionine beta-synthase; cgl, cystathionine gamma-lyase; dnmt1, DNA methyltransferase 1; gnmt, glycine N-methyltransferase; mtr, methionine synthase; sahh, adenosylhomocysteinase.

### 4.6. Statistical Analysis on Metabolic Analysis and Gene Expression Data

Results are given as the mean ± SEM. Statistical analysis was performed using statistical software R (R Core Team). All data were tested for normality and homogeneity. Gene expression data were rank transformed before further analysis. Principle component analysis (PCA) was performed on the free amino acid dataset in search for biological clusters and outliers (R: factoextra [58]). One-way ANOVA was used to identify differences between Se treatments or sex. Tukey's HSD was used as a post hoc test in case a significant difference (*p* < 0.05) was detected.

### 4.7. DNA Extraction, RRBS Library Preparation and Sequencing

DNA extraction on swim-up fry livers was performed using a QIAGEN DNeasy Blood and Tissue Kit (cat. no. 69504), following the manufacturer's instruction. DNA quantity was measured using Qubit fluorometric quantitation (Life Technologies, Carlsbad, California, USA), ensuring that the sample contained a minimum of 200 ng of DNA. The DNA extract was stored at -20 °C before DNA methylation was measured by reduced representation bisulfate sequencing (RRBS) performed at the Biomedical Sequencing Facility BSF in Vienna, Austria.

The RRBS library preparation was performed, as previously described [23], on 100 mg genomic DNA including DNA digestion (Msp1 20 U, 16h at 37 °C), enzymatic adapter ligation (T4 DNA
Ligase rapid), quantification and pooling. Bisulfite conversion was performed using EZ DNA Methylation-Direct Kit D5020, Zymo Research, but conversion reagent was used at 0.9× concentration with incubation for 20 cycles of 1 min at 95 °C, 10 min at 60 °C and a desulphonation time of 30 min to increase the number of CpG nucleotides covered. Enrichment PCR was performed after AMPure XP clean up and library concentrations were quantified with the Qubit Fluorometric Quantitation system (Life Technologies) and size distribution by a Bioanalyzer High Sensitive DNA Kit (Agilent). Sequencing was performed on Illumina HiSeq 3000/4000 instruments. The data have been stored in SRA [59] under the accession number PRJNA629594.

## 4.8. Rainbow Trout Genome and Genomic Annotation

The reference genome data of rainbow trout (Omyk_1.0) were downloaded from the NCBI assembly site (https://www.ncbi.nlm.nih.gov/assembly/ GCF_002163495.1).

For genes with multiple RefSeq sequences, only the longest sequence was kept after eliminating overlapped isoforms. All the CpG sites in the genome were identified and split into four regions—gene body (GB), promoter (P), flanking regions around mRNA (flanks), and intergenic. Gene body was further divided into two sub-regions, intron and exon, whereas promoter was also divided based on the distance from the transcriptional start site (TSS) as P250 (1 bp–250 bp), P1K (251 bp–1000 bp) and P6K (1001 bp–6000 bp). Flanks were defined as a combination of 4K upstream from the 5' end of P6Ks (equivalently 6001–10,000 bp from TSS) and 10K downstream of the 3' end of mRNA. All the regions outside of gene bodies, promoters and flanks were annotated as intergenic. Each CpG site was defined as a unique and non-redundant region or sub-region according to the precedence of exon > intron > P250 > P1K > P6K > flanks > intergenic.

## 4.9. RRBS Data Processing

Illumina2bam tools (1.17.3; https://github.com/wtsi-npg/illumina2bam) were used to de-multiplex pooled samples. SAMtools [60] was used to convert BAM files into FASTQ, before quality check by FastQC (Babraham Institute; https://www.babraham.ac.uk) and MultiQC [61]. Adapters and low-quality reads in the RRBS mode based on Cutadapt [62] were removed with Trim Galore! (Babraham Institute). Long reads were trimmed to 50 bp, and reads were selected by in-house python scripts to keep only those digested by MspI and TaqI.

Reads were aligned to the rainbow trout genome by Bismark [63] with Bowtie 1 [64]. Two Bismark tools, bismark_methylation_extraction and coverage2Cystosine, were used to retrieve methylation calls at CpG sites. Reads were filtered by methylKit tool [65] when either the number of reads was above 99.9th percentile or less than or equal to 10.

Cluster analysis was performed by Rtsne [66] for t-SNE [67], with perplexity = 2 and factoextra [58] for PCA, scree plot and hierarchical clustering with Ward's method.

Prior to differential methylation calculation, the unite function of methylKit was used to form SS:NC and SO:NC with NC as control and SO:SS with SS as control. Methylation differences were calculated by methylKit for all the CpG sites with methylation calls as a percentage and *p*-values by logistic regression. The SLIM method [68] was used to calculate *q*-values. CpG sites with a *q*-value of < 0.01 and ≥ 20% methylation difference were defined as differentially methylated cytosines (DMCs). Genes with at least one DMC in the gene body or promoter region are considered to be differentially methylated genes (DMGs).

In-house R and Python scripts were coordinated in a pipeline by using Snakemake [69].

## 4.10. Functional Annotation and Statistical Analysis of DMGs

To find Kyoto Encyclopedia of Genes and Genomes (KEGG) [70] orthologues that correspond to rainbow trout genes, the results of BLASTKoala, GhostKoala [71] and KEGG Automatic Annotation Server (KAAS) [72] were merged. The precedence of BLASTKoala > GhostKoala > KASS was applied when conflicting annotation occurred. A total of 22501 orthologues along with 168 KEGG pathways
were identified. Over representation analysis (ORA) on KEGG pathways and Gene Ontology (GO) terms [73] was performed on DMGs by the R package clusterProfiler [74,75].

The Wilcoxon signed-rank test (Wilcox) was used to test the differences of methylation rates between two groups in a pair-wise manner for KEGG pathways.

A bootstrap version of the Kolmogorov–Smirnov test (KS-boot; the number of iterations: 1000) was used to test the methylation differences that are associated with a KEGG pathway against the methylation differences of the whole CpG sites in a region. All three methods of enrichment analysis were performed for all the defined regions, and the *p*-values were adjusted by the Benjamini–Hochberg procedure.

## 5. Conclusions

Our results demonstrate that in rainbow trout, parental Se nutrition decreased transsulfuration and modified the methionine cycle, as summarized in Figure 7. A decrease in the methyl donor SAM was noticed in parental fish and their offspring by Se supplementation. In the offspring, significant changes in the DNA methylation pattern were identified, especially for genes related to signal transmission and immune function, by parental Se supplementation with organic and inorganic Se forms. It could be suspected that such epigenetic changes might persist during subsequent growth and development of the fish, leading to long-term molecular and metabolic alterations in the progeny, which deserves further investigation.

**Author Contributions:** Conceptualization, S.F.-D., S.J.K. and P.A.J.P.; methodology, C.H. and K.H.S.; software, T.S.; formal analysis, P.W. and T.S.; investigation, P.W., C.H., P.A.J.P., S.F.-D. and K.H.S.; resources, C.H.; data curation, T.S.; writing—original draft preparation, P.W.; writing—review and editing, S.F.-D. and K.H.S.; visualization, T.S. and P.W.; supervision, B.F., P.A.J.P., S.F.-D. and K.H.S.; project administration, S.J.K., P.A.J.P. and S.F.-D.; funding acquisition, S.J.K., P.A.J.P., S.F.-D. and K.H.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by I-SITE E2S: ENERGY AND ENVIRONMENTAL SOLUTIONS from the UNIVERSITY OF PAU AND PAYS ADOUR, UPPA (contract number 2017-17 and contract number 2018–178), INSTITUTE OF MARINE RESEARCH (ParSel project number 15329), ADISSEO FRANCE SAS (SelGen contract number 22001062 labelled by the Institut Carnot France Futur Elevage IC-F2E).

**Acknowledgments:** The authors wish to thank P. Maunas and N. Turonnet for the care of the fish and F. Terrier, F. Sandres and A. Lanuque for the preparation of diets. We are grateful to M. Parailloux for RNA extraction and qPCR analysis and the IMR technical staff for analytical assistance. We also thank Amelie Nemc, Bekir Ergüner and Christoph Bock at the Biomedical Sequencing Facility at CeMM, Research Centre for Molecular Medicine (Vienna Austria), for RRBS library preparation, sequencing and bi-sulfite read pre-processing.

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

## Abbreviations

<table><tr><td>amd</td><td>Adenosylmethionine decarboxylase</td></tr><tr><td>bhmt</td><td>Betaine-homocysteine S-methyltransferase</td></tr><tr><td>cis</td><td>Cystathionine beta-synthase</td></tr><tr><td>cgl</td><td>Cystathionine gamma-lyase</td></tr><tr><td>CpG</td><td>Cytosine phosphate guanine</td></tr><tr><td>DMCs</td><td>Differentially methylated cytosines</td></tr><tr><td>DMGs</td><td>Differentially methylated genes</td></tr><tr><td>DNMT</td><td>DNA methyltransferase</td></tr><tr><td>gnmt</td><td>Glycine N-methyltransferase</td></tr><tr><td>mtr</td><td>Methionine synthase</td></tr><tr><td>NC</td><td>Non-supplemented control treatment</td></tr><tr><td>OH-SeMet</td><td>Hydroxy-selenomethionine</td></tr><tr><td>RRBS</td><td>Representation bisulfite sequencing</td></tr><tr><td>SAH</td><td>S-adenosylhomocysteine</td></tr><tr><td>sahh</td><td>Adenosylhomocystenase</td></tr><tr><td>SAM</td><td>S-adenosylmethionine</td></tr><tr><td>SO</td><td>OH-SeMet treatment</td></tr><tr><td>SS</td><td>Sodium selenite treatment</td></tr></table>
Appendix A

<table><caption>Table A1. Alignment statistics of 12 RRBS samples.</caption><thead><tr><th>No</th><th>Name</th><th>Diet</th><th>Total Reads</th><th>Uniquely Mapped</th><th>(%)</th><th>Multi-Mapped</th><th>(%)</th><th>Non-Mapped</th><th>(%)</th></tr></thead><tbody><tr><td>1</td><td>PSL1</td><td>NC</td><td>53 987 762</td><td>26 310 305</td><td>48.7</td><td>19 820 505</td><td>36.7</td><td>7 856<br>952</td><td>14.6</td></tr><tr><td>2</td><td>PSL2</td><td>NC</td><td>39 364 790</td><td>19 543 117</td><td>49.6</td><td>14 534 198</td><td>36.9</td><td>5 287<br>475</td><td>13.4</td></tr><tr><td>3</td><td>PSL3</td><td>NC</td><td>54 833 812</td><td>27 173 708</td><td>49.6</td><td>20 233 149</td><td>36.9</td><td>7 426<br>955</td><td>13.5</td></tr><tr><td>4</td><td>PSL4</td><td>NC</td><td>63 481 758</td><td>30 305 331</td><td>47.7</td><td>24 622 757</td><td>38.8</td><td>8 553<br>670</td><td>13.5</td></tr><tr><td>5</td><td>PSL5</td><td>SS</td><td>59 570 698</td><td>27 527 855</td><td>46.2</td><td>23 489 555</td><td>39.4</td><td>8 553<br>288</td><td>14.4</td></tr><tr><td>6</td><td>PSL6</td><td>SS</td><td>45 143 094</td><td>21 814 078</td><td>48.3</td><td>17 312 386</td><td>38.4</td><td>6 016<br>630</td><td>13.3</td></tr><tr><td>7</td><td>PSL7</td><td>SS</td><td>67 495 806</td><td>31 151 354</td><td>46.2</td><td>27 389 024</td><td>40.6</td><td>8 955<br>428</td><td>13.3</td></tr><tr><td>8</td><td>PSL8</td><td>SO</td><td>54 472 443</td><td>25 214 607</td><td>46.3</td><td>22 345 586</td><td>41.0</td><td>6 912<br>250</td><td>12.7</td></tr><tr><td>9</td><td>PSL9</td><td>SO</td><td>69 558 057</td><td>33 562 782</td><td>48.3</td><td>26 840 717</td><td>38.6</td><td>9 154<br>558</td><td>13.2</td></tr><tr><td>10</td><td>PSL10</td><td>SO</td><td>79 827 373</td><td>38 293 256</td><td>48.0</td><td>31 322 516</td><td>39.2</td><td>10 211<br>601</td><td>12.8</td></tr><tr><td>11</td><td>PSL11</td><td>SO</td><td>55 391 941</td><td>25 516 081</td><td>46.1</td><td>21 796 848</td><td>39.4</td><td>8 079<br>012</td><td>14.6</td></tr><tr><td>12</td><td>PSL12</td><td>SS</td><td>51 950 394</td><td>23 789 256</td><td>45.8</td><td>20 152 309</td><td>38.8</td><td>8 008<br>829</td><td>15.4</td></tr></tbody></table>

A CpGs with high variances (>50% quantile)

<figure><img src="image_15.png" /></figure>

B CpGs with high variances (>95% quantile)

<figure><img src="image_16.png" /></figure>

C Scree plot (>50% quantile)

<figure><img src="image_17.png" /></figure>

D Scree plot (>95% quantile)

<figure><img src="image_18.png" /></figure>

E Pairwise correlations (Pearson)

<figure><img src="image_19.png" /></figure>

F Hierarchical clustering (k = 2)

<figure><img src="image_20.png" /></figure>

**Figure A1.** Clustering of the 12 RRBS samples by different algorithms. (A + B) PCA biplot by dietary treatment. (C + D) Scree plot with percentage of explained variance within the top 10 dimensions of the PCA. (E) Heatmap with sample-sample distance calculated by normalized Pearson's correlation coefficient in a range between 0 and 1, with *d* = 0 as *r* = 1 and *d* = 1 as *r* = -1. (F) Dendrogram with hierarchical clustering.
Appendix B

<figure><img src="image_21.png" /><figcaption>B: SO:NC</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>C: SO:SS</figcaption></figure>

<figure><img src="image_23.png" /></figure>

**Figure A2.** Enriched KEGG pathways in SS:NC (A), SO:NC (B) and SO:SS (C), with the number of differentially methylated genes (DMGs).
## Appendix C

<table><caption>Table A2. Top five DMGs in the dataset SS:NC, with the highest number of CpG in each sub-region.</caption><thead><tr><th></th><th>Gene Symbol</th><th>Gene ID</th><th>Gene Name</th><th>Total DMC</th><th>DMC in the Region</th><th>Hyper-/Hypomethylated CpG</th></tr></thead><tbody><tr><td rowspan="5">Exon</td><td>LOC110496419</td><td>110496419</td><td>Glycoprotein endo-alpha-1,2-mannosidase-like protein</td><td>10</td><td>10</td><td>0/10</td></tr><tr><td>LOC110503414</td><td>110503414</td><td>TCDD-inducible poly [ADP-ribose] polymerase-like</td><td>10</td><td>8</td><td>8/0</td></tr><tr><td>LOC110497587</td><td>110497587</td><td>SAM and SH3 domain-containing protein 1-like, transcript variant X2</td><td>10</td><td>7</td><td>0/7</td></tr><tr><td>LOC110520294</td><td>110520294</td><td>Von Willebrand factor C domain-containing protein 2-like, transcript variant X1</td><td>7</td><td>7</td><td>5/2</td></tr><tr><td>LOC110529233</td><td>110529233</td><td>MAM domain-containing glycosylphosphatidylinositol anchor protein 2-like</td><td>6</td><td>5</td><td>5/0</td></tr><tr><td rowspan="4">Intron</td><td>CTNNA2 alk</td><td>110534326 110506268</td><td>Catenin alpha-2 ALK receptor tyrosine kinase</td><td>18 17</td><td>18 17</td><td>11/7 13/4</td></tr><tr><td>Isamp</td><td>110507545</td><td>Limbic system-associated membrane protein, transcript variant X5</td><td>15</td><td>15</td><td>9/6</td></tr><tr><td>NRXN2-like</td><td>110531840</td><td>Neurexin-2-like</td><td>16</td><td>12</td><td>5/7</td></tr><tr><td>CDH4-like</td><td>110532012</td><td>Cadherin-4-like</td><td>13</td><td>12</td><td>5/7</td></tr><tr><td rowspan="5">P250</td><td>LOC110498119</td><td>110498119</td><td>Radical S-adenosyl methionine domain-containing protein 2-like</td><td>5</td><td>5</td><td>0/5</td></tr><tr><td>LOC110493563</td><td>110493563</td><td>Septin-9-like</td><td>4</td><td>2</td><td>0/2</td></tr><tr><td>COX6A2</td><td>100335037</td><td>Cytochrome c oxidase subunit VIa</td><td>3</td><td>2</td><td>0/2</td></tr><tr><td>TIMP2-like</td><td>110487076</td><td>Metalloproteinase inhibitor 2-like</td><td>2</td><td>2</td><td>0/2</td></tr><tr><td>LOC110490026</td><td>110490026</td><td>Mitogen-activated protein kinase kinase kinase 8-like</td><td>2</td><td>2</td><td>2/0</td></tr><tr><td rowspan="6">P1K</td><td>LOC110534101</td><td>110534101</td><td>Matrin-3-like</td><td>4</td><td>4</td><td>4/0</td></tr><tr><td>LOC110489756</td><td>110489756</td><td>Transmembrane protein 14C-like</td><td>4</td><td>3</td><td>0/3</td></tr><tr><td>LOC110486159</td><td>110486159</td><td>Proline-rich protein 15-like protein A</td><td>3</td><td>3</td><td>2/1</td></tr><tr><td>TNFR11B-like</td><td>110506163</td><td>Tumor necrosis factor receptor superfamily member 11B-like</td><td>3</td><td>3</td><td>0/3</td></tr><tr><td>TPPP3X2-like</td><td>110518569</td><td>Tubulin polymerization-promoting protein family Member 3-like, transcript variant X2</td><td>3</td><td>3</td><td>0/3</td></tr><tr><td>LOC110503024</td><td>110503024</td><td>Cytochrome c oxidase subunit 4 isoform 2, mitochondrial-like</td><td>8</td><td>8</td><td>8/0</td></tr><tr><td rowspan="4">P6K</td><td>LOC110523211</td><td>110523211</td><td>Oocyte zinc finger protein X1COF6-like, transcript Variant X2</td><td>8</td><td>8</td><td>0/8</td></tr><tr><td>LOC110498688</td><td>110498688</td><td>Fatty acid-binding protein, liver-type-like</td><td>7</td><td>7</td><td>7/0</td></tr><tr><td>LOC110505815</td><td>110505815</td><td>Gamma-aminobutyric acid receptor subunit rho-2-like</td><td>6</td><td>5</td><td>0/5</td></tr><tr><td>LOC110503024</td><td>110503024</td><td>Ras-related C3 botulinum toxin substrate 3-like</td><td>5</td><td>4</td><td>0/4</td></tr></tbody></table>

Genes in bold are commonly highly affected also in SO:NC and SO:SS.

<table><caption>Table A3. Top five DMGs in the dataset SO:NC, with the highest number of CpG in each sub-region.</caption><thead><tr><th></th><th>Gene Symbol</th><th>Gene ID</th><th>Gene Name</th><th>Total DMC</th><th>DMC in the Region</th><th>Hyper-/Hypomethylated CpG</th></tr></thead><tbody><tr><td rowspan="5">Exon</td><td>LOC110490066</td><td>110490066</td><td>E3 ubiquitin-protein ligase rifylin-like, transcript variant X2 von Willebrand factor C</td><td>6</td><td>6</td><td>5/1</td></tr><tr><td>LOC110520294</td><td>110520294</td><td>domain-containing protein 2-like, transcript variant X1</td><td>6</td><td>6</td><td>3/3</td></tr><tr><td>LOC110531157</td><td>110531157</td><td>Serine/threonine-protein kinase WNK2-like</td><td>7</td><td>5</td><td>0/5</td></tr><tr><td>LOC110497587</td><td>110497587</td><td>SAM and SH3 domain-containing protein 1-like, transcript variant X2</td><td>6</td><td>5</td><td>1/4</td></tr><tr><td>LOC110506316</td><td>110506316</td><td>Muscarinic acetylcholine receptor M4-like, transcript variant X1</td><td>6</td><td>5</td><td>0/5</td></tr></tbody></table>
<table><caption>Table A3. *Cont.*</caption><thead><tr><th>Gene Symbol</th><th>Gene ID</th><th>Gene Name</th><th>Total DMC</th><th>DMC in the Region</th><th>Hyper-/Hypomethylated CpG</th></tr></thead><tbody><tr><td>lsamp</td><td>110507545</td><td>Limbic system-associated membrane protein, transcript variant X5</td><td>24</td><td>24</td><td>9/15</td></tr><tr><td rowspan="4">Intron</td><td>LOC110505581</td><td>Placenta growth factor-like</td><td>15</td><td>15</td><td>9/6</td></tr><tr><td>LOC110501635</td><td>Serine/threonine-protein kinase BRSK2-like</td><td>14</td><td>14</td><td>10/4</td></tr><tr><td>LOC110506270</td><td>Protein kinase C-binding protein NELL1-like, transcript variant X1</td><td>13</td><td>13</td><td>5/8</td></tr><tr><td>Catenin alpha-2</td><td>Catenin alpha-2</td><td>13</td><td>13</td><td>6/7</td></tr><tr><td rowspan="5">P250</td><td>LOC110501919</td><td>Alkyldihydroxyacetonephosp hatesynthase, peroxisomal-like, transcript variant X2</td><td>5</td><td>3</td><td>2/1</td></tr><tr><td>LOC110494831</td><td>Complement C1q-like protein 2</td><td>4</td><td>3</td><td>0/3</td></tr><tr><td>LOC110508922</td><td>MARVEL domain-containing protein 2-like, transcript variant X2</td><td>4</td><td>3</td><td>3/0</td></tr><tr><td>CRIP2-like</td><td>Cysteine-rich protein 2-like</td><td>3</td><td>3</td><td>3/0</td></tr><tr><td>LOC110497531</td><td>Uncharacterized LOC110497531, transcript variant X2</td><td>4</td><td>2</td><td>2/0</td></tr><tr><td rowspan="5">P1K</td><td>LOC110493345</td><td>Gastrula zinc finger protein</td><td>4</td><td>4</td><td>4/0</td></tr><tr><td>LOC110521247</td><td>XICGF17.1-like, transcript variant X1</td><td>4</td><td>4</td><td>4/0</td></tr><tr><td>LOC110533598</td><td>Lactadherin-like, transcript variant X1</td><td>4</td><td>4</td><td>0/4</td></tr><tr><td>LOC110498119</td><td>Radical S-adenosyl methionine domain-containing protein 2-like</td><td>3</td><td>3</td><td>3/0</td></tr><tr><td>NFATC3-like</td><td>Nuclear factor of activated T-cells, cytoplasmic 3-like</td><td>3</td><td>2</td><td>0/2</td></tr><tr><td rowspan="5">P6K</td><td>LOC110505815</td><td>Gamma-aminobutyric acid receptor subunit rho-2-like</td><td>8</td><td>5</td><td>4/1</td></tr><tr><td>LOC110519930</td><td>Uncharacterized LOC110519930, transcript variant X2</td><td>6</td><td>5</td><td>0/5</td></tr><tr><td>LOC110520086</td><td>Collagen alpha-1(XXVIII) chain-like</td><td>5</td><td>5</td><td>3/2</td></tr><tr><td>tafl</td><td>TATA-box binding protein associated Factor 6 like, transcript variant X1</td><td>5</td><td>5</td><td>5/0</td></tr><tr><td>LOC110537362</td><td>Glutamate receptor 3, transcript variant X3</td><td>5</td><td>5</td><td>5/0</td></tr></tbody></table>

Genes in bold commonly highly affected also in SS:NC and SO:SS.

<table><caption>Table A4. Top five DMGs in the dataset SO:SS, with the highest number of CpG in each sub-region.</caption><thead><tr><th>Gene Symbol</th><th>Gene ID</th><th>Gene Name</th><th>Total DMC</th><th>DMC in the Region</th><th>Hyper-/Hypomethylated CpG</th></tr></thead><tbody><tr><td>CDH2-like</td><td>110506386</td><td>Neural-cadherin-like</td><td>10</td><td>10</td><td>7/3</td></tr><tr><td rowspan="4">Exon</td><td>PLEKHG7-like</td><td>Pleckstrin homology domain-Containing family G member 7-like</td><td>9</td><td>9</td><td>0/9</td></tr><tr><td>NRXN2-like</td><td>Neurexin-2-like</td><td>21</td><td>8</td><td>2/6</td></tr><tr><td>LOC110496419</td><td>endo-alpha-1,2-mannosidase-like protein</td><td>8</td><td>8</td><td>8/0</td></tr><tr><td>LOC110496815</td><td>Glutamate receptor ionotropic, kainate5-like</td><td>8</td><td>7</td><td>7/0</td></tr><tr><td>lsamp</td><td>110507545</td><td>Limbic system-associated membrane protein, transcript variant X5</td><td>22</td><td>22</td><td>9/13</td></tr><tr><td rowspan="4">Intron</td><td>LOC110500600</td><td>Adhesion G protein-coupled receptor L3-like, transcript variant X5</td><td>23</td><td>20</td><td>17/3</td></tr><tr><td>FBXL17X1</td><td>F-box and leucine rich repeat protein 17, transcript variant X2</td><td>18</td><td>16</td><td>13/3</td></tr><tr><td>LOC110535694</td><td>Glutamate receptor ionotropic, delta-2, transcript variant X3</td><td>17</td><td>15</td><td>5/10</td></tr><tr><td>ZNF407-like</td><td>Zinc finger protein 407-like</td><td>16</td><td>15</td><td>3/12</td></tr></tbody></table>
<table><caption>Table A4. *Cont.*</caption><thead><tr><th></th><th>Gene Symbol</th><th>Gene ID</th><th>Gene Name</th><th>Total DMC</th><th>DMC in the Region</th><th>Hyper-/Hypomethylated CpG</th></tr></thead><tbody><tr><td rowspan="5">P250</td><td>LOC110488021</td><td>110488021</td><td>Calcitonin gene-related peptide type 1 receptor-like</td><td>8</td><td>7</td><td>7/5</td></tr><tr><td>GATA 2-like</td><td>110494514</td><td>GATA-binding factor 2-like</td><td>7</td><td>6</td><td>6/0</td></tr><tr><td>LOC110498119</td><td>110498119</td><td><strong>Radical S-adenosyl methionine domain-containing protein 2-like</strong></td><td><strong>4</strong></td><td><strong>3</strong></td><td><strong>3/0</strong></td></tr><tr><td>LOC110508922</td><td>110508922</td><td>MARVEL domain-containing protein 2-like, transcript variant X2</td><td>3</td><td>3</td><td>3/0</td></tr><tr><td>CRIP2-like</td><td>110505831</td><td>Cysteine-rich protein 2-like</td><td>4</td><td>2</td><td>2/0</td></tr><tr><td rowspan="5">P1K</td><td>LOC110527134</td><td>110527134</td><td>Methyltransferase-like protein 7A</td><td>5</td><td>5</td><td>0/5</td></tr><tr><td>SEMA5A-like</td><td>110489566</td><td>Semaphorin-5A-like</td><td>6</td><td>4</td><td>0/4</td></tr><tr><td>LOC110493345</td><td>110493345</td><td>Gastrula zinc finger protein</td><td>4</td><td>4</td><td>4/0</td></tr><tr><td>LOC100135939</td><td>100135939</td><td>XICGF17.1-like, transcript variant X1</td><td>4</td><td>4</td><td>1/3</td></tr><tr><td>GRTP1a-like</td><td>110527156</td><td>Proteoglycan 4, transcript variant X2</td><td>4</td><td>4</td><td>0/4</td></tr><tr><td rowspan="5">P6K</td><td>LOC110523211</td><td>110523211</td><td>Oocyte zinc finger protein</td><td>8</td><td>8</td><td>8/0</td></tr><tr><td>MARVELD2-like</td><td>110523471</td><td>XICOF6-like, transcript variant X2</td><td>8</td><td>8</td><td>8/0</td></tr><tr><td>MED12-like</td><td>110488993</td><td>MED12-like, transcription subunit 12-like</td><td>9</td><td>7</td><td>7/0</td></tr><tr><td>LOC110522593</td><td>110522593</td><td>F-box only protein 31-like, transcript variant X1</td><td>7</td><td>7</td><td>0/7</td></tr><tr><td>LOC110505815</td><td>110505815</td><td><strong>Gamma-aminobutyric acid receptor subunit rho-2-like</strong></td><td><strong>6</strong></td><td><strong>6</strong></td><td><strong>6/0</strong></td></tr></tbody></table>

Genes in bold are commonly highly affected also in SS:NC and SO:NC.

## References

1. Rayman, M.P. Selenium intake, status, and health: A complex relationship. *Hormones* **2020**, *19*, 9–14.
[CrossRef] [PubMed]

2. Mariotti, M.; Ridge, P.G.; Zhang, Y.; Lobanov, A.V.; Pringle, T.H.; Guigo, R.; Hatfield, D.L.; Gladyshev, V.N. Composition and evolution of the vertebrate and mammalian selenoproteomes. *PLoS ONE* **2012**, *7*, e33066.
[CrossRef] [PubMed]

3. Gatlin, D.M., III; Barrows, F.T.; Brown, P.; Dabrowski, K.; Gaylord, T.G.; Hardy, R.W.; Herman, E.; Hu, G.; Krogdahl, Å.; Nelson, R.; et al. Expanding the utilization of sustainable plant products in aquafeeds: A review. *Aquac. Res.* **2007**, *38*, 551–579.
[CrossRef]

4. Glencross, B.D.; Baily, J.; Berntssen, M.H.G.; Hardy, R.; MacKenzie, S.; Tocher, D.R. Risk assessment of the use of alternative animal and plant raw material resources in aquaculture feeds. *Rev. Aquac.* **2020**, *12*, 703–758.
[CrossRef]

5. Fontagné-Dicharry, S.; Godin, S.; Liu, H.; Prabhu, P.A.J.; Bouyssiere, B.; Bueno, M.; Tacon, P.; Médale, F.; Kaushik, S.J. Influence of the forms and levels of dietary selenium on antioxidant status and oxidative stress-related parameters in rainbow trout (*Oncorhynchus mykiss*) fry. *Br. J. Nutr.* **2015**, *113*, 1876–1887.
[CrossRef] [PubMed]

6. Wischhusen, P.; Parailloux, M.; Geraert, P.-A.; Briens, M.; Bueno, M.; Mounicou, S.; Bouyssiere, B.; Antony Jesu Prabhu, P.; Kaushik, S.J.; Fauconneau, B.; et al. Effect of dietary selenium in rainbow trout (*Oncorhynchus mykiss*) broodstock on antioxidant status, its parental transfer and oxidative status in the progeny. *Aquaculture* **2019**, *507*, 126–138.
[CrossRef]

7. Steinbrenner, H.; Speckmann, B.; Klotz, L.-O. Selenoproteins: Antioxidant selenoenzymes and beyond. *Arch. Biochem. Biophys.* **2016**, *595*, 113–119.
[CrossRef] [PubMed]

8. Hill, K.E.; Burk, R.F. Effect of selenium deficiency and vitamin E deficiency on glutathione metabolism in isolated rat hepatocytes. *J. Biol. Chem.* **1982**, *257*, 10668–10672.
[PubMed]

9. Hill, K.E.; Burk, R.F. Effect of selenium deficiency on the disposition of plasma glutathione. *Arch. Biochem. Biophys.* **1985**, *240*, 166–171.
[CrossRef]

10. Lu, S.C. Regulation of glutathione synthesis. *Mol. Asp. Med.* **2009**, *30*, 42–59.
[CrossRef]

11. Bunk, M.J.; Combs, G.F. Evidence for an impairment in the conversion of methionine to cysteine in the selenium-deficient chick. *Proc. Soc. Exp. Biol. Med.* **1981**, *167*, 87–93.
[CrossRef] [PubMed]

12. Halpin, K.M.; Baker, D.H. Selenium deficiency and transsulfuration in the chick. *J. Nutr.* **1984**, *114*, 606–612.
[CrossRef] [PubMed]
13. Dalto, D.; Matte, J.-J. Pyridoxine (vitamin B6) and the glutathione peroxidase system; a link between one-carbon metabolism and antioxidation. *Nutrients* **2017**, *9*, 189. [CrossRef] [PubMed]

14. Uthus, E.O.; Ross, S.A. Dietary selenium affects homocysteine metabolism differently in Fisher-344 rats and CD-1 mice. *J. Nutr.* **2007**, *137*, 1132–1136. [CrossRef] [PubMed]

15. Uthus, E.O.; Yokoi, K.; Davis, C.D. Selenium deficiency in Fisher-344 rats decreases plasma and tissue homocysteine concentrations and alters plasma homocysteine and cysteine redox status. *J. Nutr.* **2002**, *132*, 1122–1128. [CrossRef]

16. Speckmann, B.; Schulz, S.; Hiller, F.; Hesse, D.; Schumacher, F.; Kleuser, B.; Geisel, J.; Obeid, R.; Grune, T.; Kipp, A.P. Selenium increases hepatic DNA methylation and modulates one-carbon metabolism in the liver of mice. *J. Nutr. Biochem.* **2017**, *48*, 112–119. [CrossRef] [PubMed]

17. Jones, P.A.; Takai, D. The role of DNA methylation in mammalian epigenetics. *Science* **2001**, *293*, 1068–1070. [CrossRef] [PubMed]

18. Razin, A.; Cedar, H. DNA methylation and gene expression. *Microbiol. Mol. Biol. Rev.* **1991**, *55*, 451–458. [CrossRef]

19. Speckmann, B.; Grune, T. Epigenetic effects of selenium and their implications for health. *Epigenetics* **2015**, *10*, 179–190. [CrossRef]

20. Martín, I.; Gibert, M.J.; Pintos, C.; Noguera, A.; Besalduch, A.; Obrador, A. Oxidative stress in mothers who have conceived fetus with neural tube defects: The role of aminothiols and selenium. *Clin. Nutr.* **2004**, *23*, 507–514. [CrossRef]

21. Xu, J.; Sinclair, K.D. One-carbon metabolism and epigenetic regulation of embryo development. *Reprod. Fertil. Dev.* **2015**, *27*, 667–676. [CrossRef] [PubMed]

22. Skjærven, K.H.; Jakt, L.M.; Dahl, J.A.; Espe, M.; Aanes, H.; Hamre, K.; Fernandes, J.M. Parental vitamin deficiency affects the embryonic gene expression of immune-, lipid transport-and apolipoprotein genes. *Sci. Rep.* **2016**, *6*, 34535. [CrossRef] [PubMed]

23. Skjærven, K.H.; Jakt, L.M.; Fernandes, J.M.O.; Dahl, J.A.; Adam, A.-C.; Klughammer, J.; Bock, C.; Espe, M. Parental micronutrient deficiency distorts liver DNA methylation and expression of lipid genes associated with a fatty-liver-like phenotype in offspring. *Sci. Rep.* **2018**, *8*, 1–16. [CrossRef] [PubMed]

24. Antony Jesu Prabhu, P.; Schrama, J.W.; Kaushik, S.J. Mineral requirements of fish: A systematic review. *Rev. Aquac.* **2016**, *8*, 172–219. [CrossRef]

25. Schrauzer, G.N. Nutritional selenium supplements: Product types, quality, and safety. *J. Am. Coll. Nutr.* **2001**, *20*, 1–4. [CrossRef]

26. Burk, R.F.; Hill, K.E. Regulation of selenium metabolism and transport. *Annu. Rev. Nutr.* **2015**, *35*, 109–134. [CrossRef]

27. Dalto, B.D.; Tsoi, S.; Audet, I.; Dyck, M.K.; Foxcroft, G.R.; Matte, J.J. Gene expression of porcine blastocysts from gilts fed organic or inorganic selenium and pyridoxine. *Reproduction* **2015**, *149*, 31–42. [CrossRef]

28. Dalto, D.B.; Audet, I.; Lapointe, J.; Matte, J.J. The importance of pyridoxine for the impact of the dietary selenium sources on redox balance, embryo development, and reproductive performance in gilts. *J. Trace Elem. Med. Biol.* **2016**, *34*, 79–89. [CrossRef]

29. Weekley, C.M.; Harris, H.H. Which form is that? The importance of selenium speciation and metabolism in the prevention and treatment of disease. *Chem. Soc. Rev.* **2013**, *42*, 8870. [CrossRef]

30. Sbodio, J.I.; Snyder, S.H.; Paul, B.D. Regulators of the transsulfuration pathway. *Br. J. Pharmacol.* **2019**, *176*, 583–593. [CrossRef]

31. Geillinger, K.E.; Rathmann, D.; Köhrle, J.; Fiamoncini, J.; Daniel, H.; Kipp, A.P. Hepatic metabolite profiles in mice with a suboptimal selenium status. *J. Nutr. Biochem.* **2014**, *25*, 914–922. [CrossRef] [PubMed]

32. Wolf, N.M.; Mueller, K.; Hirche, F.; Most, E.; Pallauf, J.; Mueller, A.S. Study of molecular targets influencing homocysteine and cholesterol metabolism in growing rats by manipulation of dietary selenium and methionine concentrations. *Br. J. Nutr.* **2010**, *104*, 520–532. [CrossRef] [PubMed]

33. Uthus, E.O.; Ross, S.A.; Davis, C.D. Differential effects of dietary selenium (Se) and folate on methyl metabolism in liver and colon of rats. *Biol. Trace Elem. Res.* **2006**, *109*, 201–214. [CrossRef]

34. Bakke, A.M.; Tashjian, D.H.; Wang, C.F.; Lee, S.H.; Bai, S.C.; Hung, S.S.O. Competition between selenomethionine and methionine absorption in the intestinal tract of green sturgeon (*Acipenser medirostris*). *Aquat. Toxicol.* **2010**, *96*, 62–69. [CrossRef] [PubMed]
35. Fernandes, A.P.; Wallenberg, M.; Gandin, V.; Misra, S.; Tisato, F.; Marzano, C.; Rigobello, M.P.; Kumar, S.; Björnstedt, M. Methylselenol formed by spontaneous methylation of selenide is a superior selenium substrate to the thioredoxin and glutaredoxin systems. *PLoS ONE* **2012**, *7*, e50727. [CrossRef] [PubMed]

36. Nakamuro, K.; Okuno, T.; Hasegawa, T. Metabolism of selenoamino acids and contribution of selenium methylation to their toxicity. *J. Health Sci.* **2000**, *46*, 418–421. [CrossRef]

37. Davis, C.D.; Uthus, E.O.; Finley, J.W. Dietary selenium and arsenic affect DNA methylation in vitro in Caco-2 cells and in vivo in rat liver and colon. *J. Nutr.* **2000**, *130*, 2903–2909. [CrossRef]

38. King, W.D.; Ho, V.; Dodds, L.; Perkins, S.L.; Casson, R.I.; Massey, T.E. Relationships among biomarkers of one-carbon metabolism. *Mol. Biol. Rep.* **2012**, *39*, 7805–7812. [CrossRef]

39. Fiala, E.S.; Staretz, M.E.; Pandya, G.A.; El-Bayoumy, K.; Hamilton, S.R. Inhibition of DNA cytosine methyltransferase by chemopreventive selenium compounds, determined by an improved assay for DNA cytosine methyltransferase and DNA cytosine methylation. *Carcinogenesis* **1998**, *19*, 597–604. [CrossRef]

40. Zeng, H.; Yan, L.; Cheng, W.-H.; Uthus, E.O. Dietary selenomethionine increases exon-specific DNA methylation of the p53 gene in rat liver and colon mucosa. *J. Nutr.* **2011**, *141*, 1464–1468. [CrossRef]

41. Xiang, N.; Zhao, R.; Song, G.; Zhong, W. Selenite reactivates silenced genes by modifying DNA methylation and histones in prostate cancer cells. *Carcinogenesis* **2008**, *29*, 2175–2181. [CrossRef] [PubMed]

42. Jabłońska, E.; Reszka, E. Chapter Eight—Selenium and Epigenetics in Cancer: Focus on DNA Methylation. In *Advances in Cancer Research*; Tew, K.D., Galli, F., Eds.; Selenium and Selenoproteins in Cancer; Academic Press: Cambridge, MA, USA, 2017; Volume 136, pp. 193–234.

43. Ahmed, M.Y.; Al-Khayat, A.; Al-Murshedi, F.; Al-Futaisi, A.; Chioza, B.A.; Pedro Fernandez-Murray, J.; Self, J.E.; Salter, C.G.; Harlalka, G.V.; Rawlins, L.E.; et al. A mutation of EPT1 (SELENOI) underlies a new disorder of Kennedy pathway phospholipid biosynthesis. *Brain* **2017**, *140*, 547–554. [CrossRef] [PubMed]

44. Horibata, Y.; Elpeleg, O.; Eran, A.; Hirabayashi, Y.; Savitzki, D.; Tal, G.; Mandel, H.; Sugimoto, H. EPT1 (selenoprotein I) is critical for the neural development and maintenance of plasmalogen in humans. *J. Lipid Res.* **2018**, *59*, 1015–1026. [CrossRef] [PubMed]

45. Solovyev, N.D. Importance of selenium and selenoprotein for brain function: From antioxidant protection to neuronal signalling. *J. Inorg. Biochem.* **2015**, *153*, 1–12. [CrossRef] [PubMed]

46. Fanning, T.G.; Hu, W.S.; Cardiff, R.D. Analysis of tissue-specific methylation patterns of mouse mammary tumor virus DNA by two-dimensional Southern blotting. *J. Virol.* **1985**, *54*, 726–730. [CrossRef]

47. Duque-Guimarães, D.E.; Ozanne, S.E. Nutritional programming of insulin resistance: Causes and consequences. *Trends Endocrinol. Metab.* **2013**, *24*, 525–535. [CrossRef] [PubMed]

48. Pacitti, D.; Lawan, M.M.; Feldmann, J.; Sweetman, J.; Wang, T.; Martin, S.A.M.; Secombes, C.J. Impact of selenium supplementation on fish antiviral responses: A whole transcriptomic analysis in rainbow trout (*Oncorhynchus mykiss*) fed supranutritional levels of Sel-Plex®. *BMC Genom.* **2016**, *17*, 116. [CrossRef] [PubMed]

49. Jesu, P.P.A.; Holen, E.; Espe, M.; Silva, M.S.; Holme, M.-H.; Hamre, K.; Lock, E.-J.; Waagbø, R. Dietary selenium required to achieve body homeostasis and attenuate pro-inflammatory responses in Atlantic salmon post-smolt exceeds the present EU legal limit. *Aquaculture* **2020**, *526*, 735413. [CrossRef]

50. Elumalai, P.; Kurian, A.; Lakshmi, S.; Faggio, C.; Esteban, M.A.; Ringø, E. Herbal immunomodulators in aquaculture. *Rev. Fish Sci. Aquac.* **2020**, 1–25. [CrossRef]

51. Vacchina, V.; Dumont, J. Total selenium quantification in biological samples by inductively coupled plasma mass spectrometry (ICP-MS). In *Selenoproteins*; Springer: Berlin, Germany, 2018; pp. 145–152.

52. ECP; INRAE. *Ecology and Fish Population Biology Facility (ECP)—Catalogue des Infrastructures de Recherche; INRAE: Paris, France, 2018.* [CrossRef]

53. Espe, M.; Lemme, A.; Petri, A.; El-Mowafi, A. Can Atlantic salmon (*Salmo salar*) grow on diets devoid of fish meal? *Aquaculture* **2006**, *255*, 255–262. [CrossRef]

54. Toyooka, T.; Imai, K. New fluorogenic reagent having halogenobenzofurazan structure for thiols: 4-(aminosulfonyl)-7-fluoro-2,1,3-benzoxadiazole. *Anal. Chem.* **1984**, *56*, 2461–2464. [CrossRef]

55. She, Q.-B.; Nagao, I.; Hayakawa, T.; Tsuge, H. A simple HPLC method for the determination of S-adenosylmethionine and S-adenosylhomocysteine in rat tissues: The effect of vitamin B6 deficiency on these concentrations in rat liver. *Biochem. Biophys. Res. Commun.* **1994**, *205*, 1748–1754. [CrossRef]

56. Albrektsen, S.; Hagve, T.A.; Lie, Ø. The effect of dietary vitamin B6 on tissue fat contents and lipid composition in livers and gills of Atlantic salmon (*Salmo salar*). *Comp. Biochem. Phys. A* **1994**, *109*, 403–411. [CrossRef]
57. Mæland, A.; Rønnestad, I.; Fyhn, H.J.; Berg, L.; Waagbø, R. Water-soluble vitamins in natural plankton (copepods) during two consecutive spring blooms compared to vitamins in *Artemia franciscana nauplii* and *metanauplii*. *Mar. Biol.* **2000**, *136*, 765–772. [CrossRef]

58. Kassambara, A.; Mundt, F. Factoextra: Extract and Visualize the Results of Multivariate Data Analyses, R Package Version 1.0.6. Available online: https://CRAN.R-project.org/package=factoextra (accessed on 24 July 2020).

59. National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA). Available online: https://www.ncbi.nlm.nih.gov/sra (accessed on 24 July 2020).

60. Li, H.; Handsaker, B.; Wysoker, A.; Fennell, T.; Ruan, J.; Homer, N.; Marth, G.; Abecasis, G.; Durbin, R. The sequence alignment/map format and SAMtools. *Bioinformatics* **2009**, *25*, 2078–2079. [CrossRef] [PubMed]

61. Ewels, P.; Magnusson, M.; Lundin, S.; Käller, M. MultiQC: Summarize analysis results for multiple tools and samples in a single report. *Bioinformatics* **2016**, *32*, 3047–3048. [CrossRef] [PubMed]

62. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. *EMBnet J.* **2011**, *17*, 10–12. [CrossRef]

63. Krueger, F.; Andrews, S.R. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. *Bioinformatics* **2011**, *27*, 1571–1572. [CrossRef] [PubMed]

64. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. *Genome Biol.* **2009**, *10*, R25. [CrossRef]

65. Akalin, A.; Kormaksson, M.; Li, S.; Garrett-Bakelman, F.E.; Figueroa, M.E.; Melnick, A.; Mason, C.E. methylKit: A comprehensive R package for the analysis of genome-wide DNA methylation profiles. *Genome Biol.* **2012**, *13*, R87. [CrossRef]

66. Krijthe, J. T-Distributed Stochastic Neighbor Embedding Using Barnes-Hut, R Package Version 0.15. Available online: https://CRAN.R-project.org/package=Rtsne (accessed on 24 July 2020).

67. van der Maaten, L.; Hinton, G. Visualizing data using t-SNE. *J. Mach. Learn. Res.* **2008**, *9*, 2579–2605.

68. Wang, H.-Q.; Tuominen, L.K.; Tsai, C.-J. SLIM: A sliding linear model for estimating the proportion of true null hypotheses in datasets with dependence structures. *Bioinformatics* **2011**, *27*, 225–231. [CrossRef] [PubMed]

69. Köster, J.; Rahmann, S. Snakemake—A scalable bioinformatics workflow engine. *Bioinformatics* **2012**, *28*, 2520–2522. [CrossRef] [PubMed]

70. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. *Nucleic Acids Res.* **2000**, *28*, 27–30. [CrossRef] [PubMed]

71. Kanehisa, M.; Sato, Y.; Morishima, K. BlastKOALA and GhostKOALA: KEGG tools for functional characterization of genome and metagenome sequences. *J. Mol. Biol.* **2016**, *428*, 726–731. [CrossRef] [PubMed]

72. Moriya, Y.; Itoh, M.; Okuda, S.; Yoshizawa, A.C.; Kanehisa, M. KAAS: An automatic genome annotation and pathway reconstruction server. *Nucleic Acids Res.* **2007**, *35*, W182–W185. [CrossRef]

73. Ashburner, M.; Ball, C.A.; Blake, J.A.; Botstein, D.; Butler, H.; Cherry, J.M.; Davis, A.P.; Dolinski, K.; Dwight, S.S.; Eppig, J.T. Gene ontology: Tool for the unification of biology. *Nat. Genet.* **2000**, *25*, 25–29. [CrossRef]

74. Yu, G.; Wang, L.-G.; Han, Y.; He, Q.-Y. clusterProfiler: An R package for comparing biological themes among gene clusters. *Omics* **2012**, *16*, 284–287. [CrossRef]

75. Yu, G.; Wang, L.-G.; Yan, G.-R.; He, Q.-Y. DOSE: An R/Bioconductor package for disease ontology semantic and enrichment analysis. *Bioinformatics* **2015**, *31*, 608–609. [CrossRef]

<figure><img src="image_24.png" /></figure>