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

Background The transport and storage of samples in temperatures of minus 80 °C is commonly considered as the gold standard for microbiome studies. However, studies conducting sample collection at remote sites without a reliable cold-chain would benefit from a sample preservation method that allows transport and storage at ambient temperature. Methods In this study we compare alpha diversity and 16S microbiome composition of 20 fecal sample replicates from Damaraland mole-rats (Fukomys damarensis) preserved in a minus 80 °C freezer and transported on dry ice to freeze-dried samples that were stored and transported in ambient temperature until DNA extraction. Results We found strong correlations between relative abundances of Amplicon Sequence Variants (ASVs) between preservation treatments of the sample, no differences in alpha diversity measures between the two preservation treatments and minor effects of the preservation treatment on beta diversity measures. Our results show that freeze-drying samples can be a useful method for cost-effective transportation and storage of microbiome samples that yields quantitatively almost indistinguishable results in 16S microbiome analyses as those stored in minus 80 °C.


INTRODUCTION
The transport and storage of microbiome samples on dry ice and in minus 80 • C freezers (i.e., super freezers) is commonly considered the gold standard (Vandeputte et al., 2017). However, for research carried out in remote locations, at field stations or when the analytical laboratory work is carried out at laboratories far away from the collection site, it can be challenging to provide super freezers and maintain uninterrupted cold chains (Kim et al., 2017;Vandeputte et al., 2017). Super freezers require a reliable source of electricity and are expensive to purchase and maintain. Ordinary freezers may be a cost efficient alternative and increasing evidence suggests that DNA storage of microbiome samples in those freezers yield similar 16S rRNA gene sequencing data profiles as samples stored in super freezers (Song et al., 2016;Gavriliuc et al., 2021). However, because the microbiota composition in frozen samples is sensitive to thawing (Cardona et al., 2012;Song et al., 2016), freezing requires an uninterrupted cold chain from sampling to library preparation and samples are therefore often transported on dry ice or in liquid nitrogen from the field sites to the laboratory. Transportation of frozen samples hence introduces the risk of thawing and incurs additional costs of cold chain logistics. This is especially problematic when the sampling location is in a remote area or in another country than the laboratory facility and commercial couriers are needed.
When freezing is not an available option, the most commonly used and tested alternative storage of samples include ethanol, RNAlater or FTA-cards (Hale et al., 2015;Song et al., 2016;Vogtmann et al., 2017;Wang et al., 2018). Although not requiring immediate freezing, these sample treatments are commonly frozen after a couple of hours or days at ambient temperatures for long-term storage prior to library preparation. Out of these alternative methods, ethanol is the cheapest and has been shown to be the most effective for preserving the microbiome community structure when freezers are not available (Hale et al., 2015). However, the use of ethanol comes with drawbacks, for example the additional step of ethanol removal prior to DNA extraction, as well as risk of evaporation. Furthermore, ethanol is both volatile and flammable, which increases shipping costs and restricts export transport options (Song et al., 2016;Vandeputte et al., 2017). Consequently, methods where samples can be transported at room temperature and without hazardous liquids could be of great value to research projects working at remote locations.
In this study, we evaluate freeze-drying as an alternative method for the preservation of fecal samples prior to transport, storage and 16S sequencing. Although freeze-drying has rarely been used in microbiome studies (Gavriliuc et al., 2021), it offers considerable advantages. No chemicals are added and the samples can be stored and transported at ambient temperature, eliminating the need for cold chain logistics and hazardous material handling. Previous studies have offered some support for this preservation technique by showing an overlap in 16S data from freeze-dried human stool samples with data from frozen samples (Kia et al., 2016) and little influence of freeze-drying on homogenized neonatal fecal samples (Shen et al., 2021). However experimental work evaluating the effects of freeze-drying on field samples and samples from non-human animals is still lacking. To evaluate whether freeze-drying, subsequent long-distance transport, and storage at ambient conditions affects the results of 16S sequencing, we conducted an experiment where we investigated the gut microbiome diversity and composition of 20 replicate fecal samples from Damaraland mole-rats (Fukomys damarensis) by comparing freeze-drying with super freezers.

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 TM 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.

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 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).

