Freeze-drying can replace cold-chains for transport and storage of fecal microbiome samples

Sample collection and preservation treatment

Fecal samples were collected from 20 wild-caught Damaraland mole-rats shortly after capture at the Kuruman River Reserve, South Africa, between 1st of April and 15th of October 2019 as part of a long-term study on social behaviours and ecology of this subterranean rodent (e.g., Zöttl et al., 2016; Thorley et al., 2021). The animals were from five randomly selected family groups (four individuals per group). Each sample was split into two replicate tubes with 1–3 fecal pellets per tube and stored in a minus 80 °C freezer at the field site. After all samples had been collected, one replicate per sample was thawed, and subsequently freeze-dried for 48 h at <−40 °C using an ALPHA 1-2 LDplus Freeze Dryer following the manufacturer’s protocol. Freeze-dried samples were then stored and transported at ambient temperature while frozen replicates were transported on dry ice with a commercial courier service to Sweden.

The animal captures and sample collection protocol used in this study was approved by the animal ethics committee of the University of Pretoria (EC050-16) and by Northern Cape Nature Conservation (ODB 1859-2016).

Library preparation and sequencing

The 20 samples were randomly placed on three 96-well plates together with other fecal samples from Damaraland mole-rats (as part of a larger microbiome study). Freeze-dried samples had been stored at room temperature (4 to 10 months) while frozen samples were thawed shortly before extraction. Each plate included four negative controls without a fecal sample and one mock community standard (25 µl ZymoBIOMICS Microbial Community DNA Standard). DNA from samples were extracted inside a UV-hood using the DNeasy PowerSoil Pro Kit (Qiagen) following the manufacturer’s protocol. DNA concentration was quantified using a NanoDrop™1000 spectrophotometer (Thermo Fisher Scientific), and 2.5 ng extracted DNA from each sample was amplified in an initial PCR using the primer 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′) targeting the V3-V4 region of the 16S rRNA gene (Herlemann et al., 2011; Hugerth et al. 2014). The primers included adapter sequences for Illumina n5/n7 index primers used in a second PCR. DNA samples were amplified in 25 µl reactions containing 0.5 µM of each primer and 12.5 µl Phusion High-Fidelity PCR Master Mix (Thermo Fisher Scientific). Thermal PCR cycling conditions used were as follows: 30 s at 98 °C, followed by 20 cycles consisting of 10 s at 98 °C, 30 s at 58 °C and 15 s at 72 °C, and a final 2 min elongation step at 72 °C. The PCR products were purified using AMPure XP magnetic beads (Becker Coulter, USA) and used as template in the second PCR where each sample within a sequencing plate was amplified with a unique combination of Illumina n5/n7 index primers. The 50 µl reaction mix contained 23 µl purified DNA from the first PCR, 0.2 µM index primers and 25 µl Phusion High-Fidelity PCR Master Mix. For the PCR, the following thermal cycling steps were used: 30 s at 98 °C, followed by 12 cycles consisting of 10 s at 98 °C, 30 s at 62 °C and 5 s at 72 °C, and a final 2 min elongation at 72 °C. After the second PCR and another round of purification, DNA concentrations were measured using a Qubit fluorometer (Thermo Fisher Scientific), and equimolar amounts of each sample library from individual sample plates combined into pools with a final concentration of 4 ng/µl. The pools were 300-bp paired end sequenced following standard Illumina sequencing protocols on an Illumina MiSeq platform at the Swedish National Genomics Infrastructure (NGI) at SciLifeLab in Uppsala, Sweden.

Bioinformatics and sequencing filtering

