Optimized DNA extraction and purification method for characterization of bacterial and fungal communities in lung tissue samples

Human lungs harbor a scarce microbial community, requiring to develop methods to enhance the recovery of nucleic acids from bacteria and fungi, leading to a more efficient analysis of the lung tissue microbiota. Here we describe five extraction protocols including pre-treatment, bead-beating and/or Phenol:Chloroform:Isoamyl alcohol steps, applied to lung tissue samples from autopsied individuals. The resulting total DNA yield and quality, bacterial and fungal DNA amount and the microbial community structure were analyzed by qPCR and Illumina sequencing of bacterial 16S rRNA and fungal ITS genes. Bioinformatic modeling revealed that a large part of microbiome from lung tissue is composed of microbial contaminants, although our controls clustered separately from biological samples. After removal of contaminant sequences, the effects of extraction protocols on the microbiota were assessed. The major differences among samples could be attributed to inter-individual variations rather than DNA extraction protocols. However, inclusion of the bead-beater and Phenol:Chloroform:Isoamyl alcohol steps resulted in changes in the relative abundance of some bacterial/fungal taxa. Furthermore, inclusion of a pre-treatment step increased microbial DNA concentration but not diversity and it may contribute to eliminate DNA fragments from dead microorganisms in lung tissue samples, making the microbial profile closer to the actual one.


Results
Quality and quantity of genomic DNA from five extraction protocols. Lung tissue specimens from seven autopsied individuals were processed for DNA extraction following the five methods described (Table 1) and the concentration and quality measured (Fig. 1A). Protocol 2 including the bead-beating and 3 including the Phenol:Chloroform:Isoamyl alcohol step, had slightly increased DNA yields compared to protocol 1; however these findings were not significant. Combining both additional steps, protocol 4 significantly increased DNA concentrations (1,104.8 ± 522.7 ng/µl) compared to protocols 1 (525.4 ± 460.2 ng/µl, p = 0.003) and 2 (557.4 ± 245.2 ng/µl, p = 0.036). Protocol 5 (686.2 ± 172.3 ng/µl) also showed a relative, but non-significant, increase in DNA compared to protocol 1.
Evaluation of human, bacterial, and fungal DNA in the five extraction protocols. The evaluation of the content of human, bacterial, and fungal genomes through qPCR amplification of human β-actin, 16S rRNA and 18S rRNA genes respectively showed that the DNA extraction method did not affect the Ct obtained for the β-actin gene, suggesting that the level of human DNA was equivalent for all extraction protocols ( Fig. 2A). However, the Ct of the 16S rRNA gene (mean ± SD) was significantly decreased with protocol 5 (22.9 ± 2.2) compared to all others (26.3 ± 2.2 for protocol 1; 25.3 ± 2.5 for protocol 2; 26.3 ± 1.9 for protocol 3; and 26.1 ± 2.4 for protocol 4) (Fig. 2B).
We detected low quantities of fungi in lung tissue, as the Ct values of 18S RNA gene ranged from 37.3 to 40 (Fig. 2C). Although there were no significant differences, the Ct tended to be decreased in samples processed with protocol 5 (38.9 ± 1.0) compared to others (39.6 ± 0.7 for protocol 1; 39.2 ± 0.9 for protocol 2; 40.0 ± 0.0 for protocol 3; and 39.7 ± 0.5 for protocol 4).
To evaluate the contamination level of our samples, we identified reads mapping 100% identity with those of negative controls. For bacteria, the number of reads removed was 337,935 (61.7% of the reads before clustering), reaching 9655 ± 13,707 (mean ± SD) reads per sample. For fungi, figures were lower, with 168,683 (8.7% of the reads) and 4820 ± 7295 reads per sample. The negative controls contained 1580 reads (6.9%) showing 100% identity with at least one bacterial sample, whereas in fungi this figure was lower, 258 reads (0.4%). Those reads were removed from the samples. Despite the relatively high number of reads in the negative controls, most of www.nature.com/scientificreports/ them were actually unique and restricted, as a consequence the impact of the removal of reads due to potential contamination on the samples was relatively low. To ensure that the microorganisms detected in our samples were not completely a result of contamination, we compared the microbial community of our samples with that of negative controls by performing a CCA based on Bray-Curtis distances for bacterial and fungal communities ( Supplementary Fig. S1A,B). We observed that the majority of samples and negative controls clustered separately for both bacteria and fungi, which was confirmed by the Adonis test (p = 0.003 for bacteria and p = 0.017 for fungi).
Lung microbial community diversity. Sequencing of the V3-V4 regions of the 16S rRNA gene carried out on the 35 samples from the 7 individuals using 5 protocols yielded 1,005,778 raw reads, averaging 28,737 reads per sample, and ranging from 4717 to 137,755 reads. Similarly, the sequencing of the ITS regions in these  After removal of sequences that clustered with negative control reads at 100% identity the composition and abundance of the remaining 209,533 bacterial and 1,776,971 fungal "contaminant-free" reads are summarized in Fig. 4 (see Supplementary Fig. S2 for individual data).
Impact of DNA extraction protocol on the bacterial community. We investigated the effect of the extraction protocols on bacterial community detection in lung tissue samples. In the Fig. 5, CCA analysis shows that DNA extraction protocols did not drive significant differences. Rather, samples clustered by individual rather than by extraction protocol (Fig. 5). Furthermore, there were no significant differences in the Shannon index between the different DNA extraction protocols ( Supplementary Fig. S3).

