Large-Scale Multi-Omics Genome-Wide Association Studies (Mo-GWAS): Guidelines for Sample Preparation and Normalization

1. Experimental design and plant cultivation

NOTE: Set up the experiment depending on the experimental hypothesis, e.g., using a large-scale GWAS population decreases the necessity of multiple replicates, as statistical testing will be performed based on the haplotypes of all the individual SNPs instead of the accession. In contrast, multiple replicates are indispensable in other experimental approaches. The following points must be considered while preparing the experiment.

1. Include enough biological replicates, depending on the experimental hypothesis.
2. Randomize the biological replicates block-wise to reduce local environmental bias during cultivation, e.g., greenhouse, field.
3. Ensure proper maintenance of the plant during growth. Treat plants homogeneously to reduce bias.

2. Preparation of biological plant material

1. Harvest preparation
   1. Label harvesting tubes (20 mL) containing two 5 mm and two 8 mm diameter metal beads for homogenizing. Fill up a dewar with liquid nitrogen.
      NOTE: Plants should be in the vegetative stage for fresh leaf and root tissue harvesting.
2. Harvest biological samples by flash-freezing in liquid nitrogen. Harvest as quickly as possible to exclude the influence of circadian oscillation on metabolism during prolonged harvesting durations^12,13. Store the harvested fresh leaf and root tissues for further processing at -80 °C.
   NOTE: Leaf-cutting to flash-freezing should not take longer than a few seconds as after leaf cleavage, active biological processes would alter metabolic profiles due to wounding. For roots, preclean the roots by washing with water before flash-freezing in liquid nitrogen. Excess water on the root surface should be soaked up with paper tissue. Dried seeds can be stored at room temperature; no freezing in liquid nitrogen is required.
3. Grind the tissue using a tissue mixer mill.
   1. Precool the tube holders in liquid nitrogen for a couple of minutes to maintain a low temperature while grinding the tissue.
   2. Transport the biological samples in a nitrogen-containing dewar after taking them out from the -80 °C freezer.
   3. Grind the tissues to obtain homogeneous powder; use 25 Hz for 1 min and repeat after freezing in liquid nitrogen if the tissue is not homogeneously ground.
4. For grinding dried seeds, place the seeds in a grinding jar with a 15 mm diameter metal bead. Use the same frequency and time as mentioned in 2.3.3.
   NOTE: Clean and precooled mortars and pestles can be used if a tissue mixer mill is unavailable.
5. Precool labeled 2 mL safe-lock microcentrifuge tubes. Weigh 50 mg with an error of ±5 mg of fresh plant material by using an analytical scale. Precool the tools used for transferring plant material in liquid nitrogen. Ensure that plant material stays frozen during the weighing process.
   NOTE: Do not expose fresh plant material too long to room temperature as biological processes are activated by increasing temperature, altering metabolic profiles¹⁴.
6. Generate additional quality control (QC) samples by pooling a proportion of each sample and weighing 50 mg with an error of ±5 mg of pooled fresh plant material into precooled 2 mL safe-lock microcentrifuge tubes.
   NOTE: At least three QC samples are advised for every 60 samples. The QC samples are essential for the downstream correction, normalization, and analyses.

3. Extraction reagents

1. Fresh tissue, e.g., leaves and roots
   NOTE: Sample extraction is based on a previously described protocol¹⁵. This protocol has been modified based on present needs, e.g., multiple tissues, different internal standards, and large-scale experiments. Additionally, all volumes and instrument settings mentioned below are adjusted to in-house analytical units. Protocol users should adjust these according to their analytical unit and biological samples, based on test samples.
   1. Extraction mixture 1 (EM1): methyl tert-butyl ether (MTBE)/methanol (MeOH) (3:1 v/v)
      1. Prepare a mixture of MTBE/MeOH in a 3:1 ratio. For 100 mL of extraction solvent, mix 75 mL of MTBE with 25 mL of MeOH in a clean glass bottle.
         NOTE: Solvents should be handled carefully in the fume hood with proper safety equipment.
      2. Add 45 µL of 1,2-diheptadecanoyl-sn-glycero-3-phosphocholine (1 mg/mL in chloroform) as an internal standard for the UHPLC-MS based lipid analysis, 400 µL of ribitol (1 mg/mL in water) as internal standard for the GC-MS based analysis, and 125 µL of isovitexin (1 mg/mL in MeOH/water (1:1 v/v)) for UHPLC-MS-based metabolite analysis.
         NOTE: The addition of internal standards is necessary for the post-analysis normalization according to analytical needs. As 1 mL of EM1 is needed for each sample, prepare a stock solution according to the experimental sample size, which should be used for the entire experiment. EM1 must be stored at -20 °C. Check for the absence of the used internal standard and overlapping with other compounds in the investigated species. Several standards can be used; the selection of the internal standards in this protocol was based on previous tests using common bean extracts¹⁶.
   2. Extraction mixture 2 (EM2) water/ methanol (MeOH) (3:1 v/v)
      1. For 100 mL EM2, add 75 mL of double-distilled water and 25 mL of MeOH in a clean glass bottle.
      2. Add 500 µL of EM2 per sample, and prepare a stock solution according to the experimental sample size, which should be used for the entire experiment. Store EM2 at 4 °C.