DISCUSSION
Fecal samples from wild animal populations can yield important insights into the role of the gut microbiome for the fitness, health and survival of their hosts (Suzuki, 2017). However, the transport of samples in stable cold chains from remote field collection sites to laboratories can be challenging and ways of preserving fecal samples for microbiome analyses without freezing or chemical preservatives would facilitate the development of this research field. In this study we show that frozen and freeze-dried samples do not differ in alpha diversity measures and the variation induced by the treatment amounted to less than 2.2% of the variation across the samples. Alpha diversity measures did not differ between the two preservation methods, similar to other studies comparing alternative sample methods to deep-frozen microbiome samples (Song et al., 2016;Gavriliuc et al., 2021). Because freeze-drying did not bias any of the three alpha diversity measures we investigated, contamination or shifts in the sample community composition during the additional freeze-drying procedure are likely rare. Compared to freeze-drying, some other preservation methods that successfully eliminate the necessity of freezing sometimes come with the drawback of introducing the risk of altering sample alpha diversity through either contamination by additional bacterial exposure or change how easily different taxa are extracted. For example, Song et al. (2016) found that FTA-cards had increased Shannon diversity compared to other sample methods and fresh samples. Furthermore, the library sizes in our study did not differ systematically between the treatment types ( Fig. S4) which suggests that rare ASVs are sequenced with similar probabilities (Wu et al., 2010), and both methods perform equally using the same DNA-extraction and library preparation protocol. Our data suggests that none of these possible biases were detectable in the comparison between frozen and freeze-dried samples. For all beta diversity measures, sample preservation treatment explained a minor proportion of the variation (1.5-2.1%) within the microbiome community composition whereas sample identity explained a much larger part of the variation (73.2 -86.6%). Sample identity explained more of the variation in beta diversity calculated on Hellingertransformed than CLR-transformed data. Hellinger transformation puts weight on abundant ASVs (Legendre & Legendre, 2012) and might hence be less sensitive to random sampling of rare ASVs which might explain the higher impact of sample identity for this transformation. Alternatively, the CLR-transformed data was better at picking up differences between pellets from the same animal as each replicate tube from an animal contained one to two pellets and replicate samples did not originate from a homogenized mix of pellets. While the variation explained by treatment was still small compared to the sample identity, it suggests that taxonomic groups are most often prevalent in both sample pairs but that some taxonomic groups vary slightly in abundance between treatments (see further discussion below). Sample library size explained a similar proportion of the variation as sample treatment, and because sample treatment did not bias library sizes, this suggest that the amount of variation caused by sample treatment was similar to random variation in library sizes between samples. The relative abundances of ASVs were overall strongly correlated between the two treatments, suggesting that freeze-died samples reliably reflect the prevalence and abundance of bacterial populations in deep-frozen samples. Notably, one sample pair was dominated by the same single ASV in both sample treatments, showing that both treatments succeeded at capturing ''extremes''. Given that this ASV had unusually high abundance within both treatments of the sample and was prevalent in other samples too, it was likely a truly dominating ASV within the gut microbiome of the individual and not a contamination. Furthermore, both treatments had similar proportions of sample-pair unique ASV abundances of the more common phyla and these were in general lower in relative abundance than ASVs shared within sample pairs, demonstrating that both treatments do equally well with sequencing less dominant ASVs. The phylum with the lowest correlation was Lentisphaerae, in which ASVs unique to one of the samples were overrepresented within freeze-dried samples (Fig. 4)  ASVs were in general strong within the dominating phyla, we found that among families with >1% mean relative abundance, three families within the second most common phylum Firmicutes differed significantly between treatments. It is likely that differences in abundance of these taxa were responsible for the small effect of treatment on variation in in microbiome community structure. Given that this was still small compared to the effect of sample identity in all four beta diversity measures, our conclusion is that freeze-drying has only minor influences on the composition and reflect the overall diversity of the fecal microbiome well. It is however important to note that fecal microbiome samples are widely known to slightly change in composition with time and sample preservation method and the importance of not mixing different preservation methods to avoid batch effects is widely acknowledged (Song et al., 2016). Because our results suggests that there may be a minor but significant bias in the community composition of freeze dried samples of a magnitude of 1-2% percent of the variation, we suggest that researcher considering mixing both frozen and freeze dried samples should consider whether the benefits of mixing outweigh the costs of introducing this small bias.

CONCLUSION
Our results show that freeze-drying can be a suitable method when microbiome samples have to be exported over long-distances and stable cold chains are not available. Alpha diversity measures were not biased between frozen and freeze-dried samples and the sample preservation treatment had only a minor influence on the community composition whereas sample identity explained a large proportion of the variance. Together, our findings suggest that already collected samples stored in freezers can be thawed and freeze-dried to be transported and stored at ambient temperatures prior to 16S microbiome analyses.