The effect of reverse transcription enzymes and conditions on high throughput amplicon sequencing of the 16S rRNA

Study site and sample collection scheme

Soil samples were collected at the central Negev Desert highlands, Israel (Zin Plateau, 30°86′N, 34°80′E) at an established ecological research site. The mean annual precipitation at the sampling site is 90 mm and the mean annual temperature is 30° (LTER data). Samples were collected under the canopy of perennial shrub Hammada scoparia in October 2015 at the end of the dry season as previously described (Baubin et al., 2019). Briefly, sampling was conducted in seven random blocks, each providing two technical replicates resulting in 14 samples. Samples were collected from the top five cm of the soil, following the removal of crust and debris. The soil samples were processed within 24 h of collection. Samples were homogenized using two mm sieve and the duplicates from each block were composited. This resulted in seven final replicates.

RNA preparation

Total RNA was extracted from each of the seven samples using a phenol-chlorophorm extraction previously described by Angel (2012). The reaction buffer pH was adjusted to 5. The total RNA was subsequently purified with the MagListo™ Total RNA Extraction Kit (Bioneer, Daejeon, Republic of Korea). Contaminant DNA was removed using a DNAse I from the MasterPure RNA Purification Kit (Epicenter, Madison, WI, USA) with two successive treatments of 30 min according to manufacturer’s instructions. The reaction mixture was purified using the MagListo™ Total RNA Extraction Kit (Bioneer, Daejeon, Republic of Korea). The absence of contaminant DNA was verified using total bacterial primers 341F (5′-CCTACGGGAGGCAGCAG-3′) and 515R (5′-TTACCGCGGCTGCTGGCAC-3′) (Klindworth et al., 2013) and DreamTaq DNA polymerase (Thermo Scientific, Waltham, MA, USA) under the following conditions: 95 °C for 5 min, followed by 26 cycles of 95 °C for 15 s, 60 °C for 30 s and 72 °C for 30 s for extension followed by 72 °C for 5 min for final extension. If amplification was detected the sample was discarded, re-extracted, purified, and tested. Only DNA-free samples were used in this study.

Reverse transcription reaction conditions

RT kits used in this experiment were chosen to represent the most commonly used enzymes in the field (Fig. 1). All seven samples were reverse transcribed by the same kit to reduce batch effects. Each RT kit used in this study originated from a different source: (I) ImProm-II Reverse Transcription System enzyme (Promega, Madison, WI, USA) originates from AMV RT, (II) SuperScript IV Reverse Transcriptase Kit enzyme (ThermoFisher Scientific, Waltham, MA, USA) originates from MMLV RT and (III) TGIRT originates from the mobile group II introns reverse transcriptase (Mohr et al., 2013) TGIRT™-III Enzyme (InGex, St. Louis, MO, USA). Each reaction consisted of 50 ng of total RNA template, measured by Quanti-iT™ RNA Assay Kit (ThermoFisher, Waltham, MA, USA), and random hexamer primers (0.5 μg/reaction). Template and primer mix were heated to 70 °C (ImProm-II) or 65 °C (SuperScript IV). Each reaction was subsequently cooled to 4 °C for 5 min and incubated at 42 °C (ImProm-II), 55 °C (ImProm-II and Superscript) or 57 °C (TGIRT) for 60 min (ImProm-II), 120 min (TGIRT), or 10 min (SuperScript IV). All reactions were terminated and DNA was removed by alkaline lysis using 2 μL of 1M NaOH, incubating for 12 min at 70 °C. After which the reaction was neutralized using 4 μL of 0.5M acetic acid (Table 2).

Illumina sequence preparation

The V3 and V4 regions of the resulting cDNA were amplified using 341F (5′-CCTACGGGAGGCAGCAG-3′) and 806R (5′-GGTCTGGACTACHVGGGTWTCTAAT-3′) (Klindworth et al., 2013) primers. Each reaction was performed in triplicate and consisted of 1 mM bovine serum albumin (Takara, Kusatsu, Japan), 2.5 μL of 10× standard buffer, five μM primers, 0.8 mM dNTPs, 0.4 μL DreamTaq DNA polymerase (Thermo Scientific, Waltham, MA, USA), and 4 μL of template cDNA. The reaction mixtures were subsequently amplified using the following PCR conditions: 95 °C for 30 s, 27 cycles of 95 °C for 15 s, 50 °C for 30 s, 68 °C for 30 s, and 68 °C for 5 min. Resulting amplicon presence was verified using 1.5% agarose gel electrophoresis. Resulting technical triplicates were combined, and the sequencing libraries were constructed using the TruSeq^® DNA Sample Preparation Kit (Illumina, San Diego, CA, USA) following the manufacturer’s recommendations. The amplicon libraries were sequenced (250 × 2 base pairs, pair-end) on the Illumina MiSeq System platform at the Research Resources Centre at the University of Illinois.

Sequence analysis