Scientific Reports
| (2020) 10:17377 | https://doi.org/10.1038/s41598-020-74137-2 www.nature.com/scientificreports/  www.nature.com/scientificreports/ Impact of DNA extraction protocol on the fungal community. We investigated whether the DNA extraction protocols affected the composition of fungal community. CCA analysis showed that samples tend to cluster by extraction protocol, however this grouping was not statistically significant, suggesting that the extraction method affects the fungal community (Fig. 5). In addition, we observed that samples cluster according to the subjects (Fig. 5). Furthermore, there were no significant differences in the Shannon index between the different DNA extraction protocols ( Supplementary Fig. S3). Compositional comparison revealed that the fungal families identified ranged from 9 in protocol 5 to 24 in protocol 2, with figures for the other protocols of 13 (protocol 4), 15 (protocol 3), and 17 (protocol 1). We noted that the bead-beating step improved fungal family retrieval (protocol 2), although this increase was not found for protocol 4, which also the Phenol:Chloroform:Isoamyl alcohol steps.
To determine the effect of bead beating, we compared changes in the relative abundance of taxa between protocols 2 and 1. The following families showed log 2 fold-change increases in protocol 2 ( Fig. 7): Malasseziaceae  Table S3 for individual data). In contrast, the Phenol:Chloroform:Isoamyl alcohol step appeared to induce fewer changes (protocol 3) compared to protocol 1, as similar number of families were recovered (15 and 17 detected, respectively). However, log2 fold-change increases were observed in the phylum Basidiomycota, including the family Malasseziaceae (+ 5.8) and in the phylum Ascomycota, including Hypocreales_fam_Incertae_sedis (+ 5.2) and Cladosporiaceae (+ 0.9).

