Metagenomic Analysis of Silage

1. Site Location

1. Collect the silage sample from an appropriate site such as a farm. Here, the farm was located in Ballydulea, Co. Cork, Ireland (51°51'58.4"N 8°16'48.7"W).

2. DNA Extraction

NOTE: DNA extraction was performed using a commercial kit following the manufacturer's instructions. A negative control, which contained no sample, was used throughout the library preparation method.

1. Add 100 – 400 mg of sample to 978 µL sodium phosphate buffer and 122 µL soil lysis buffer in the supplied lysis tubes.
2. Homogenize samples by placing the lysis tubes into the homogenizer for 40 s at a speed of 6.0 m/s.
3. Centrifuge lysates at 14,000 x g for 15 min and transfer the supernatant to a clean micro-centrifuge tube containing 250 µL of Protein Precipitate Solution (PPS). Mix the solution by inverting 10 times and centrifuge at 14,000 x g for 5 min.
4. Add the supernatant to 1 mL DNA binding matrix in a clean 15 mL centrifuge tube. Mix the solution by inverting the tube constantly for 3 min. Allow the mixture to settle for 3 min, then discard 500 µL of supernatant. Mix the remaining supernatant.
5. Transfer 600 µL of the suspension to a spin filter and centrifuge at 14,000 x g for 1 min. Discard the filtrate and repeat the process with the remaining suspension.
6. Add 500 µL of wash buffer to the DNA binding matrix within the spin filter, mix by pipetting, then centrifuge at 14,000 x g for 1 min.
7. Discard the filtrate and centrifuge the spin filter again at 14,000 x g for 2 min to ensure all wash buffer is removed. Dry the spin filter at 23 ˚C for 5 min.
8. Pre-warm (70 ˚C) the DNase-free water (DES) and re-suspend the DNA binding matrix in 100 µL of DES within the spin filter. Transfer the spin filter to a clean 1.5 mL micro-centrifuge tube and centrifuge at 14,000 x g for 1 min to elute DNA. Store the purified DNA at -20 ˚C until further analysis is performed.

3. DNA Purification Using DNA Purification Beads

NOTE: Prior to metagenomic library preparation the extracted DNA was purified using purification beads to ensure a pure DNA sample was obtained.

1. Incubate the beads at 23 ˚C for 30 min before use. Add 2 volumes of beads to the DNA sample and incubate the solution at 23 ˚C for 5 min.
2. Place the samples onto a separation magnet for 5 min and then discard the supernatant. Wash the beads twice with 200 µL fresh 80% ethanol (EtOH). Air dry the beads for 10 min.
3. Remove the samples from the separation magnet and add 50 µL of elution buffer (EB), mix by pipetting.
4. Incubate the suspension at 23 ˚C for 5 min, after which place the samples back onto the separation magnet for 3 min.
5. Transfer the supernatant, which contains the DNA, to a clean tube. Discard the beads.
6. Quantify the purified DNA as per section four.

4. Quantification of Purified DNA

NOTE: Purified DNA was quantified using a fluorometer and double-stranded (dsDNA) High Sensitivity (HS) assay kit following the manufacturer's instructions.

1. Prepare a working solution using 199:1 ratio of buffer to reagent.
2. Add 10 µL of each DNA standard to 190 µL of working solution.
3. Add 10 µL of purified DNA to 190 µL of working solution. The final volume should be 200 µL. Incubate standard and DNA samples at 23 ˚C for 2 min.
4. Analyze standards before the DNA samples on the fluorometer using the on-screen instructions.

5. Shotgun Sequencing Library Preparation

NOTE: The shotgun sequencing library was prepared using a commercial library preparation kit using the manufacturer's instructions.

1. Dilute the DNA samples to 0.2 ng/µL using EB. Any sample which is already below this concentration, i.e. the negative control, is left at its current concentration.
2. Mix 5 µL of the purified DNA with 10 µL buffer and 5 µL enzyme mix. Incubate samples at 55 ˚C for 5 min.
3. Add 5 µL of neutralizing buffer and incubate the solution at 23 ˚C for 5 min.
4. Add 5 µL of each of the sample specific sequencing indices and 15 µL of PCR master mix.
5. In a thermocycler, incubate the samples at 72 ˚C for 3 min, 95 ˚C for 30 s, before 12 cycles of 95 ˚C for 10 s, 55 ˚C for 30 s and 72 ˚C for 30 s. Incubate samples finally at 72 ˚C for 5 min.
6. Purify the prepared DNA using the bead purification as before but with a final elution of 30 µL of EB.

6. Library Quantity and Quality Check

NOTE: The quantity and quality of the prepared libraries were assessed using a commercial kit and instrumentation.

1. Incubate the kit components at 23 ˚C for 30 min prior to use.
2. Add 2 µL of DNA to 2 µL of buffer and vortex for 1 min at 2,000 rpm.
3. Spin down the sample to ensure it is at the bottom of the tube.
4. Insert the sample tubes, analysis tape and tips into the instrument, and perform analysis as directed by the software.

7. DNA Sequencing

1. Transfer the prepared and quantified DNA sequencing libraries samples to a sequencing service and sequence using 300 bp paired end sequencing ⁴⁵.

8. Analysis of Raw Sequence Data

NOTE: The commands for each program using a Linux operating system are shown below the protocol step. The pipeline used for sequence data analysis is shown in Figure 1. The programs are to be installed by the user prior to analysis. This process should be performed individually for each sample.

1. Analyze and visualize DNA sequence data using FastQC ⁴⁶ by typing in to the command line /path-to-file/fastqc, followed by the forward and reverse raw reads raw_read1.fastq raw_read2.fastq.
2. Specify an output folder by typing -o output_fastqc and the file format of the raw read files by -f fastq.
3. View the output file (Figure 2).
   path-to-file/fastqc raw_read1.fastq raw_read2.fastq -o output_directory -f fastq.

9. Quality Control Trimming and Filtering Sequence Data

1. Run the trimming program, Trimmomatic ²⁸ by typing into the command line java -jar /path-to-file/ trimmomatic-0.35.jar.
2. Specify the files are paired end files by typing 'PE'. State that 16 central processing units (CPUs) should be used by the program by typing -threads 16.
3. List the two files to QC check by typing the names of the raw forward and reverse reads. The prefix of the output files is determined by typing -baseout silage.
4. Define the options for the program by typing ILLUMINACLIP:NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 CROP:200 HEADCROP:15 MINLEN:36.
5. Once complete, analyze the trimmed sequences using FastQC as before and compare the output to the raw sequence data to ensure trimming has been performed successfully.
   NOTE: The software tool, Trimmomatic, trimmed reads further by removing leading low quality or N bases (below quality 3), removing trailing low quality or N bases (below quality 3) and scanning each read with a 4-base wide sliding window. The parameters were set for cutting when the average quality per base drops below 20 and then to drop any reads below 36 bases long. Finally, 15 bases were cropped from the head of each read and reads were cropped to keep 200 bases from the start of the read. This final step was performed to overcome some quality issues when sequencing long (> 200 bp) reads. These can be adjusted for specific samples ²⁸.
   java -jar /path-to-file/trimmomatic-0.35.jar PE -threads 16 raw_read1.fastq raw_read2.fastq -baseout silage ILLUMINACLIP:NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 CROP:200 HEADCROP:15 MINLEN:36

10. Metagenome Assembly

1. Merge the unpaired, trimmed reads by typing cat followed by the unpaired reads; silage_read1_unpaired.fastq silage_read2_unpaired.fastq. Write the files to a new file by typing > silage_merged_unpaired.fastq
   cat silage_read1_unpaired.fastq silage_read2_unpaired.fastq > silage_merged_unpaired.fastq
2. To de novo assemble the sequenced DNA, use SPAdes (St. Petersburg genome assembler) ³⁰ by typing /path-to-file/spades.py. Specify that 16 CPUs are to be used by typing -t 16 and that the metagenomic parameter should applied by typing –meta.
3. Identify the trimmed forward reads using -1 silage_read1_paired.fastq and the reverse reads by -2 silage_read2_paired.fastq. The merged unpaired reads are specified by -s silage_merged_unpaired.fastq.
4. Define the output folder by typing -o silage_spades.
   path-to-file/spades.py -t 16 –meta -1 silage_read1_paired.fastq -2 silage_read2_paired.fastq -s silage_merged_unpaired.fastq -o silage_spades

11. Paired-end Read Overlap

1. Merge pairs of DNA sequence reads using FLASH (Fast Length Adjustment of Short Reads) ²⁹ by typing into command line /path-to-file/flash. Specify that 16 CPUs should be used by using -t 16 and the output prefix by typing -o silage.
2. Identify trimmed reads by typing silage_trimmed_R1.fastq silage_trimmed_R2.fastq
   path-to-file/flash -t 16 -o FLASHed silage_read1_paired.fastq silage_read2_paired.fastq

12. Taxonomic Classification

1. Type /path-to-file/kraken and specify the database by typing –db /path-to-file/standard.
2. Define that 16 CPUs should be used by typing –threads 16 and identify an output folder by using –output FLASHed_silage_extendedFrags_kraken.txt. Type the input file name; FLASHed_silage.extendedFrags.fastq
   path-to-file/kraken –db standard –thread 16 –output FLASHed_silage_extendedFrags_kraken.txt FLASHed_silage.extendedFrags.fastq
   NOTE: Classification of the assembled DNA sequence scaffolds using Kraken ⁷ was completed against the most recent, standard Kraken database that contained all available Prokaryote genome sequences.
3. Transfer columns 2 and 3 from the output file and to a new file by typing cut -f2,3 FLASHed_silage_extendedFrags_kraken.txt > FLASHed_silage_extendedFrags_kraken.int
4. Open the output file in web browser.
   cut -f2,3 FLASHed_silage_extendedFrags_kraken.txt > FLASHed_silage_extendedFrags_kraken.int
5. Import the new file into Krona ¹² by typing ktImportTaxonomy. Specify the input file by typing FLASHed_silage_extendedFrags_kraken.int. Identify the output file by typing -o FLASHed_silage_extendedFrags_kraken.out.html.
   path-to-file/ktImportTaxonomy FLASHed_silage_extendedFrags_kraken.int -o FLASHed_silage_extendedFrags_kraken.out.html

13. Functional Annotation

1. Go to the MG-RAST ⁴⁷ website, http://metagenomics.anl.gov/. Register as a new user if required. After logging in, Click on the “Upload” button. Upload the assembled scaffolds from Step 10.
2. Once the files have uploaded, click on "Submit" and follow the instructions and await the completion of analysis.
3. After the analysis is complete, view the link sent via email from MG-RAST, or alternatively, click on "Progress". There is a list of completed jobs. Click on the relevant job id and then on the link to the "download page".
4. On the download page, under the heading "Protein Clustering 90%", click on the protein button to download the predicted protein file, 550.cluster.aa90.faa.
5. To classify the proteins as putatively belonging to a particular CAZy enzyme class, compare the downloaded proteins to the CAZy database ⁴⁸. Download the Carbohydrate-Active enZYmes Database (CAZy) from files are: AA.zip, CE.zip, GH.zip, GT.zip and PL.zip. These files represent the following enzyme classes respectively: Auxiliary Activities (AA), Carbohydrate Esterases (CE), Glycoside Hydrolases (GH), Glycosyl Transferases (GT) and Polysaccharide Lyases (PL).
6. Unzip the database files and annotate the proteins by determining the protein similarity to the CAZy database proteins using the USEARCH UBLAST algorithm ⁴⁹. To use a bash loop (for i in *.txt) to iterate through the 5 database .txt files type "for i in *.txt; do".
7. Run USEARCH by typing /path-to-file/usearch8 with the parameter -ublast in order to use the ublast algorithm. Then type in the name of the protein sequence file downloaded from MG-RAST, "mgmXXXXXX.3.550.cluster.aa90.faa".
8. To indicate the database file to be used type "-db $i" and to specify the E-value threshold at 1e^-5, type "-evalue 1e-5".
9. To terminate the search after the discovery of a target sequence and therefore classifying that protein sequence as belonging to the target enzyme class, e.g. GH, type "-masaccepts 1".
10. To define that 16 CPUs should be used type "-threads 16" and to specify the format of the output file as atab-separated text type "-blast6out". To identify the output file type "$i.ublast". To terminate the bash loop, type "; done"
    for i in *.txt;
    do /path-to-file/usearch8 -ublast ../mgmXXXXXX.3.550.cluster.aa90.faa -db $i -evalue 1e-5 -maxaccepts 1 -threads 16 -blast6out $i.ublast;
    done

14. Visualizing CAZy Annotation

1. To visualize the output from the CAZy annotation as a Venn diagram, generate protein ID lists for each enzyme class using a bash loop. Type "for i in *.ublast; do".
2. To transfer column 1 from the output file and to a new file, type "cat $i | cut -f 1 >$i.list".
3. Terminate the loop and type "; done".
4. Open the .list files in a text editor. Go to the website , select the number of sets as 5 and paste the content of each list file in a separate box. Download the resulting diagram as a .SVG file.
   for i in *.ublast;
   do cat $i | cut -f 1 >$i.list;
   done