2. Dried seeds
   1. Extraction mixture 3 (EM3) methanol (MeOH)/ water (7:3 v/v)
      1. For 100 mL of EM3, add 70 mL of MeOH and 30 mL of double-distilled water in a clean glass bottle. Prepare 1 mL of EM3 for each sample.
      2. Add 400 µL of ribitol (1 mg/mL in water) as internal standards for the GC-MS-based analysis and 125 µL of Isovitexin (1 mg/mL in MeOH/water (1:1 v/v)) for UHPLC-MS-based metabolite analysis.
         NOTE: Prepare a stock solution according to the experimental sample size and use it for the entire experiment. Store EM3 at 4 °C.

4. Sample extraction

1. Fresh tissue, e.g., leaves and roots
   1. Prepare three 1.5 mL safe-lock microcentrifuge tubes for each sample. Keep EM1 in a -20 °C liquid cooling system. Transfer the fresh samples from the -80 °C freezer to dry ice or liquid nitrogen for transportation. Add 1 mL of precooled EM1 to each 50 mg aliquot and vortex briefly before keeping on ice.
   2. Incubate the samples on an orbital shaker at 800 × g for 10 min at 4 °C.
   3. Sonicate the samples in an ice-cooled sonication bath for 10 min.
   4. Add 500 µL of EM2 using a multichannel pipette to avoid variation in added volumes.
   5. Vortex the samples briefly to mix the extraction mixtures before centrifuging at 11,200 × g for 5 min at 4 °C.
   6. After phase separation occurs, transfer 500 µL of the upper lipid-containing phase to a prelabeled 1.5 mL safe-lock microcentrifuge tube. Remove the rest of the upper phase.
      NOTE: Take care while transferring as this upper phase has a high vapor pressure and tends to leak out from the pipette.
   7. Transfer 150 µL and 300 µL of the lower polar and semi-polar metabolite-containing phases in two 1.5 mL safe-lock microcentrifuge tubes used for GC-MS and UHPLC-MS analysis, respectively.
   8. Concentrate all the extracted fractions by letting the solvents evaporate without heating using a vacuum concentrator and store at -80 °C.
2. Dried seeds
   1. Prepare two 1.5 mL safe-lock microcentrifuge tubes for each sample. Keep EM3 on ice. Place a 5 mm diameter metal bead in the sample aliquots.
   2. Add 1 mL of EM3 in each 50 mg aliquot and homogenize the samples at 25 Hz for 2-3 min before putting them on ice.
   3. Sonicate the samples in an ice-cooled sonication bath for 10 min.
   4. Vortex the samples briefly before centrifuging at 11,200 × g for 5 min at 4 °C.
   5. Transfer 150 µL and 300 µL of the supernatant in two 1.5 mL safe-lock microcentrifuge tubes used for GC-MS and UHPLC-MS analysis, respectively.
   6. Concentrate all extracted fractions by letting the solvents evaporate without heating using a vacuum concentrator and store at -80 °C.
      NOTE: Based on experience, users are advised to perform step 4.2 for semi-polar metabolites and derivatized metabolite analysis in dried seeds. Perform extraction step 4.1 for dried seed lipid analysis.

5. Analysis of lipids using UHPLC-MS