The quality of the reads was checked with FastQC v0.11.8 and MultiQC v1.9 (Andrews, 2010; Ewels et al., 2016). Raw reads from FastQ inputs were processed using the Ampliseq workflow v1.2.0dev (https://nf-co.re/ampliseq/1.2.0, Straub et al., 2020) which uses Cutadapt v.2.8 (Marin, 2011) to identify sequences with primers, remove sequences without primers and delete primers. Sequences passing the primer identification are denoised with the QIIME2 v2019.10.0 (Bolyen et al., 2019) implementation of DADA2 v.1.10.0 (Callahan et al., 2016) and Amplicon sequence variants (ASVs) are created with taxonomy assigned using the SILVA database v.132 (Quast et al., 2013) and QIIME2′s Bayesian classifier. We used the default parameters, besides specifying our own primer sequences and trimming lengths (259 for forward, and 199 for reverse reads) so that sequences were trimmed to equal lengths before the actual denoising, as suggested by the DADA2 protocol. A phylogenetic tree of the ASV sequences was estimated using SEPP on the Greengenes 16S rRNA gene reference tree (McDonald et al., 2012; Mirarab, Nguyen & Warnow, 2012; Janssen et al., 2018).

Quality check and filtering of NGS data

Analyses post Ampliseq were conducted in R version 4.0.1, using mainly functions within the tidyverse, phyloseq and vegan packages (McMurdie & Holmes, 2013; Wickham et al., 2019; Oksanen et al., 2020), and figures were created with ggplot2 (Wickham & Chang, 2016). We combined reads of samples on plates that had been sequenced twice (plate 2 and 3) and removed all contaminant ASVs identified as more prevalent in the negative control samples than true fecal samples by the prevalence method in the decontam package v1.8.0, with a threshold of 0.5 and plate number as batch argument (Davis et al., 2018).

Statistical analysis

ASV richness, Shannon index and Faith’s phylogenetic diversity (PD) were calculated on subsampled ASV counts rarefied to the minimum library size (39323 reads) with phyloseq (McMurdie & Holmes, 2013) which removed 144 ASVs. The effect of sample preservation treatment on library sizes was analysed using a linear mixed model, fitting library size as response variable, including the sample preservation treatment as fixed factor and sample identity and plate number as random factors (Bates et al., 2015). Subsequently, we tested if any of the alpha diversities were associated with sample library size using Pearson correlations. For those alpha diversity measures with a significant association between library size and alpha diversity (ASV richness and PD), we controlled for the effects of library size by calculating the residual alpha diversity from a linear regression between alpha diversity and library size. We assessed the effect of sample preservation treatment on (residual) alpha diversity measures fitting a linear mixed model with (residual) alpha diversity as response variable, the sample preservation treatment as fixed factor and sample identity and plate number as random factors. All mixed models were fitted using the R package lme4 (Bates et al., 2015). Finally, we tested Pearson correlations between alpha diversity measures from sample replicates.

To investigate beta diversity of samples, we performed two principal component analyses (PCA) with the prcomp function on two phylogeny-independent Euclidean distance matrices based on two different transformation methods of raw counts of ASVs: Hellinger transformation (Rao, 1995) using the decostand function in the vegan package (Oksanen et al., 2020) and centered log-ratio (CLR) transformation (Aitchison, 1982) using the clr function in the compositions package. These methods deal with skewed abundance distribution within amplicon microbiome data in different ways. The Hellinger transformation reduces the influence of uncommon ASVs and weighs heavier on the more common ASVs (McMurdie & Holmes, 2013), while the CLR weighs heavier on the rare ASVs (Aitchison, 1982). Additionally, we calculated weighted and unweighted UniFrac distances (Lozupone & Knight, 2005), based on a phylogeny built with SEPP (Janssen et al., 2018). Non-metric multidimensional scaling (NMDS) using the ordinate-function in the phyloseq-package was applied on UniFrac distances (McMurdie & Holmes, 2013). One sample pair which was dominated by the family Enterobacteriaceae (sample 4, see Fig. S1) strongly reduced distances between the other samples in the NMDS (see Fig. S2) and was excluded for further analysis and visualizations.

The amount of variation explained by preservation treatment and sample identity for each beta diversity measure was tested with a Permutational Multivariate Analyses of Variance (PERMANOVA) with the adonis2-function from the vegan package with sample preservation treatment, sample library size and sample identity as fixed factor, by-argument set as “margin”, and testing with and without plate number as the strata argument to control for variation between sequencing plates. Likewise for the NMDS on phylogenetic beta diversity measures, we excluded the outlier sample pair (sample 4, see Fig. S1) in the PERMANOVAs to resemble that of the visualizations of the NMDS. To evaluate multivariate homogeneity of group dispersion of sample preservation treatment and sequencing plates we ran beta disperser tests on each of the beta diversity measures.

Correlations of relative abundances of ASVs within sample pairs were analysed with Pearson correlations across all ASVs and separated by the 8 most common phyla. Lastly, we tested for differences in relative abundances of families between sample pairs, and ran paired Wilcoxon signed-rank tests on families with a mean relative abundance of >1% as relative abundances were non-normally distributed, and p-values were adjusted with Bonferroni correction for multiple testing.

Preservation treatment did not bias library size

We obtained a total of 5,010,719 raw sequence reads from our 40 samples containing 1914 unique ASVs. After removing 146 contaminant ASVs from our samples, we obtained a total of 3,626,584 sequence reads and 1768 unique ASVs. The mean number of reads per sample was 90665 (range freeze-dried = 47510 to 141211, range frozen = 39323 to 183409) and treatment types did not differ significantly in number of reads per sample (p = 0.63).

Preservation treatment did not bias alpha diversity

Although there was substantial between-sample variation, the preservation treatment did not significantly bias measures of alpha diversity. For both sample treatments, we found large variation in ASV richness (freeze-died: mean = 304.55, range = 202–405; frozen: mean = 303.30, range = 202–402), Shannon index (freeze-dried: mean = 3.96, range = 3.20–4.55; frozen: mean = 3.97, range = 2.42–4.54) and Faith’s PD (freeze-dried: mean = 38.24, range = 27.41–47.63; frozen: mean = 37.13, range = 24.86–47.44) among samples. Some of that variation was explained by library size (ASV richness: R = 0.68, p = <0.001; Shannon index: R = 0.12, p = 0.48; Faith’s PD: R = 0.6, p < 0.001). However, after controlling for the effect of library size, the preservation treatment of the sample replicates did not significantly bias ASV richness, Shannon index and Faith’s PD (Fig. 1A residual ASV richness LMM (Estimate +- Std. Error): Intercept = 1.73 +- 8.09, estimate Treatment = 5.05 ± 7.542, p = 0.09; Fig. 1B Shannon index LMM: Intercept = 3.94 ± 0.10, estimate Treatment = 0.004 ± 0.07, p = 0.96; Fig. 1C residual Faith’s PD LMM: Intercept = 0.69 ± 0.98, estimate Treatment = −1.46 ± 0.87, p = 0.10), and the correlations between replicate samples were positive (Fig. 1D residual ASV richness: R = 0.64, p =0.003; Fig. 1E Shannon index R = 0.7, p < 0.001; Fig. 1F residual Faith’s PD R = 0.62, p < 0.001).

Preservation treatment has minor effect on beta diversity

Compositional differences between samples of the preservation treatments were minor. PCAs of Euclidean distances on both Hellinger and CLR transformed counts of ASVs revealed no clear clustering of sample treatment (Fig. 2). Similarly, NMDS on phylogenetic distance measures did not reveal clustering by sample treatment (Fig. 3).

Consistent with the interpretation of the PCAs and the NMDS, a PERMANOVA analysis revealed that sample treatment explained only a small proportion of the variation among the samples (1.5–2.2%) whereas the sample identity was identified as the main source of variation (73.2–86.6%, Table 1). Library size was significant for Euclidean distances on CLR-transformed data and unweighted UniFrac, but only explained 1.8−3.1% (Table 1). The PERMANOVA analysis further confirmed that plate number did not significantly explain any of the variation in microbiome composition of any of the four distance measures (Table 1). Finally, we did not detect any significant differences between sample treatments with Beta dispersion tests (p = 0.08−0.94, Table 1).

Preservation treatment has a minor effect on compositional differences

Out of the 1768 ASVs, 98.5%, 91.9% and 73.7% were assigned to a phylum, family or genus, respectively. The dominating phyla were Bacteroidetes (mean relative abundance 72.9%; range 31.4–88.9%) and Firmicutes (mean relative abundance 16.4%; range 4.9–47.8%). However, within one sample, one ASV of the family Enterobacteriaceae dominated the community composition in both sample treatments with as much as 43.7% and 58.6% of the relative abundance of the freeze-dried and frozen replicates, respectively (see sample 4, Fig. S1). The ASV was not unique to the sample replicates, but prevalent in five other samples but in much lower abundance. Among the families with >1% mean relative abundance, only three families within the phylum Firmicutes, Christensenellaceae, Ruminococcaceae and Lachnospiraceae, was significantly different in relative abundances between treatments and was underrepresented within freeze-dried samples (Fig. 4, adjusted p < 0.001, Table S1).

About half (54%) of the ASVs were shared between the two sample treatments (Fig. S3) and these reads summed to a total of 99.5% of the reads of the full dataset. The two preservation treatment types had very similar numbers of unique ASV, and the majority of these were also unique to a single sample (Fig. S3). Overall, we found a strong correlation between relative abundances of ASVs between the two types of treatments (R = 0.87, p < 0.001). When analysing the phyla separately, the more common phyla still showed strong correlations (Fig. 5). The phyla Lentisphaerae and Spirochaetes showed weaker trends between treatment types (Fig. 5), and Lentisphaerae was the only phylum that had overrepresentation of unique ASVs in one of the sample pairs (the freeze-dried treatment, Fig. 5).