Synchronized shift of urine, faeces and saliva

Supplementary Material

Synchronized shift of urine, faeces and saliva microbiotas in bats and natural infection dynamics during seasonal reproduction Muriel Dietrich, Teresa Kearney, Ernest C. J. Seamark, Janusz T. Paweska and Wanda Markotter

ROYAL SOCIETY OPEN SCIENCE

Text S1. DNA extraction protocol Text S2. Illumina sequencing and bioinformatics Text S3. Analysis of bacterial and viral shedding Text S4. Phylogenetic analyses of infectious agents

Table S1. Summary of statistical models used to analyse infection dynamic and microbiota diversity Table S2. Model selection tables. Table S3. Details of the number of bat samples collected at different sample periods

Figure S1. Variation of microbiota diversity

References

1

Text S1. DNA extraction protocol We used a previously validated DNA extraction protocol developed for the Human Microbiome Project [1,2], recently used for bat samples [3]. Care was taken to avoid bacterial DNA contamination by utilizing DNA-free reagents when applicable, filter sterilizing all solutions through a 0.2 µM filter, and working in a PCR-clean hood. Before DNA extraction, faeces samples were weighted and the volume of urine samples was measured. To control for the introduction of contaminating DNA, negative extraction controls (reagents only) were included in the extraction procedure. The protocol includes the addition of the peptidoglycan-degrading enzymes mutanolysin and lysozyme and the use of the QIAmp DNA Micro kit (Qiagen, Valencia, CA), as previously described [3]. We used a 3 hours incubation time at 56°C, the addition of carrier RNA, and faeces samples were centrifuged (2 min. at 604g) and the supernatant collected before the transfer to the Qiagen column. DNA was eluted into 32 µl of buffer AE and stored at -20°C.

2

Text S2. Illumina sequencing and bioinformatics We subjected the DNA extracts and the negative controls (from the field and the extraction) to V3-V4 region 16S rRNA PCR and barcoded Illumina sequencing, following the standardized and optimized Metabiote® protocol developed by Genoscreen (Lille, France). Samples from M. natalensis and R. aegyptiacus were processed separately, with the production of two distinct libraries. Additional negative and positive PCR controls, in the form of an ultrapure water sample and a mock community respectively, were provided by Genoscreen and included in the preparation of each library. The indexed DNA libraries were equimolar pooled and diluted at a final concentration of 4 nM. Paired-end (2 x 250) sequencing was performed using a MiSeq Reagent Nano kit on two different Illumina MiSeq runs, using two lanes each, using PhiX DNA (15%) as a spike-in control for the estimation of the error rate during sequencing. Libraries were de-multiplexed using Casava version 1.8 and raw reads were recovered as FASTQ files. Forward and reverse primers were removed at 100% nucleotide identity by Genoscreen and sequences were quality trimmed when Q score < 30. Alignment of PhiX reads to the reference genome yielded an error rate during sequencing comprised between 1.13-1.42% for read 1 and 1.101.33% for read 2 (depending on the lane). The correct taxonomic assignation of the 11 bacterial species present in the positive control was validated by Genoscreen, using the Greengenes reference database (DeSantis et al. 2006). Samples were then analyzed using MOTHUR v.1.33.3 following the MiSeq SOP Pipeline [4,5]. Assembled reads were quality trimmed based on their length prior to alignment against the MOTHUR-formatted SILVA database. Preclustering of the data was performed using a 4-bp difference, following by the detection and removal of chimeras using the UCHIME algorithm [6]. We then classified sequences using the MOTHUR-formatted version of the RDP training set (v.9), and any unknown, chloroplast, mitochondrial, archaeal, or eukaryotic sequences were removed. Remaining sequences were clustered into phylotypes using a 97% identity threshold. Using both the negative and positive controls, we identified potential exogenous DNA from laboratory reagents, as previously described from low-biomass bat samples [3]. We therefore produced, for each bat species, different datasets corresponding to different levels of exogenous phylotypes removal, following the approach described in Dietrich et al. [3]. Briefly, in the first dataset (D1) all phylotypes were included; in the second (D2) and third (D3) datasets, exogenous phylotypes with a relative abundance > 10% and > 0.1% in the controls respectively were removed; and in the fourth dataset (D4), all exogenous phylotypes were removed. Preliminary analyses confirmed that both bacterial community diversity (ANOVA: p = 0.001) and composition (PERMANOVA: p = 0.001) were different among body habitats in M. natalensis [3]. Therefore, we analyzed separately the microbiota of each body habitat.