1. Re-suspend the dried lipid fractions in 250 µL of acetonitrile:2-propanol (7:3, vol/vol).
2. Sonicate the lipid phase for 5 min, centrifuge at 11,200 × g for 1 min.
3. Transfer 90 µL of the supernatant to a glass vial for LC-MS.
4. Inject 2 µL of the extracts into the LC-MS.
5. Perform lipid fractionation on a reversed-phase C₈ column held at 60 °C running with a flow of 400 µL/min with gradual changes of eluent A and B as shown in Table 1. Acquire the mass spectra in positive ionization mode with a mass range of 150-1,500 m/z.
6. Include several QC samples in all daily batches and a blank to ensure correction for analytical variation. Randomize samples block-wise in sequential order.

6. Analysis of polar and semi-polar metabolites using UHPLC-MS

1. Re-suspend the dried polar phase in 180 µL of UHPLC-grade methanol: water (1:1 v/v).
2. Sonicate the polar phase for 2 min, centrifuge at 11,200 × g for 1 min.
3. Transfer 90 µL of the supernatant to a glass vial for LC-MS.
4. Inject 3 µL of the extracts into the LC-MS.
5. Perform metabolite fractionation on a reverse phase C₁₈ column held at 40 °C running with a flow of 400 µL/min with gradual changes of eluent A and B as shown in Table 1. Acquire the mass spectra in a mass range of 100-1,500 m/z in a full MS scan and all ion fragmentation (AIF) induced by high-energy collisional dissociation (HCD) of 40 keV.
   NOTE: Use both ionization modes. However, due to limited capacity while running large numbers of samples, run test samples in both ionization modes to determine the preferred ionization mode.
6. Include several QC samples in all daily batches and a blank to ensure correction for analytical variation. Randomize samples block-wise in sequential order.
7. Run a pooled QC in data-dependent MS² in both negative and positive ionization modes. Use the obtained mass spectra in a later step (8.5) for annotation.

7. Analysis of derivatized Metabolites using GC-MS^17,18

NOTE: The analysis of derivatized metabolites is based on a previously described protocol¹⁷. Handle all derivatization reagents in the fume hood. Ensure that N-methyl-N-(trimethylsilyl)trifluoracetamide (MSTFA) does not get in contact with water and humidity.

1. Derivatization reagent 1 (DR1)
   1. Dissolve methoxyamine hydrochloride in pyridine to obtain a concentration of 30 mg/mL of DR1. Use 40 µL of DR1 for each sample. Prepare a stock solution according to the sample size, and store at room temperature.
2. Derivatization reagent 2 (DR2)
   1. Dissolve MSTFA with 20 µL of fatty acid methyl esters (FAMEs) per 1 mL of MSTFA. Use 70 µL of DR2 for each sample. Prepare a stock solution according to the sample size. Store MSTFA at 4 °C and the FAMEs at -20 °C.
      NOTE: FAMEs include methylcaprylate, methyl pelargonate, methylcaprate, methyllaurate, methylmyristate, methylpalmitate, methylstearate, methyleicosanoate, methyldocosanoate, lignoceric acid methyl ester, methylhexacosanoate, methyloctacosanoate, and triacontanoic acid methylester, which are dissolved in CHCl₃ at a concentration of 0.8 µL/mL or 0.4 mg/mL for liquid or solid standards, respectively.
3. Re-dry the pellet from the polar phase (stored at -80 °C) using a vacuum concentrator for 30 min to avoid any interference of H₂O originating during the storage with the solvents used for the downstream derivatization.
4. Add 40 µL of DR1.
5. Shake the samples at 950 × g for 2 h at 37 °C using an orbital shaker, followed by a short spin-down of the liquid.
6. Add 70 µL of DR2.
7. Shake again at 950 × g for 30 min at 37 °C using an orbital shaker.
8. Centrifuge briefly at room temperature before transferring 90 µL into glass vials for GC-MS analysis.
9. Inject 1 µL to GC-MS splitless mode, depending on the metabolite concentrations, with a constant helium carrier gas flow of 2 mL/min. Injection temperature is set to 230 °C using a 30-m MDN-35 capillary column.
   NOTE: Additional information, e.g., temperature gradient, can be found in Table 1. The mass range is set to 70-600 m/z with 20 scans/min. Include split modes to enable the quantification of putative overloading compounds, saving costs and time for extract re-derivatization in such cases.
