Micronutrient supplementation affects DNA methylation in male gonads with potential intergenerational epigenetic inheritance involving the embryonic development through glutamate receptor-associated genes

Background DNA methylation has an important role in intergenerational inheritance. An increasing number of studies have reported evidence of germline inheritance of DNA methylation induced by nutritional signals in mammals. Vitamins and minerals as micronutrients contribute to growth performance in vertebrates, including Atlantic salmon (Salmo salar), and also have a role in epigenetics as environmental factors that alter DNA methylation status. It is important to understand whether micronutrients in the paternal diet can influence the offspring through alterations of DNA methylation signatures in male germ cells. Results Here, we show the effect of micronutrient supplementation on DNA methylation profiles in the male gonad through a whole life cycle feeding trial of Atlantic salmon fed three graded levels of micronutrient components. Our results strongly indicate that micronutrient supplementation affects the DNA methylation status of genes associated with cell signalling, synaptic signalling, and embryonic development. In particular, it substantially affects DNA methylation status in the promoter region of a glutamate receptor gene, glutamate receptor ionotropic, NMDA 3A-like (grin3a-like), when the fish are fed both medium and high doses of micronutrients. Furthermore, two transcription factors, histone deacetylase 2 (hdac2) and a zinc finger protein, bind to the hyper-methylated site in the grin3a-like promoter. An estimated function of hdac2 together with a zinc finger indicates that grin3a-like has a potential role in intergenerational epigenetic inheritance and the regulation of embryonic development affected by paternal diet. Conclusions The present study demonstrates alterations of gene expression patterns and DNA methylation signatures in the male gonad when Atlantic salmon are fed different levels of micronutrients. Alterations of gene expression patterns are of great interest because the gonads are supposed to have limited metabolic activities compared to other organs, whereas alterations of DNA methylation signatures are of great importance in the field of nutritional epigenetics because the signatures affected by nutrition could be transferred to the next generation. We provide extensive data resources for future work in the context of potential intergenerational inheritance through the male germline. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-08348-4.

reach 36% protein and 34% lipid (~24 MJ) in the largest pellet size in seawater. The experimental feeds were a low FM (fish meal) / FO (fish oil) formulation (initially 15% and 8% in freshwater, decreasing to 5% and 3% respectively in seawater). Feeds were supplemented with a nutrient package (NP) at one of three inclusion levels to produce three dietary treatments: L1, 100% NP; L2, 200% NP; L3, 400% NP (Table S1), the assumption being that the 100% NP package should contain 100 % of assumed requirement based on the given requirement levels reported for Atlantic salmon at the time (NRC 2011) and modified according to earlier trials as part of the EU-funded ARRAINA project Hemre et al. 2016). Specifically, the NP contained 24 nutrients in total these being vitamins (A, D3, E, K3, C, thiamine, riboflavin, B6, B12, niacin, pantothenic acid, folic acid and biotin), minerals (Ca, Co, I, Se, Fe, Mn, Cu and Zn), crystalline amino acids (L-histidine and taurine) and cholesterol. Total and available phosphorus were fixed in all diets at 12.0 g/kg and 9.0 g/kg respectively, and magnesium at 1.5 g/kg, and were not part of the NP. Pellet size was adjusted according to fish weight (2mm, 3.5mm, 5mm, 7mm, 9mm). All non-oil ingredients were mixed, and pellets produced by extrusion to produce three base pellets that had oil added by vacuum coating. All feeds were produced at the BioMar Tech-Centre (Brande, Denmark).

Micronutrient analysis of experimental diets
Micronutrient analysis were performed by several different technologies and methods as described previously (Taylor et al. 2019). Vitamins were determined by microbiological methods and highperformance liquid chromatography (HPLC), whereas minerals were determined by inductively couple plasma mass spectrometry (ICP-MS).