3

Text S3. Analysis of bacterial and viral shedding The second half of DNA extracts from saliva and faeces samples were used as templates for herpesvirus and adenovirus detection respectively, using previously published nested-PCR protocols both targeting the DNA polymerase gene [7,8]. Briefly, both PCRs were ran in 25µl total volume containing 12.5µl of DreamTaq Green PCR Master Mix (Thermo Fisher Scientific) and 1 µM of each primer. After electrophoresis on a 2% agarose gel, PCR products of the approximate anticipated size were submitted to Sanger sequencing (Inqaba). For the detection of Leptospira bacteria in urine samples, a probespecific real-time PCR, targeting the 16S rRNA gene, was performed as previously described [9]. All PCRs were ran with a negative (PCR mix only) and positive controls (i.e. DNA from a Human herpesvirus 1 cell culture, DNA from an adenovirus-positive bat (Neoromicia spp.) sample, DNA from a Leptospira interrogans culture). For Leptospira positive samples, a partial fragment of the 16S rRNA gene was then amplified, as described in Dietrich et al. [10] and sequenced.

4

Text S4. Phylogenetic analyses of infectious agents Nucleotide sequences of Leptospira, adenoviruses and herpesviruses were checked individually with ChromasLite 2.01 (Technelysium Pty, South Brisbane, Australia) and then assembled and aligned with reference sequences using CLC Sequence Viewer 7.8.1 (CLC Bio, Aarhus, Denmark). Careful manual checking of chromatograms did show double peaks for herpesviruses found in R. aegyptiacus, thus we retained only the mono-infected samples for the phylogenetic analyses. Phylogenetic trees were constructed using BEAST v.1.8.4. [11], using a GTR+I+G substitution model, a strict clock and a constant population size coalescent tree prior. Analyses were run for 100 x106 (for viruses) and 300 x 106 (for Leptospira) generations, sampling every 1,000 generations, with the initial 10% discarded as burn-in. TRACER v.1.635 was then used to verify that the effective sample size of each parameter was higher than 200. The sampled posterior trees were summarized using TREEANNOTATOR v.1.8.4 to generate a maximum clade credibility tree. Sequences produced in this study have been submitted to GenBank database under accession numbers MG680317-MG608402.

5

Table S1. Summary of statistical models used to analyse infection dynamic and microbiota diversity. N gives the number of samples included in each analysis. Full model includes all the tested variables. “k” is the number of model’s coefficients.

Microbiota diversity

Infection dynamics

Variable of interest Leptospira shedding prevalence in M. natalensis Adenovirus shedding prevalence in M. natalensis Herpesvirus shedding prevalence in M. natalensis Herpesvirus shedding prevalence in R. aegyptiacus Inverse Simpson index in urine in M. natalensis Inverse Simpson index in faeces in M. natalensis Inverse Simpson index in saliva in M. natalensis Inverse Simpson index in saliva in R. aegyptiacus

Model nb.

N

Full model

k

Significant variables

GLM1

73

Session+Sex+Repro+Age

6

Age

GLM2

90


6

Session

GLM3

103


6

Session+Age

GLM4

276


7

Session+Age+ Repro

GLM5

73

7

Lepto

GLM6

90

7

Sex

GLM7

103

7

Session

GLM8

276

8

Session+Sex

Session+Sex+Repro+Age +Lepto Session+Sex+Repro+Age +AdV Session+Sex+Repro+Age +Herpes Session+Sex+Repro+Age +Herpes

6

Table S2. Model selection tables. Output were obtained with the MuMin package in R (dredge function) for the eight GLMs analyzed (in Table S1). Selected models (after verification of the effect of each variable) are highlighted in red.