10. Include several QC samples in all daily batches and a blank to ensure correction for analytical variation. Randomize samples properly block-wise in sequential order.

8. Chromatogram processing and compound annotation

1. Filter chemical noise by defining intensity thresholds. Include all QC samples while processing the chromatograms.
   NOTE: For large-scale data, noise filtering is crucial to decrease computing time and processing power.
2. Align the chromatograms by defining a retention time shift window. Check the chromatograms from each batch to assess the intra- and inter-batch variation.
3. Perform peak detection depending on the peak shape, e.g., height and width for full width at half-maximum (FWHM) calculations.
4. Cluster isotopes to reduce redundant signals and filter out singletons.
   NOTE: See the Table of Materials for details on software used for chromatogram processing. In-depth protocols on how to process chromatograms using various freely available software tools, e.g., MS-DIAL, MetAlign, MzMine, and Xcalibur^19,20,21, are provided.
5. Use the ddMS² data of a pooled QC sample for compound annotation. Assess the molecular structure by determining the monoisotopic mass and observing common neutral losses, known charged aglycones, and different types of cleavages, e.g., homolytic or heterolytic^16,22.
6. For reporting metabolite data, follow the recommendation described in Fernie et al. 2011²³.
   NOTE: Different computational metabolomics approaches can be used to analyze metabolomics data^24,25,26.

9. Normalization of large-scale metabolomics dataset

1. Check the distribution of the internal standard(s) and normalize by correcting for the response of single or multiple internal standards.
2. Correct the peak intensities obtained from the chromatogram over the exact sample weight by dividing the peak intensities by the aliquoted homogenized sample weight from step 2.5.
3. Correct for intensity drift across multi-batch series. Perform QC-based correction methods such as locally estimated scatterplot smoothing (LOESS)²⁷ using R.
   NOTE: Several tools and packages are available to address the drift of the MS performance during the acquisition of the whole batches^28,29.
4. Ensure normal distribution of traits by data transformation, e.g., Box-Cox transformation³⁰ using the boxcox () function from the R package MASS for carrying out GWAS.
5. Perform data scaling, e.g., Pareto scaling, for multivariate analysis to ensure proper weighing of low abundant compounds³¹.
   NOTE: If possible, perform a recovery assay to avoid matrix effects, e.g., ion suppression¹⁴.

10. Genome-wide association studies (GWAS)³²

1. Call single nucleotide polymorphism (SNP) or structural variants (SV) from the sequencing data^33,34.
2. Filter genotypic data for minor allele frequency (MAF) < 5% and missing rate of >10% to avoid low-frequency bias using Tassel³⁵.
3. Calculate best linear unbiased predictions (BLUPs) for each normalized feature over the experimental repetitions to eliminate bias originating from environmental factors (random effects) using the R package Ime4³⁶.
4. Use BLUPs of each feature individually to perform GWAS using the rMVP package in R³⁷.
   NOTE: Each metabolomics feature is viewed here as an individual stand-alone phenotype.
5. While performing GWAS, correct for population structure using principal component analysis (PCA) and identity by state (IBS) or vanRaden to minimize confounding effects. Furthermore, consider using a mixed linear model (MLM) or a multi-locus mixed model (MLMM), as mixed models contain fixed and random effects.

11. QTL detection

1. Check the SNPs showing significant association, taking the Manhattan plots into consideration, for linkage disequilibrium (LD) calculations to determine the underlying genetic region. Perform LD calculations using the R package LD heatmap or Tassel 5.
2. Check the associated SNPs for the effect size over the trait by examining the trait levels for statistical changes between haplotypes to find potential causal SNPs, e.g., SNPs leading to an amino acid change in the protein-coding sequence, which could explain the phenotypic variation.
   NOTE: As SNP-trait associations do not necessarily yield causal association, it is crucial to determine the genomic region. Compound identity by feature annotation can help immensely in finding the right candidate genes in a specific genomic region. We suggest to combine all detected QTL associated with certain compounds in a pleiotropic map to underline the genetic regions³⁸, as shown in Figure 4. For validation of candidate genes, several approaches can be performed (see the discussion).