Discussion
Here, using near-death autopsied lungs from people dying in accidents, we provide the first description of the bacterial and fungal community inhabiting the lungs in healthy human subjects. Furthermore, we evaluated the effect of five different DNA extraction protocols on the quantity and quality of nucleic acids, as well as on the bacterial and fungal community composition. We show that additional steps, including bead-beating, Phenol:Chloroform:Isoamyl alcohol, and pre-treatment affect the characterization of the microbial community in lung tissue samples. Over the last decade, the development of culture-independent techniques for microbiological analysis has uncovered the existence of a microbial community in the lung, which was previously considered to be sterile. The lower respiratory tract of healthy subjects was mainly colonized by anaerobic bacteria, such as the phylum Bacteroidetes including the family Prevotellaceae and the phylum Firmicutes including the families Veillonellaceae and Streptococcaceae, as well as aerobic bacteria, such as the phylum Proteobacteria including the families Pseudomonaceae, Neisseriaceae, and Pasteurellaceae 14,15 . However, these data are based on studies whose samples were contaminated by microorganisms of the upper respiratory tract (i.e. sputum, bronchial aspirate, lung biopsy and bronchoalveolar lavage), as during their collection they transited through the oropharynx; thus, they do not represent the commensal communities inhabiting the lung 16 . Except for sputum, the collection of these samples is ethically difficult to perform in healthy subjects, as they require invasive procedures, making the characterization of the healthy lung microbiome extremely difficult.
In this study, we choose to use lung tissue samples obtained near death from subjects for characterizing the lung microbiome, as the initial postmortem microbiome would be a reflection of the host microbiome preceding death. The microbiome composition of distinct body habitats during the first 48 h postmortem showed still differences and remained stable in the first two days after death compared to those analyzed at a later time 17 . These data suggest that anatomic locations were not still contaminated by postmortem transmigration. In addition, the lung, like the brain, would be less sensitive to postmortem biological changes, compared to other organs. In effect it was shown that the lung was less sensitive to postmortem mRNA degradation 18,19 . As the microcirculation continues working for some time after death, it is hypothesized that tissues with abundant small blood vessels and microcirculation such as brain and lung would be less sensitive to postmortem changes. Therefore, despite the fluidic process of decomposition, the body can maintain its microbial communities with strong niche differentiation among anatomic locations. Together these data highlight the valuable biological resource of the human postmortem tissue (within 48 h of death) for characterizing the lung microbiome.
Our study shows that lung tissue samples from healthy subjects were composed mainly of the phylum Firmicutes, largely represented by the families Anaerobacillaceae, Streptococcaceae, Clostridiaceae and Veillonellaceae. We also detected the phyla Bacteroidetes, mainly represented by the family Prevotellaceae, and Proteobacteria, mainly the families Neisseriaceae and Pasteurellaceae, in lower concentrations. We observed great diversity in the microbiome inhabiting the lung. Although many fungal sequences remain unidentified, due to a lack of reference genomes 20 , the families we identified came primarily from the phyla Ascomycota and Basidiomycota, as well as Saccharomycetae, Saccharomycetales family Incertae sedis, and Aspergillaceae. Interestingly, species of some identified families such as Candida, Malassezia and Saccharomyces, have been previously reported to cause lung infection, suggesting that they may be present in the lungs as commensals 21 .
The study of the lung microbiota remains a technical challenge for researchers due to the low microbial biomass inhabiting the human lung, estimated at approximately 2.2 × 10 3 bacterial genomes per cm 214,22 . Consequently, these samples are more sensitive to contamination during the DNA extraction process, affecting sequencing data and interpretation of the microbiota. This is even more of a concern for fungi, present in lower concentrations in our samples compared to bacteria. This problem was well illustrated in a study by Lauder and collaborators, who characterized the microbiota of the placenta, which is also composed of a low microbial biomass 23 . The authors were unable to detect significant differences between the placental samples and contamination controls. In contrast, our study shows a significant separation between lung tissue samples and environmental contamination.
However, we identified potential contaminants in our samples with bacterial contamination more common than fungal one (67.7% vs 8.7% of the reads). Contaminants were identified through a computational approach, based on 100% identity sequences present in both negative controls and samples. This strategy appears more successful in removing sequences, which is vital as the lung microbiota can be contaminated by the bacterial content of inhaled air. It is likely that some bacteria of the lung microbiota are the same found in contaminant indoor air. Nevertheless, it is probable that some taxa from the lung microbiota were removed through our strategy. As the length of sequencing reads was only about 460 bp, distinct microbial strains present in the lung microbiota and in the contaminants could have 100% identity sequences. Thus, we probably overestimated the level of contamination in our samples.
In microbiome studies, it is well known that DNA extraction protocols are directly related to the quality of sequencing and taxonomic identification of microorganisms. The need for high concentration and quality of extracted DNA is of special concern in low biomass samples, such as lung tissue samples. In addition to the low microbial concentration, lung tissue contains mostly human DNA, making it harder to amplify and sequence microbial DNA. In this context, choosing an adequate DNA extraction protocol ensures an efficient recovery of all microorganisms present in the samples and the quality of the extracted DNA.
In our study, the addition of a bead-beating step in protocol 2 did not significantly increase the recovery of extracted DNA compared with protocol 1. However, the mechanical action of the beads on the microbial wall improves microbial disruption favoring the detection of Gram-positive bacteria from the phyla Firmicutes (e.g. families Bacillaceae and Clostridiaceae) and Actinobacteria (e.g. Actinomycetaceae and Atopobiaceae) that have many layers of peptidoglycan in their thick cell wall, which is not easily disrupted. Moreover, although the Scientific Reports | (2020) 10:17377 | https://doi.org/10.1038/s41598-020-74137-2 www.nature.com/scientificreports/ Proteobacteria were reduced, we also retrieved increased abundances in some Gram-negative bacteria from the phyla Bacteroidetes (e.g. families Weeksellaceae and Porphyromonadaceae) and Proteobacteria (e.g. Halomonadaceae, Rhizobiaceae and Sphingomonadaceae). Previous DNA extraction methods based on enzymatic treatment without physical disruption showed reduced recovery of Gram-positive bacteria and relatively elevated levels of Gram negatives 24,25 . It is noteworthy that fungi often have a cell wall that is harder to lyse than bacterial cell walls, and DNA extraction kits are generally not optimized for fungal DNA extractions. We observed a higher abundance of the families Malasseziaceae and Aspergillaceae in samples processed with bead beating. It has been previously shown that extraction protocols that employ bead beating increase DNA yields of genera from the Aspergillaceae 26 and the Malasseziaceae families 27 . In addition to these two families, our study revealed that the bead-beating step also increased the recovery of other fungal taxa, such as the families Cladosporiaceae and Dipodascaceae.
On the other hand, the Phenol:Chloroform:Isoamyl alcohol step alone in protocol 3 leads to modest changes in the bacterial and fungal communities. However, we observed that this step ensures a better quality of the extracted DNA (ratio 260/280-absorbance) compared to protocol 1. While the phenol can aid in disrupting the cell wall by denaturing protein and lipids, it also permits the removal of PCR inhibitors by separating nucleic acids from other compounds.
When the bead-beating and phenol-chloroform steps were combined (protocol 4), they resulted in the highest total DNA concentration and the best quality of extracted DNA (Fig. 1). Nevertheless, the bacterial and fungal DNA, estimated by qPCR, was not significantly elevated compared to protocol 1. Moreover, the majority of bacterial and fungal taxa increased in protocol 4 were also increased in protocol 2 and/or the protocol 3. Subsequently, we added a pre-treatment step using PBS 1× under agitation and filtering on gauze and centrifugation to harvest microbes. Compared to protocol 1, we did not observe a loss of total DNA in lung tissue samples processed with the protocol 5. In contrast, higher concentrations of bacterial and fungal DNA, determined by qPCR, were observed compared to the remaining protocols. These results were expected, as the filtering and centrifugation steps in pre-treatment may diminish the concentration of human DNA and concentrate microorganisms. Nevertheless, the estimation of human DNA actin was not significantly different from other protocols. Although higher concentration of microbial DNA was observed in protocol 5, fewer bacterial and fungal taxa were identified, suggesting that the addition of a pre-treatment resulted in the loss of microorganisms. On the other hand, this may have resulted from the removal of microbial DNA from dead and/or broken-down microorganisms. No methodologies based on sequence analysis of amplicons from microbial markers (such as 16S rRNA and ITS genes) can differentiate living, dead, or ruptured bacteria, as all of these generate the same positive signals 28 . Here, the centrifugation performed during the pre-treatment step contributed to pellet the living microorganisms and remove other particles in the supernatant such as free DNA, resulting in the enrichment of living microbes. Although we did not show that the pre-treatment step decreases the proportion of dead microorganisms, the protocol 5 including the pre-treatment step could be a very important issue. In effect, controversy exists regarding the presence of living commensal bacteria in the lung since the detected DNA sequences may result of the breakdown of microorganisms and not living/reproducing microbial community members.