Infection dynamics GLM1 - Candidate models Age Repro + Age Age + Sex Age + Repro + Sex Session + Age Session + Age + Repro null Session + Age + Sex Session + Age + Sex + Repro Session Repro Sex Sex + Repro Session + Sex Session + Repro Session + Repro + Sex

GLM2 - Candidate models Session + Sex Session + Sex + Repro Session Session + Repro Session + Sex + Age Session + Sex + Age + Repro Session + Age Session + Age + Repro Age null Age + Sex Age + Repro Sex Age + Repro + Sex Repro Repro + Sex

k 2 3 3 4 4 5 1 5 6 3 2 2 3 4 4 5

k 4 5 3 4 5 6 4 5 2 1 3 3 2 4 2 3

logLik -32.257 -31.594 -31.996 -31.103 -31.711 -30.994 -35.682 -31.466 -30.535 -34.407 -35.577 -35.682 -34.789 -34.359 -34.367 -33.375

logLik -21.747 -20.791 -23.200 -22.512 -21.747 -20.791 -23.200 -22.512 -54.266 -55.799 -53.842 -53.969 -55.502 -53.510 -55.680 -54.847

AICc 68.7 69.5 70.3 70.8 72.0 72.9 73.4 73.8 74.3 75.2 75.3 75.5 75.9 77.3 77.3 77.6

ΔAIC 0.00 0.85 1.66 2.11 3.33 4.20 4.74 5.14 5.66 6.48 6.64 6.85 7.24 8.62 8.64 8.96

AICc 52.0 52.3 52.7 53.5 54.2 54.6 54.9 55.7 112.7 113.6 114.0 114.2 115.1 115.5 115.5 116.0

ΔAIC 0.00 0.33 0.71 1.53 2.24 2.63 2.91 3.77 60.71 61.68 62.00 62.25 63.18 63.53 63.53 64.01

weight 0.317 0.207 0.139 0.110 0.060 0.039 0.030 0.024 0.019 0.012 0.011 0.010 0.008 0.004 0.004 0.004

weight 0.251 0.212 0.175 0.117 0.082 0.067 0.059 0.038 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

7

GLM3 - Candidate models Session + Age Session + Age + Sex Session + Age + Repro Session + Age + Repro + Sex Sex + Age Age Age + Repro Age + Repro + Sex Session Session + Repro Session + Sex Session + Sex + Repro null Sex Repro Repro + Sex

k 4 5 5 6 3 2 3 4 3 4 4 5 1 2 2 3

logLik -61.374 -60.760 -61.113 -60.760 -64.417 -65.798 -64.957 -64.417 -67.459 -66.736 -67.051 -66.736 -71.156 -70.975 -715155 -70.721

AICc 131.2 132.1 132.8 134.4 135.1 165.7 136.2 137.2 141.2 141.9 142.5 144.1 144.4 146.1 146.4 147.7

ΔAIC 0.00 0.98 1.69 3.24 3.92 4.56 5.00 6.09 10.00 10.72 11.35 12.93 13.20 14.91 15.27 16.53

weight 0.380 0.232 0.163 0.075 0.054 0.039 0.031 0.018 0.003 0.002 0.001 0.001 0.001 0.000 0.000 0.000

GLM4 - Candidate models Session + Age + Repro Session + Repro Session + Age + Repro + Sex Session + Repro + Sex Session + Age Session + Age + Sex Session Session + Sex Age Repro + Age Age + Sex Repro Repro + Age + Sex Repro + Sex null Sex

k 6 5 7 6 5 6 4 5 2 3 3 2 4 3 1 2

logLik -161.889 -163.845 -161.858 -163.772 -164.903 -164.724 -170.214 -169.854 -180.261 -180.071 -180.256 -181.788 -180.059 -181.783 -183.820 -183.805

AICc 336.1 337.9 338.1 339.9 340.0 341.8 348.6 349.9 364.6 366.2 366.6 367.6 368.3 369.7 369.7 371.7

ΔAIC 0.00 1.82 2.04 3.77 3.94 5.67 12.48 13.84 13.84 30.14 30.51 31.53 32.17 33.56 33.56 35.56