Sampling and growth measurement
Fish were sampled at start and end of the freshwater phase, and then at approximately 250g, 500g, 1kg, and ~2.5kg in seawater prior to dietary pellet size/formulation changes. In freshwater at each sampling point, 50 fish per tank were anaesthetised (Tricaine/MS222, PHARMAQ, UK), individual weights (±0.1g) and fork lengths (±1.0mm) measured, while in seawater, all fish per pen were counted and individually measured.
For body weights and fork lengths, all fish per pen were individually measured. Following measurement, all fish were allowed to recover in aerated water before returning to their original experimental pens. Fulton's condition factor (K) was calculated using: K = (W / L 3 ) * 100 where W is body weight (g), and L is fork length (mm). Hepatosomatic index (HSI) was calculated as HSI (%) = (liver weight (g) / body weight (g)) * 100. Gonadosomatic-index (GSI) was calculated as GSI (%) = (gonad weight (g) / body weight (g)) * 100.
Following measurement, all fish were allowed to recover in aerated water before returning to their original experimental tanks/pens. Maturation at harvest was determined by assessment of external appearance of secondary sexual characteristics and gonad development (n = 30 per pen). Fish were classified as sexually recruited based on a threshold value of GSI >0.20% or >1.0% for males and females, respectively (Kadri et al. 1997).