Conclusion
Collectively, the most relevant findings of our study are that we to show that low microbial biomass samples such as lung tissue samples require particular care to avoid the risk of retrieving highly contaminated sequences. In addition, we found that the DNA extraction protocols affect the detection and abundance of bacterial and fungal taxa, which might have an influence on the interpretation of results, even though a small percentage of the total estimated microbial communities were affected. The bead-beating and Phenol:Chloroform:Isoamyl alcohol steps improve the detection of bacterial and fungal communities through an efficient lysis of microorganisms and a better quality of extracted DNA respectively. Additionally, we found that the addition of a pre-treatment step improves the amplification of microbial DNA, but also it might contribute to eliminate microbial DNA fragment resulting of dead microorganisms in lung tissue samples. Therefore, the microbial communities' profiles determined in samples processed with the protocol 5 that includes the pre-treatment step might be closer to those inhabiting the lung. As lung microbiome studies are needed to determine whether pathogenic relationships between microbiome and lung disease exist, this protocol, including the pre-treatment step described, represent a reproducible alternative to avoid controversy on the live microorganisms inhabiting the lung.

Ethical issues. The Ethics Commission for Studies in Human Subjects of the University of Chile School of
Medicine approved this study under protocol CEISH #092-2013. Informed consent was not obtained because this was a retrospective study of stored samples without identifiers as to the subjects of origin, which is in line with Chilean laws and regulations and assures protection of the subjects. All methods were performed in accordance with the principles of the Declaration of Helsinki for medical research involving human subjects and with the relevant local guidelines and regulations.
Study design and sample collection. Chilean law requires autopsy for accidental deaths and for unexpected infant deaths. They are conducted at the Servicio Médico Legal de Santiago, which is the coroner's office institution for the Metropolitan Area of Chile. Fresh frozen samples are kept stored after the autopsy procedure so they can be used for further analyses. Stored lung samples from legally required autopsies were selected based on unexpected death, storage at − 80 °C within less than 12 h from death, absence of previous admission to the hospital, and absence of known immunocompromising conditions and of obvious pulmonary disease at autopsy examination. Samples from four children aged 4.7 ± 2.8 months without an ascertainable cause of death and www.nature.com/scientificreports/ categorized as sudden infant deaths, and from three adults aged 42.4 ± 17.4 years whose death was categorized as road accident, were sent frozen to our laboratory without possible identifiers to link their identity. For each subject, the right upper lobe was removed using sterile equipment and stored at − 80 °C in a sterile plastic bag until processing for analysis. The samples were removed from the bag and placed on a sterile plate in a biosafety cabinet. The pleura was removed to access untouched tissue using sterile equipment. Small samples were obtained from deep lung tissue, cut into small pieces (0.4 g), and frozen at -80 °C.