weight 0.473 0.190 0.170 0.072 0.066 0.028 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

8

Microbiota diversity GLM5 - Candidate models Lepto Lepto + Sex Lepto + Repro Lepto + Age Lepto + Session + Sex Lepto + Repro + Sex Lepto + Age + Sex Lepto + Repro Lepto + Repro + Age Lepto + Session + Repro Lepto + Session + Sex + Age Lepto + Session + Sex + Repro Lepto + Sex + Repro + Age Lepto + Session + Age null Lepto + Session + Repro + Age Sex Age Lepto + Session + Repro + Age + Sex Repro Session + Sex Session Age + Sex Repro + Sex Repro + Age Session + Repro Session + Repro + Sex Session + Age + Sex Session + Age Repro + Age + Sex Session + Repro + Age Session + Repro + Age + Sex

k 3 4 4 4 6 5 5 5 5 6 7 7 6 6 2 7 3 3 8 3 5 4 4 4 4 5 6 6 5 5 6 7

logLik -247.444 -246.824 -247.045 -247.441 -245.386 -246.800 -246.806 -246.845 -247.036 -245.924 -245.132 -245.345 -246.784 -246.841 -251.548 -245.779 -250.991 -251.010 -245.096 -251.307 -249.148 -250.619 -250.650 -250.846 -250.885 -249.899 -248.957 -249.122 -250.327 -250.504 -249.797 -248.934

AICc 501.2 502.2 502.7 503.5 504.0 504.5 504.5 504.6 505.0 505.1 506.0 506.4 506.8 507.0 507.3 507.3 508.3 508.4 508.4 509.0 509.2 509.8 509.9 510.3 510.4 510.7 511.2 511.5 511.5 511.9 512.9 513.6

ΔAIC 0.00 1.00 1.44 2.23 2.81 3.26 3.27 3.35 3.73 3.88 4.75 5.18 5.60 5.72 6.03 6.04 6.04 7.13 7.20 7.73 7.95 8.59 8.65 9.04 9.12 9.46 9.95 10.28 10.31 10.67 11.63 12.35

weight 0.242 0.147 0.118 0.079 0.059 0.047 0.047 0.045 0.037 0.035 0.022 0.018 0.015 0.014 0.012 0.012 0.012 0.007 0.007 0.005 0.005 0.003 0.003 0.003 0.003 0.002 0.002 0.001 0.001 0.001 0.001 0.001

9

GLM6 - Candidate models Sex Repro Sex + Repro Sex + Age Sex + AdV Repro + AdV Repro + Age Session + Sex Sex + Repro + Age Sex + Repro + AdV Sex + AdV + Age Session + Repro Session + Sex + Age Repro + AdV + Age Session + Sex + Age Session + Sex + Repro Sex + Repro + AdV + Age Session + Repro + Age null Session + Repro + Age Session + Repro + Sex + AdV Session + Sex + Age + AdV Session Session + Repro + Sex + Age Session + Age Age AdV Session + Repro + Age + AdV Session + Repro + Age + Sex + AdV Session + AdV Session + Age + AdV Age + AdV

k 3 3 4 4 4 4 4 5 5 5 5 5 6 5 6 6 6 6 2 6 7 7 4 7 5 3 3 7 8 5 6 4

logLik -185.315 -185.955 -185.284 -185.297 -185.314 -185.945 -185.952 -185.030 -185.240 -185.284 -185.295 -185.703 -184.709 -185.943 -184.989 -184.997 -185.239 -185.545 -190.075 -185.696 -184.635 -184.685 -188.212 -184.905 -187.544 -189.956 -190.015 -185.543 -184.553 -188.190 -187.521 -187.521

AICc 376.9 378.2 379.0 379.1 379.1 380.4 380.4 380.8 381.2 381.3 381.3 382.1 382.4 382.6 383.0 383.0 383.5 384.1 384.3 384.4 384.6 384.7 384.9 385.2 385.8 386.2 386.3 386.5 386.9 387.1 388.1 388.2