DNA and RNA extraction
At the end of the feeding trial, gonad and liver tissue were dissected for RNA and DNA extraction from 6 fish per diet at the same area of each individual, the ventral posterior lobe, and snap frozen in liquid nitrogen and stored at -80°C until further processing. DNA and RNA were extracted from same fish. For each feed group, liver samples were selected from the triplicate feeding tanks for RNA sequencing and RRBS. For both RNA and DNA extraction, tissue samples were homogenized using ceramic beads CK28 and a Precellys 24 homogenizer (Bertin Technologies).
RNA was extracted using the BioRobot EZ1 and EZ1 RNA Universal Tissue kit (Qiagen) and DNase treated with Ambion DNA-free DNA removal kit (Invitrogen, USA) according to the protocols. RNA quantity, which was 2460 ± 652 ng/ml on average, were assessed using NanoDrop ND-1000 Spectrophotometer (Nanodrop Technologies). RNA integrity (RIN), which was 9.55 ± 0.15 on average, were analysed using an Agilent 2100 Bioanalyser (RNA 6000 Nano LabChip kit, Agilent Technologies).
Homogenate from six single tissue samples from each feed group were resuspended in lysis buffer and DNA isolation was performed using the DNeasy Blood &Tissue Kit (Qiagen,Cat. No. #69506) according to the manufacturers protocol, except that after homogenization, tissue for DNA extraction were pre-treated with RNase A (provided by the Qiagen kit, 50ng/µL, 10 min at room temperature) immediately followed by proteinase K treatment (New England biolabs, #8102S 20µg/µL, 1.5 h at 55°C). DNA was eluted in Milli Q H2O. Quantification of double stranded genomic DNA was done using the Qubit High Sensitivity Assay (Life Technologies #Q32854).

RNA-seq library preparation and sequencing
RNA-sequencing (RNA-seq) was performed by the DeepSeq sequencing facility at Nord University, Bodø, Norway. RNA-seq library preparation was completed using an NEBNext Ultra II Directional RNA Library Prep Kit for Illumina using the manufacturer protocol and workflow (New England Biolabs). Libraries were ligated with four primers using NEBNext Multiplex Oligos for Illumina Index Primer Sets 1, 2, 3 and 4. Samples were multiplexed into two pools of nine samples, with multiplexing barcodes ligated to each sample during the PCR amplification step, and further sequenced on the NextSeq500 machine (Illumina).

RRBS library preparation and sequencing
For RRBS, 100 ng of genomic DNA were digested for 6h at 65°C with 20 U TaqI (New England Biolabs) and 6h hours at 37°C with 20 U of MspI (New England Biolabs) in 30 μl of 1x NEB buffer 2. To retain even the smallest fragments and to minimize the loss of material, end preparation and adaptor ligation were performed in a single-tube setup. End fill-in and A-tailing were performed by addition of Klenow Fragment 3′ -> 5′ exo-(New England Biolabs) and dNTP mix (10 mM dATP, 1 mM dCTP, 1 mM dGTP). After ligation to methylated Illumina TruSeq LT v2 adaptors using T4 DNA Ligase rapid (Enzymatics), the libraries were size selected by performing a 0.75x clean-up with AMPure XP beads (Beckman Coulter).
The libraries were pooled based on qPCR data and subjected to bisulfite conversion using the EZ DNA Methylation Direct Kit (Zymo Research) with changes to the manufacturer's protocol: conversion reagent was used at 0.9x concentration, incubation performed for 20 cycles of 1 min at 95°C, 10 min at 60°C and the desulphonation time was extended to 30 min. These changes increase the number of CpG dinucleotides covered, by reducing double-strand break formation in larger library fragments. Bisulfite-converted libraries were enriched KAPA HiFi HS Uracil+ RM (Roche). The minimum number of enrichment cycles was estimated based on a qPCR experiment. After a 1x AMPure XP clean-up, library concentrations were quantified with the Qubit Fluorometric Quantitation system (Life Technologies) and the size distribution was assessed using the Bioanalyzer High Sensitivity DNA Kit (Agilent).

Figure S1. PCA biplots showing gene expression patterns affected by diets.
PCA biplots show the first and second PCA components of all the mapped genes of RNA-seq counts with VST (variance stabilization transformation) for a gonads, b liver and c G&L (gonads and liver) datasets. The area of three diet groups (L1: red, L2: green, L3: blue) and two tissue types (liver: blue, gonads: red) are outlined by convex hulls.

Figure S3. Venn diagrams showing overlapped DEGs between liver and gonads.
Ovals represent three DEG datasets of liver (orange), gonads (blue), and G&L (gonads and liver; green) with the counts of overlapped DEGs for a L2:L1 and b L3:L1 datasets. Non-redundant DEGs are a set of unique DEGs among all the three datasets -liver, gonads, and G&L.  PCA biplots show the first and second PCA components of DNA methylation rates for a six different genomic regions (exon, intron, three promoter regions: P250, P1K and P5K, and flanking regions) and b six regions around TSS (three upstream promoter regions and three downstream exon regions from TSS). The area of three diet groups (L1: red, L2: green, L3: blue) are outlined with the ellipses estimated by the Khachiyan algorithm.

Figure S6. Violin plots and bar plots showing distributions and features of DMCs and DMGs in liver.
a Two violin plots on the top show the distribution of methylation rate differences of all the mapped CpG sites (grey background) as well as hypo-methylated (blue) and hyper-methylated (red) DMCs for L2:L1 and L3:L1 datasets. The label boxes display the numbers of corresponding DMCs. b Two stacked bar plots show the proportions of DMC counts (1 DMC, 2 DMCs and >3DMCs) per DMG for L2:L1 and L3:L1 datasets. The numbers on the bars represent the numbers of corresponding DMGs. The label boxes next to the region names display the total number of DMGs per region. c Two stacked bar plots show the number of DEG:DMGs (DEGs that are also DMGs) for L2:L1 and L3:L1 datasets. The numbers on the bars represent the numbers of corresponding DEGs. The label boxes next to the region names display the ratio of DEGs and DMGs as(#DEG:DMGs)/(#DMGs). All the genes that belong to multiple datasets (for instance, liver and gonads) are categorized at "overlap".

Figure S7. Venn diagrams of overlapped DMCs between liver and gonads.
Ovals represent three DMCs datasets of liver (orange), gonads (blue), and G&L (gonads and liver; green) with the counts of overlapped DMCs for a L2:L1 and b L3:L1 datasets. Non-redundant DMCs are a set of unique DMCs among all the three datasets -liver, gonads, and G&L.

Figure S8. Violin plots and bar plots showing distributions and features of DMCs and DMGs in both gonads and liver (G&L).
a Two violin plots on the top show the distribution of methylation rate differences of all the mapped CpG sites (grey background) as well as hypo-methylated (blue) and hyper-methylated (red) DMCs for L2:L1 and L3:L1 datasets. The label boxes display the numbers of corresponding DMCs. b Two stacked bar plots show the proportions of DMC counts (1 DMC, 2 DMCs and >3DMCs) per DMG for L2:L1 and L3:L1 datasets. The numbers on the bars represent the numbers of corresponding DMGs. The label boxes next to the region names display the total number of DMGs per region. c Two stacked bar plots show the number of DEG:DMGs (DEGs that are also DMGs) for L2:L1 and L3:L1 datasets. The numbers on the bars represent the numbers of corresponding DEGs. The label boxes next to the region names display the ratio of DEGs and DMGs as(#DEG:DMGs)/(#DMGs). All the genes that belong to multiple datasets (for instance, liver and gonads) are categorized at "overlap".

Figure S9. Venn diagrams of overlapped DMGs between liver and gonads.
Ovals represent three DMGs datasets of liver (orange), gonads (blue), and G&L (gonads and liver; green) with the counts of overlapped DMGs for a L2:L1 and b L3:L1 datasets. Non-redundant DMGs are a set of unique DMGs among all the three datasets -liver, gonads, and G&L.      Visualization of genomic data of the nuclear pore complex protein Nup50-like (nup50-like; LOC106591755) locus and its vicinity provides information of genomic features along with methylation differences and methylation rates. a The main track at the bottom shows differences of methylation rates for L2:L1 and L3:L1. The green line indicates a threshold of 25%. b An enlarged view shows the part of the region indicated as a pink rectangle in a. The main track at the bottom shows average methylation rates of L1, L2, and L3. One common DMC between L2:L1 and L3:L1 is highlighted with a pink vertical bar.  Visualization of genomic data of the dnaJ homolog subfamily C member 16-like (dnajc16-like; LOC106582681) locus and its vicinity provides information of genomic features along with methylation differences and methylation rates. a The main track at the bottom shows differences of methylation rates for L2:L1 and L3:L1. Two green lines indicate thresholds of 25%. b An enlarged view shows the part of the region indicated as a pink rectangle in a. The main track at the bottom shows average methylation rates of L1, L2, and L3. Two common DMCs between L2:L1 and L3:L1 are highlighted with pink vertical bars.

Figure S18. Genomic feature view of ddx43 with one common DMC in P250.
Visualization of genomic data of the DEAD (Asp-Glu-Ala-Asp) box polypeptide 43 (ddx43) locus and its vicinity provides information of genomic features along with methylation differences and methylation rates. a The main track at the bottom shows differences of methylation rates for L2:L1 and L3:L1. The green line indicates a threshold of 25%. b An enlarged view shows the part of the region indicated as a pink rectangle in a. The main track at the bottom shows average methylation rates of L1, L2, and L3. One common DMC between L2:L1 and L3:L1 is highlighted with a pink vertical bar.  Visualization of genomic data of the phospholipase C beta 4 (plcb4) locus and its vicinity provides information of genomic features along with methylation differences and methylation rates. a The main track at the bottom shows differences of methylation rates for L2:L1 and L3:L1. Two green lines indicate thresholds of 25%. b An enlarged view shows the part of the region indicated as a pink rectangle in a. The main track at the bottom shows average methylation rates of L1, L2, and L3. Three common DMCs between L2:L1 and L3:L1 are highlighted with pink vertical bars.

Figure S21. Genomic feature view of scrib homolog with two common DMCs in P1K.
Visualization of genomic data of the protein scribble homolog (scrib homolog, LOC106578665) locus and its vicinity provides information of genomic features along with methylation differences and methylation rates. a The main track at the bottom shows differences of methylation rates for L2:L1 and L3:L1. Two green lines indicate thresholds of 25%. b An enlarged view shows the part of the region indicated as a pink rectangle in a. The main track at the bottom shows average methylation rates of L1, L2, and L3. Two common DMCs between L2:L1 and L3:L1 are highlighted with pink vertical bars.  Cholesterol Cholesterol Yes n.a n.a n.a NR † Nutrients added at graded levels to the feeds. ‡ Current NRC, 2011, minimum requirement recommendations determined in a rainbow trout and b Atlantic salmon. Non-numeric values represent; n.a (not analysed), NR (no requirement), NR* (no requirement freshwater), and NT (not tested). Superscripts denote significant differences between diets (p < 0.05, one-way ANOVA).                             Table S41. List of the DEGs that contain at least one DMC in liver.

Fields:
Gene ID Entrez Gene ID Gene symbol Gene symbol Gene description Gene name Log2 FC Output of DESeq2 Adj p-value Output of DESeq2 Region Intron, Exon, P250, P1K, P5K, Flanks Mdiff Methylation difference Q value Logistic regression and SLIM, < 0.01 Table S42. List of the DEGs that contain at least one DMC for G&L.    Table S44. List of the common DMGs between L2:L1 and L3:L1.

Note:
Common DMGs are defined as DMGs that have at least one common DMC. Common DMCs are defined as an intersection of DMCs between L2:L1 and L3:L1. In addition, hyper and hypo methylation must be matched between L2:L1 and L3:L1.