DNA extraction protocols.
In this study, we compared five DNA extraction protocols (Table 1)  Protocol 5. This protocol was performed according to the procedure of the protocol 4. However, lung tissue is processed by a "Pre-treatment step", consisting of homogenizing the lung by agitation with a magnetic stirrer in 20 ml of sterile PBS (pH 7.2) in ice pack-covered screw-capped flasks for 30 min, instead of with the Ultra Turrax homogenizer. The homogenate was filtered using sterile gauze and the stirring flasks were washed with sterile PBS to collect any remnants of the specimens. The filtrate was centrifuged and the pellet was reconstituted in 200 µl of sterile PBS, and DNA was extracted as in protocol 4.

DNA extraction controls.
Five blank samples consisting of buffer supplied with the kit were processed together with the lung tissue samples. Four were processed with the Ultra Turrax homogenizer and DNA was extracted according to the protocols 1-4, whereas the fifth was processed with the pre-treatment step according to the protocol 5.
DNA quantitation and quality assessment. DNA concentrations were initially measured using the Qubit double-stranded DNA (dsDNA) BR assay kit on a Qubit fluorometer (INVITROGEN), and dilutions were accordingly prepared and measured. In addition, we determined DNA purity by measuring the concentration of undiluted DNA and absorbance ratios at 260/280 using a NanoDrop 1000 spectrophotometer (THERMO SCIENTIFIC).
Estimation of human, bacterial and fungal DNA levels in DNA extracted from lung tissue samples. Levels of human, bacterial and fungal DNA were assessed using qPCR method by amplifying the human β-actin, the bacterial 16S rRNA, and the 18S rDNA gene of fungi, respectively. qPCRs were carried out in 10 μl reactions containing 2 μl of diluted template (~ 10 ng/μl) or water (negative template control), 2× LightCycler480 SYBR Green I Master (ROCHE DIAGNOSTICS), 0.3 μM each of the forward and reverse primers: 5′-TTG TTA CAG GAA GTC CCT TGCC-3′ and 5′-ATG CTA TCA CCT CCC CTG TGTG-3′ to amplify the human β-actin; the forward 515F 5′-GTG CCA GCMGCC GCG GTAA-3′ and reverse 5′-CTT GTG CGGKCCC CCG YCA ATT C-3′ to amplify the V4 hyper variable region of the 16S rRNA gene of bacteria 30 ; and the forward 5′-TTA GCA TGG AAT AAT RRA ATA GGA -3′ and reverse 5′-TCT GGA CCT GGT GAG TTT CC-3′ to amplify the V4 (partial) and V5 variable regions of the 18S rDNA of fungi 31 , and 2.4 μl of water. PCR reactions were performed on the Roche LightCycler 480 instrument (ROCHE DIAGNOSTIC). All samples were adjusted to 10 ng/µl using the Qubit Fluorometer (INVITROGEN), allowing for comparison of Ct values.
16S rRNA gene and ITS amplification, library construction, and sequencing. For bacterial 16S rRNA gene amplification, primers (forward primer 5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG CCT ACGGGNGGC WGC AG-3′ and reverse primer 5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GGA CTACHVGGG TAT CTA ATC C-3′) spanning the V3/V4 hypervariable regions were used. Internal controls of extraction and amplification were analysed together with the samples. As for fungi, an internal transcribed spacer (ITS) region was amplified. A pre-amplification with primers ITS1-F 5′-TAG AGG AAG TAA AAG TCG TAA-3′ and ITS2-R_KYO2 5′-TTY RCT RCG TTC TTC ATC -3′ spanning the small subunit and the 5.8S region of the rRNA operon, was carried out. A second amplification with internal primers (ITS1-FInt 5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG GGA AGT AAA AGT CGT AAC AAGG-3′, and ITS2_RInt: 5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACAG CTR YGT TCT TCA TCGDT-3′) containing the adapters sequence (in italics), was performed on 10.5 μl of the primary PCR. The resulting products were verified in a 1.4% agarose gel and purified amplicons were quantified using a Qubit Fluorometer (THERMO FISHER SCIENTIFIC). Next, dual indices were attached to both ends of the PCR products using Nextera XT Index Kit (ILLUMINA). Equimolar amounts of DNA per sample were pooled and on a MiSeq desktop sequencer (2 × 300 bp paired-end reads) (ILLUMINA).