ΔAIC 0.00 1.28 2.13 2.16 2.19 3.45 3.470 3.86 4.29 4.37 4.39 5.21 5.52 5.69 6.08 6.10 6.58 7.19 7.38 7.38 7.73 7.83 7.99 8.27 8.89 9.28 9.40 9.54 9.97 10.19 11.15 11.29

weight 0.255 0.134 0.088 0.087 0.085 0.045 0.045 0.037 0.030 0.029 0.028 0.019 0.016 0.015 0.012 0.012 0.009 0.007 0.006 0.006 0.005 0.005 0.005 0.004 0.003 0.002 0.002 0.002 0.002 0.002 0.001 0.001

10

GLM7 - Candidate models Session Session + Herpes Session + Repro Session + Sex Session + Age Session + Repro + Herpes Session + Age + Herpes Session + Sex + Herpes Session + Repro + Sex Session + Repro + Age Session + Sex + Age Session + Sex + Repro + Herpes Session + Herpes + Repro + Age Session + Herpes + Age + Sex Session + Repro + Age + Sex Session + Repro + Age + Sex + Herpes Sex + Age Sex + Age + Herpes Sex + Age + Repro Sex Sex + Age + Repro + Herpes Sex + Repro Sex + Herpes Age + Repro Sex + Repro + Herpes Repro + Herpes + Age Repro null Age Herpes + Repro Herpes Herpes + Age

k 4 5 5 5 5 6 6 6 6 6 6 7 7 7 7 8 4 5 5 3 6 4 4 4 5 5 3 2 3 4 3 4

logLik -333.571 -333.202 -333.277 -333.534 -333.566 -333.813 -333.115 -333.140 -333.182 -333.209 -333.533 -332.719 -332.812 -333.104 -333.145 -332.718 -344.167 -344.088 -344.132 -346.808 -344.054 -346.506 -346.722 -347.049 -346.440 -347.038 -349.295 -350.845 -350.007 -349.134 -350.685 -349.992

AICc 675.5 677.0 677.2 677.7 677.8 678.5 679.1 679.2 679.2 679.3 679.9 380.6 680.8 681.4 681.5 683.0 696.7 698.8 698.9 699.9 701.0 701.4 701.9 702.5 703.5 704.7 704.8 705.8 706.3 706.7 707.6 708.4

ΔAIC 0.00 1.47 1.62 2.14 2.20 2.95 3.56 3.61 3.69 3.74 4.39 5.07 5.25 5.84 5.92 7.42 21.19 23.25 23.33 24.31 25.43 25.87 26.30 26.96 27.95 29.14 29.28 30.26 30.71 31.13 32.06 32.84

weight 0.259 0.124 0.115 0.089 0.086 0.059 0.044 0.043 0.041 0.040 0.029 0.021 0.019 0.014 0.013 0.006 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

11

GLM8 - Candidate models Session + Sex Session + Sex + Herpes Session + Sex + Age Session + Sex + Repro Session Session + Sex + Herpes + Age Session + Sex + Herpes + Repro Session + Sex + Age + Repro Session + Herpes Session + Age Session + Repro Session + Sex + Age + Herpes + Repro Session + Herpes + Repro Session + Herpes + Age Session + Repro + Age Sex Sex + Age Session + Herpes + Repro + Age null Sex + Repro Sex + Herpes Age Sex + Herpes + Age Sex + Age + Repro Repro Herpes Sex + Herpes + Repro Herpes + Age Repro + Age Sex + Herpes + Repro + Age Herpes + Repro Herpes + Repro + Age

k 6 7 7 7 5 8 8 8 6 6 6 9 7 7 7 3 4 8 2 4 4 3 5 5 3 3 5 4 4 6 4 5

logLik -259.027 -258.952 -258.970 -259.026 -261.168 -258.915 -258.951 -258.958 -261.129 -261.160 -261.166 -258.895 -261.121 -261.126 -261.149 -265.718 -265.234 -261.109 -267.565 -265.608 -265.656 -266.767 -265.115 -265.211 -267.301 -267.511 -265.524 -266.640 -266.759 -265.096 -267.213 -266.635