Resulting paired end sequences were merged using the CASPER program (Kwon, Lee & Yoon, 2014), and the resulting merged reads were clustered using the UPARSE pipeline according to the recommended settings (Edgar, 2013). The resulting operational taxonomic unit (OTU) representative sequences were taxonomically assigned with the SINA incremental aligner using a lowest common ancestor algorithm (Pruesse et al., 2007) and the SILVA database version 132 (Quast et al., 2013). All sequences retrieved in this study were uploaded to European Nucleotide Archive (https://www.ebi.ac.uk/ena) submission number PRJEB32237.

Class enrichment plot data preparation

To explore whether G-C content contributed to differences in relative abundances of different taxa, we ran the following analysis: for each reaction condition we calculated pairwise comparisons at the class level: we normalized the proportional enrichment in each respective reaction conditions following Eq. (1) (Fig. 2; Fig. S2), where the A and B represent a class at the different conditions.

The error bars represent a standard deviation, which have been calculated as a standard deviation of each category and normalized according to the Eq. (2). The δA and δB represent the standard deviation of class A and B. Details of deriving this equations are specified in Equation_S1 and Code_S3.ipynb.

Statistical analysis

All data analysis was performed in R v3.4.3 (R Core Team, 2018). The dataset was sub-sampled (rarified) to an even depth of 9,000 sequences per sample using python numpy package v1.15.4 (Van Der Walt, Colbert & Varoquaux, 2011). The subsampling removed five samples from the dataset. Additionally, three more samples were removed as outliers (Code_S1.ipynb). In order to equalize the number of replicates, two random samples were removed from the TGIRT dataset. This reduced the number of replicates to four samples per experimental category. The sample diversity was analyzed using the vegan package v2.5-2 (Oksanen et al., 2018). The data was visualized using the R package ggplot2 package v2.2.1 (Wickham, 2016) and python package matplotlib (Hunter, 2007).

Sample preparation and diversity analysis

After removing obvious outliers, four samples from each condition were analyzed. The summary of the analysis and all the code used to produce each figure is included in the Files S1 and S2. The changes in bacterial diversity among the tested conditions were expressed using species count, Pielou’s evenness index (Pielou, 1966) and Shannon diversity index (Shannon & Weaver, 1949). Despite observed trends in the diversity indices, no statistically significant differences were detected (Fig. S1; Code_S2.ipynb).

Relative abundance plot

A relative abundance of major taxonomic classes is depicted in the Fig. 1. Each column is an average of four biological replicates. OTUs that were not taxonomically assigned at this level are summarized as “Unclassified” (≈2%). Various patterns (detailled below) were detected among the experimental conditions: some patterns could be attributed to differences in reaction temperature (which ranged from 42 to 55 °C and 57 °C). Other patterns could be linked to enzyme type. RT reactions with SupeScript IV and TGIRT RT enzymes yielded no significant differences in class relative abundances. However, transcription with ImProm-II RT at a similar temperature (≈55 °C) yielded different abundances: specifically, the abundances of Alphaproteobacteria, Bacteroidia, Deltaproteobacteria, Oxyphotobacteria, Rubrobacteria, and Verrucomicrobidae decreased. However, Chloflexia, Gammaproteobacteria and Thermoleophilia abundances increased when their ribosomes were transcribed with ImProm-II RT (Fig. 1; Table S1). When transcription occured at lower temperature (42 °C), relative abundances of Bacilli, Deltaproteobacteria, and Oxyphotobacteria were enriched, while Actinobacteria, Chloroflexia, and Acidimicrobia were depleted under the same conditions (Fig. 1; Table S1).

Class enrichment plot

The Fig. 2 also depicts the weighted average of the G-C content in each class. The proportional comparison is interpreted as follows: A value of zero in the proportional comparison represents the taxonomic class count that is exactly equal between the two compared groups (Fig. 2). A value of 1 or −1 is assigned when a given taxonomic class is only present in one category and absent in another, respectively. In general, there was a tendency toward lower G-C content lower temperature of ImProm-II (Fig. 2A; Figs. S2A and S2B). No statistically significant differences were detected between the profiles resulting from RT of SuprScript-IV or TGIRT (Fig. S2D). The taxa Gemmatimonadetes, Fibrobacteria, and Thermoanaerobaculia were relatively insensitive to the RT conditions. However the majority of the classes were enriched in some conditions. (I) the rate of RT reaction was only sensitive to temperature for the classes Alphaproteobacteria, Gemmatimonadetes, Fibrobacteria, Thermoanaerobaculia, TK10, and Blastocatellia (Fig. 2A). These groups tend to have extreme GC content (both high and low). (II) Gammaproteobacteria, Planctomycetacia, and Phycisphaerae are relatively insensitive to the reaction temperature (Fig. 2A), but their abundances vary with different RT enzymes (Fig. 2B; Fig. S2).

Class enrichment statistics

We calculated a linear regression, where the response variable was the relative proportion of each class between two tested categories, and the explanatory variable was the G-C content (the assumptions tests and plots can be found in the File S2). The linear regression assumptions were tested: in case of two condition pairs (ImProm-II at 42 °C & SuperScript-IV as well as ImProm-II at 55 °C & SuperScript-IV), the assumption conditions were not met (File S2). Since we cannot confidently discard the null hypothesis in these cases, we do not consider these two tests significant. Therefore, we are considering only the ImProm-II 42 °C & ImProm-II 55 °C as well as ImProm-II 42 °C & TGIRT as a significant outcome. Differences in the remaining cases cannot be explained by the weighted G-C content alone (Table 3).

We used a GC content as an explanatory variable of a the class enrichment. The rows marked with a • do not fulfill all test assumptions (see Code_S2.ipynb). The statistical significance: *, **, *** at P < 0.05, 0.01 and 0.001 respectively.