Scientific Reports
| (2020) 10:17377 | https://doi.org/10.1038/s41598-020-74137-2 www.nature.com/scientificreports/ Bioinformatic processing. Forward and reverse fastq files containing reads in matched order, free of primer, adapter and linker sequences served as input for the DADA2 pipeline 32 , which analyses the quality profiles, filters and trims to remove Ns, expected errors and low quality tails. After learning the error rates with the DADA2 algorithm (https ://benjj neb.githu b.io/dada2 /tutor ial.html), a dereplication step was used to reduce computation time. Paired reads were merged by aligning denoised forward and reverse reads, and merged reads served to construct the amplicon sequence variant table, to identify and remove chimeric sequences. Before taxonomy assignment, there was a removal of human sequences using bowtie2-2.3.4.2 33 against the reference human genome database GRCh38.p11, using very sensitive parameters. The unaligned reads were used to assign taxonomy, using the Silva reference database for bacteria and Unite database for fungi. Finally, counts were obtained for operational taxonomic units (OTUs), and collapsed to the family level.
Contaminated sequences assessment. The bacterial and fungal sequences clustering with the DADA2 pipeline at 100% identity with those present in the negative controls were removed from the analysis, for each group, respectively. The proportion of removed sequences was calculated for each sample.
Bacterial and fungal community composition, abundance, and diversity analysis. The OTU table was converted into Biom format, using the QIIME pipeline version 1.9.0 34 for composition and absolute and relative abundance analyses, as well as for ecological diversity. For alpha diversity, 1000 rarefactions with replacement were carried out and the Shannon diversity index was calculated. As for beta diversity, variation was assessed using canonical correspondence analysis (CCA), implemented in R version 3.1.0 35 , on a Bray-Curtis dissimilarity matrix.

Statistics.
Pairwise comparisons according to the method and individual were carried out using the qiime script compare_alpha_diveristy.py, with the default non-parametric t-tests. Boxplots were also generated by this script 36 . Fold change in relative abundance of family level taxa between distinct protocols was determined using the R DESEQ2 statistical software 37 within the phyloseq package 38 . Wilcoxon-Mann-Whitney non-parametrical tests for groups of samples were conducted using R version 3.5.1 38 .

Data availability
The sequence data are deposited in EBI Short Read Archive repository (https ://www.ebi.ac.uk/ena) under study accession number PRJEB31011 with accession numbers for bacteria from ERS3088332 to ERS3088370 and for fungi from ERS3088371 to ERS3088409.