AICc 530.4 532.3 532.4 532.5 532.6 534.4 534.4 534.5 534.6 534.6 534.6 536.5 536.7 536.7 536.7 537.5 538.6 538.8 539.2 539.4 539.5 539.6 540.5 540.6 540.7 541.1 541.3 541.4 541.7 542.5 542.6 543.5

ΔAIC 0.00 1.96 1.99 2.10 2.19 4.00 4.08 4.09 4.20 4.27 4.28 6.10 6.29 6.30 6.35 7.16 8.25 8.39 8.81 9.00 9.09 9.26 10.09 10.28 10.33 10.74 10.90 11.06 11.30 12.14 12.21 13.13

weight 0.286 0.107 0.106 0.100 0.095 0.039 0.037 0.037 0.035 0.034 0.034 0.014 0.012 0.012 0.012 0.008 0.005 0.004 0.003 0.003 0.003 0.003 0.002 0.002 0.002 0.001 0.001 0.001 0.001 0.001 0.001 0.000

12

Table S3. Details of the number of bat samples collected at different sample periods.

M. natalensis Oral Faeces Urine Subtotal R. aegyptiacus Oral Total

September

October/November

January

April

Total

19 23 4 46

36 39 45 120

48 28 24 100

0 0 0 0

103 90 73 266

47 93

95 215

38 138

96 96

276 542

13

Figure S1. Variation of microbiota diversity

14

References 1.

Yuan S, Cohen DB, Ravel J, Abdo Z, Forney LJ. 2012 Evaluation of Methods for the Extraction and Purification of DNA from the Human Microbiome. PLoS One 7, e33865. (doi:10.1371/journal.pone.0033865)

2.

Pearce MM et al. 2014 The Female Urinary Microbiome : a Comparison of Women with and without Urgency Urinary Incontinence. MBio 5, 1–12. (doi:10.1128/mBio.01283-14.Editor)

3.

Dietrich M, Kearney T, Seamark ECJ, Markotter W. 2017 The excreted microbiota of bats : evidence of niche specialisation based on multiple body habitats. FEMS Microbiol. Lett. 364, fnw284. (doi:10.1093/femsle/fnw284)

4.

Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. 2013 Development of a DualIndex Sequencing Strategy and Curation Pipeline for Analyzing Amplicon Sequence Data on the MiSeq. Appl. Environ. Microbiol. 79, 5112–5120. (doi:10.1128/AEM.01043-13)

5.

Schloss PD et al. 2009 Introducing mothur : Open-Source, Platform-Independent, CommunitySupported Software for Describing and Comparing Microbial Communities. Appl. Environ. Microbiol. 75, 7537–7541. (doi:10.1128/AEM.01541-09)

6.

Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. 2011 UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27, 2194–2200. (doi:10.1093/bioinformatics/btr381)

7.

VanDevanter DR, Warrener P, Bennett L, Schultz ER, Coulter S, Garber RL, Rose TM, Black J. 1996 Detection and Analysis of Diverse Herpesviral Species by Consensus Primer PCR. J. Clin. Microbiol. 34, 1666–1671.

8.

Li Y et al. 2010 Host Range, Prevalence, and Genetic Diversity of Adenoviruses in Bats. J. Virol. 84, 3889. (doi:10.1128/JVI.02497-09)

9.

Smythe LD, Smith IL, Smith GA, Dohnt MF, Symonds ML, Barnett LJ, McKay DB. 2002 A quantitative PCR (TaqMan) assay for pathogenic Leptospira spp. BMC Infect. Dis. 2, 13.

10.

Dietrich M, Wilkinson DA, Soarimalala V, Goodman SM, Dellagi K, Tortosa P. 2014 Diversification of an emerging pathogen in a biodiversity hotspot: Leptospira in endemic small mammals of Madagascar. Mol. Ecol. 23, 2783–96. (doi:10.1111/mec.12777)

11.

Guindon S, Gascuel O. 2003 A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst. Biol. 52, 696–704. (doi:10.1080/10635150390235520)

12.

Crawley MJ. 2007 The R Book. John Wiley. Chichester, UK.

15