Rabbit model and housing conditions
The study was performed on 33 rabbits from a line raised at the Polytechnic University of Valencia (Spain) resulting from a cross between New Zealand White and California rabbits. Animals were weaned at 30 days of age and housed in separate individual cages where they had ad libitum access to food and water. The food that was used in the diet was standard, antibiotic-free and coccidiostatic-free (Cunicebal Retirada, Nanta, Spain). Three independent experiments were performed. The number of animals per experiment varies due to the different number of cages or litters available. All experiments consisted on monitoring ERE signs and measuring the food consumption. During the first experiment 2 out of 9 animals developed ERE after weaning. In the second experiment 5 out of 14 rabbits developed ERE, and in the last experiment, 4 out of 10 rabbits developed ERE. Caecotroph samples from all experiments were analyzed together in the longitudinal analysis presented in the results. For each experiment, rabbits from different maternal origins were used (N = 2–4 litters per experiment).
Identification of rabbits developing ERE and defining the day of ERE onset
A sharp decrease in food consumption is considered the first sign that can be detected in rabbits that develop ERE [5, 6]. For this reason, food consumption was measured daily in order to detect the onset of the disease. The day of disease onset was considered to be that one in which a rabbit diminished food consumption > 50% compared to the previous day. The threshold of 50% was chosen since a few healthy rabbits diminished food consumption on specific days (up to 46%, Additional file 1) without developing any ERE signs. On the contrary, ERE rabbits developed ERE compatible signs (prostration, diarrhea, tympanism) following the detected sharp decrease in food intake (Additional file 2). Moreover, none of the ERE rabbits increased food intake after ERE initiation (Additional file 3), while those healthy rabbits that diminished food consumption on a particular day for an undetermined reason, later recovered their food intake (Additional file 3). For ethical reasons, rabbits developing ERE were sacrificed when they ate during two consecutive days less than 10 g of food, except for four rabbits, which died of ERE symptoms before fulfilling this criterion. Healthy rabbits were followed for an additional week before euthanasia (total of 3 weeks after weaning) to ensure that they would not develop any ERE compatible signs.
To confirm ERE development, a necropsy was performed on all rabbits. This necropsy was performed by co-author Jorge Martínez, a veterinarian specialist in anatomical pathology (Diplomate by the European College of Veterinary Pathologists). Examination was performed in order to identify gross lesions compatible with ERE: distension of the stomach and intestine, caecum filled with gas, liquid or impacted with solid/dry content, and presence of abundant quantities of mucus in the colon (Additional file 2). Gross lesions compatible with ERE were observed in all rabbits that developed a sharp decrease in food consumption and ERE compatible signs (prostration, diarrhea) (Additional file 2). However, no lesions suggesting the presence of other intestinal disorders or other systemic diseases were detected in any of the ERE rabbits.
Sample collection and DNA extraction
In order to assess microbiota composition of rabbits, caecotrophs were collected the day of weaning, the day after weaning and then every other day during the next 20 days. Caecotrophs are soft feces excreted by rabbits that have a similar chemical and microbiological composition to that of the caecal contents, where major pathological changes occur during ERE [13]. Caecotroph samples were collected using polyvinyl chloride (PVC) neck collars to prevent their ingestion by rabbits. Caecotrophs were kept at 4 °C during approximately 3 h and aliquots of these caecotroph samples (approximately 200 mg) were stored at −80 °C afterwards. Due to the symptomatology or rabbit death, we could only obtain caecotroph samples after ERE onset from 7 out of the 11 ERE rabbits. For comparison, we did also analyze caecotroph samples from 10 healthy littermate controls. For comparison with a mature adult microbiota (Figure 2), caecotroph samples from 4 mothers were collected the day of weaning and 7 days after weaning. DNA was extracted using a QIAamp® DNA Stool Mini kit (QIAGEN) with a previous step of mechanical disruption to improve cell lysis as we have previously described [14]. Briefly, cells were resuspended in 1.4 mL of ASL buffer and 500 μL of 0.1 mm glass beads and tubes were vortex at maximum speed for 5 min prior to the initial incubation for heat and chemical lysis at 95 °C for 7 min. Subsequent steps of the DNA extraction followed the QIAamp kit protocol.
16S rRNA gene amplification and high-throughput sequencing
The V3–V4 region of the 16S rRNA gene was amplified and sequenced using the MiSeq platform from Illumina, as described in the manual for “16S Metagenomic Sequencing Library Preparation” of the MiSeq platform (Illumina). Briefly, for each sample, a 25 μL reaction was prepared containing 12.5 ng of DNA, 12.5 μL 2× KAPA HiFi Hot Start Mix, and 0.2 mM of primers. Water was added to complete the volume of the reaction. Cycling conditions were 95 °C for 3 min, and 25 cycles of 95 °C for 30 s, 55 °C for 30 s and 72 °C for 30 s, and a final elongation cycle at 72 °C for 5 min. The amplification was confirmed through electrophoresis by loading 4 μL of the PCR reaction on a 1.6% agarose gel. Subsequently, the PCR product was purified with the AMPure XP beads as described in the Illumina protocol. Next, a limited-cycle PCR reaction was performed to amplify the DNA and add index sequences on both ends of the DNA, thus enabling dual-indexed sequencing of pooled libraries. Index PCR consisted of a 50 μL reaction containing 5 μL of the DNA obtained from the previous PCR, 25 μL of 2× KAPA HiFi Hot Start Mix, and 5 μL of forward and reverse indexed primers. Temperature conditions were the same as for the first reaction, but the number of cycles was reduced to 8. The obtained PCR product was purified with the AMPure XP beads following the manufacture’s protocol. An equal amount of the purified DNA was taken from each sample for pooling. Each pool of samples (N = 96) was sequenced following Illumina recommendations.
16S rRNA sequencing analysis
Sequences were processed using Mothur v1.35 [15] as we have previously described [14], with some modifications. Initial trimming by quality was performed on paired ends of sequences before joining them into a single read. Parameters used for trimming included elimination of sequences shorter than 200 bp or that contained homopolymers longer than 8 bp or undetermined bases. Using the base quality scores, which range from 0 to 40 (0 being ambiguous base), sequences were trimmed using a sliding-window technique, such that the minimum mean quality score over a window of 50 bases never dropped below 25. Sequences were trimmed from the 3′ end until this criterion was met. Sequences were aligned to the 16S rRNA gene using as a template the SILVA reference alignment. Potential chimeric sequences were removed using Uchime algorithm. To minimize the effect of pyrosequencing errors in overestimating microbial diversity [16], rare abundance sequences that differ in up to two nucleotides from a high abundant sequence were merged to the high abundant sequence using the pre.cluster option in Mothur. An average of 17,297 sequences were obtained per sample. Since different number of sequences per sample could lead to a different diversity (i.e. more Operational Taxonomic Units—OTUs could be obtained in those samples with higher coverage), in order to compare the diversity of different caecotroph samples, we rarefied all samples to the number of sequences obtained in the sample with the lowest number of sequences (10 997). In other words, 10 997 sequences were randomly selected from each sample for subsequent analysis: taxonomic characterization and OTUs identification. Sequences with distance-based similarity of 97% or higher were grouped into the same OTU using the VSEARCH abundance based greedy clustering (AGC) method. OTU-based microbial diversity was estimated by calculating the Shannon diversity index. Each sequence was classified using the Bayesian classifier algorithm with the bootstrap cutoff of 60% [17]. In most cases classification could be assigned to the genus level. When it was not possible to classify a sequence to a certain taxonomic level, it was assigned as “Unclassified” followed by the upper taxonomic level. Sequences have been deposited in NCBI under the accession number PRJNA509120.
Principal coordinate of analysis
In order to compare the overall microbiota similarity between different caecotroph samples, we calculated for every pair of samples the UniFrac phylogenetic distance. This distance was calculated by generating a phylogenetic tree containing all the 16S rRNA representative sequences of all the OTUs identified in the samples under comparison. Subsequently, the UniFrac distance was calculated as the fraction of the total branch length of the tree, which is not shared between the two samples under comparison. UniFrac distance values range from 0 to 1, being 1 the distance obtained between a pair of samples that do not share any branch of the phylogenetic tree (their microbiota is totally different and do not share any bacterial lineage), and 0 the distance obtained between a pair of samples that have exactly the same microbiota. To perform UniFrac analysis we first inferred a phylogenetic tree using clearcut, on the 16S rRNA sequence alignment generated by Mothur. Both unweighted UniFrac, which tabulates the presence or absence, but not the proportion, of different bacterial lineages within the tree, and weighted UniFrac, which takes into account the proportion of each bacterial lineage in each sample, were run using Mothur. Principal Coordinates of Analysis were performed on the resulting distance matrices with Mothur.
Isolation of Clostridium cuniculi
The C. cuniculi strain was isolated from a −80 °C frozen aliquot of a caecotroph sample collected (as described in the section “Sample collection and DNA extraction”) from one of the rabbits suffering from ERE indicated in Figure 1. The chosen sample was collected after ERE onset, and was characterized by 16s rRNA high-throughput sequencing as the one that contains the highest abundance of the OTU172 (the Clostridium OTU of interest, later defined as C. cuniculi). An aliquot of the chosen caecotroph sample was thawed under anaerobic conditions. 10-fold dilutions of the sample in pre-reduced phosphate buffer saline (PBS) were plated on TSA blood agar plates and incubated for 48 h under anaerobic conditions at 37 °C. Subsequently, a colony PCR using specific primers for the OTU172 was performed, as described below.
Specific primers were designed using the Primrose program. Briefly, a FASTA file with 16S rDNA sequences that are representative for every OTU detected in the caecotroph sample used for isolation were utilized to construct a database, which is used for the design of primers. This approach allows to select primers specific for the 16S rRNA sequence of the desired OTU, and that will not amplify any other OTUs present in the sample used for isolation. The maximum length of the primers was set to 20 nucleotides, while no ambiguous nucleotides were accepted. This analysis led to the design of the pair of primers: OTU172-F: AGGCGGACTTTTAAGTGAGAT and OTU172-R: CGTCAGTTACAGTCCAGAGAAT.
For each colony PCR, a 25 µL reaction was prepared containing: 1 μL of 1 bacterial colony resuspended in 10 μL of PBS, 2.5 µL 10× Standard Taq Reaction Buffer (New England BioLabs), 0.25 mM of deoxynucleotide triphosphates (dNTPs), 2.5 U of Taq DNA Polymerase (New England BioLabs) and 0.2 mM of primers. Cycling conditions were 94 °C for 5 min, and 35 cycles of 94 °C for 30 s, 56 °C for 30 s and 68 °C for 30 s, and a final elongation cycle at 68 °C for 5 min. Positive colonies were re-amplified using the primers used for 16S high-throughput sequencing, as described above. Sequences of amplified products were obtained through Sanger reaction and capillary electrophoresis. The obtained sequences were compared with the sequences that define the OTU172 using blastn in order to verify that the isolated colony was indeed the OTU172 (100% similarity to sequences contained within the OTU172). The colony that was identified as the OTU172 was regrown on TSA blood agar plates for 48 h under anaerobic conditions. Subsequently, the bacterial culture was resuspended in PBS containing 15% glycerol and stored at −80 °C.
Experimental reproduction of ERE
Intestinal contents (caecum and small intestine) from 4 rabbits with ERE were diluted 1:2 in phosphate buffer saline containing 15% glycerol and 0.1% cysteine (previously reduced overnight in an anaerobic chamber) to preserve bacterial viability. Resuspended intestinal contents were frozen at −80 °C and thawed the day of inoculation. C. cuniculi was grown on TSA blood agar plates for 48 h and resuspended, the day of the inoculation, in pre-reduced PBS-glycerol-cysteine at a concentration of 107 CFUs/mL. 5–6 weeks old SPF rabbits (N = 10 per group) were nasogastrically inoculated with either 1 mL of the resuspended intestinal contents, C. cuniculi or PBS-glycerol-cysteine (uninfected control). 2 rabbits were housed per cage. Rabbits from different groups were located in different rooms in order to avoid contamination.
Sequencing of the C. cuniculi genome
To determine genes and functions encoded by C. cuniculi, genomic sequencing was performed. The isolated strain was grown on TSA blood agar in anaerobic conditions at 37 °C. Colonies were collected and resuspended in 1 mL of PBS. Bacterial cells were centrifuged at 13 000 g during 20 min and supernatant was discarded. DNA was extracted with the QIAamp ® Fast DNA Stool Mini kit as described previously.
Genomic DNA was prepared for sequencing using Nextera XT DNA Sample Preparation Kit. Briefly, 5 μL of Amplicon Tagment Mix was added to 10 μL of Tagment DNA Buffer. The amount of DNA added to the mix was 1 ng. The mix was incubated at 55 °C during 5 min. Neutralization of the reaction was achieved by adding 5 μL of Neutralize Tagment Buffer. For the step of amplification, 15 μL of Nextera PCR Mix was added. The cycling conditions were 72 °C during 3 min, 95 °C 30 s followed by 12 cycles of 95 °C during 10 s, 55 °C 30 s and 72 °C 30 s. A final elongation step was performed at 72 °C during 5 min. The obtained PCR product was purified with the AMPure XP Beads as described in the manufacturer’s protocol. Sequencing was performed using the MiSeq platform following Illumina recommendations. The genome sequence has been deposited in NCBI under the accession number PRJNA509119.
Analysis of the C. cuniculi genome
Characterization of the strain of interest included the assembly and annotation of the data obtained by genome sequencing, evolutionary placement of the strain into the Clostridiales order, identification of the species to which it belongs, and functional annotation of the identified genes. First, sequences obtained by Illumina were processed in order to eliminate bad quality sequences. Briefly, adapter sequences were removed from the raw data using Cutadapt v. 1.12 [18]. Quality filtering was performed using UrQt v. 1.0.18 [19]. Only reads with a size of 75 bp or higher were further processed to avoid possible misclassification of short reads. Cleaned genomic data was assembled with SPAdes v. 3.7.1 [20] using the “careful” algorithm to improve the contig reconstruction. SPAdes resulted in a draft genome assembly of 160 contigs, of which 106 had a contig length larger than 1000bps. The average base pair coverage was 27.17 ± 2.97 X. Open-reading frames (ORFs) were identified and annotated using PROKKA v. 1.12 [21]. To maximize the annotation, ORFs were translated into amino acids and queried against 3 independent databases: Pfam v.27.0 [22], EggNOG v.4.5 [23], and KEGG [24]. Annotation was performed using HMMer v.3.1.2 (E value < 0.05, minimum coverage 0.5) [25].
Based on a recently published article [26] the genome draft of the isolate was compared against the reference genomes of the Clostridiales order, obtained from the GenBank RefSeq complete genome database (last accessed July 2018). To do so, we first constructed the core-genome of 98 genomes classified as Clostridiales. To build it, we started by assigning relations of orthology between pairs of genomes. Orthologs were defined as the bidirectional best hits between pairs of genes, using end-gap free global alignment (as in Touchon et al. [27]). Hits with less than 40% similarity in amino acid sequence or more than 20% difference in protein length were discarded. A total of 114 core-genes were identified for the whole Clostridiales order. A sequence concatenate of the core genome at amino acid level was built and aligned using Mafft v.7.15 [28]. Predicted genes from the isolate were aligned using Mafft against the Multiple Sequence Alignment concatenate. The Clostridiales phylogeny was reconstructed with IQ-tree using the concatenate multiple protein sequence alignment of the core-genes [29]. An outgroup of unrelated species was included to root the tree. Evolutionary placement algorithm was performed to locate the most probable phylogenetic location of the isolated strain. The best phylogenetic location was considered only when the Maximum Likelihood Weight was higher than 0.95 (which translates into a probability of misplacement lower than 0.05).
Subsequently, the genome draft of the isolate was compared against the reference sequences of the Clostridium genus, the most likely phylogenetic location. The core-genome of the Clostridium genus was constructed as previously described, resulting in a total of 229 core-genes associated to the Clostridium genus, including the isolate. A sequence concatenate of the core genome at amino acid level was built and aligned using Mafft. Core-genes were then translated to amino acids. The Clostridium genus phylogeny was reconstructed with IQ-tree, as previously described, using the LG + G + I model and 100 bootstrap iterations [29]. The average nucleotide identity of the core genome (cgANI) was calculated using an in-house script. In brief, the core genome of the isolate was fragmented in non-overlapping DNA sequences of 1000 bp and queried against the core-genome of all other references in that dataset. The percentage of identity was calculated for each fragment. The cgANI was then calculated as the mean identity of all fragments. We considered 2 genomic sequences to belong to the same species if the cgANI between them was higher than 95%.
We also defined the pan-genome of the Clostridium genus as the repertoire of gene families present in the clade. The pan-genomic reconstruction was also performed following the approach from Touchon et al. [27].
To identify possible toxins, the exotoxin database DBETH was used [30]. Pan-genome family-representative proteins were queried against DBETH using BlastP v.2.2.26. Significant hits (E-value < 0.01, coverage > 0.5) encoded by C. cuniculi were further characterized using BlastP against the Blast non-redundant protein database searching for specific domains and superfamily domains. Only hits that were confirmed as putative toxins using both the DBETH and the Blast nr database are shown (Additional file 4).
In addition, analysis of the Clostridium pan-genome, allow us to identify proteins encoded by C. cuniculi that are also encoded by other Clostridium species (i.e. presence/absence of potential virulence factors shown in Figure 4A).
Statistical analysis
In order to determine statistically significant differences in the relative abundance of different taxa and OTUs between the group of healthy rabbits and those that developed ERE, the non-parametric Wilcoxon test was applied using wilcox.test function in the “stats” R package. Taxa and OTUs with less than 5 counts in both groups were not included in the analysis. To adjust for multiple hypothesis testing, we used the FDR approach by Benjamini and Hochberg [31] implemented in the fdr.R package. Taxa with a p < 0.05 and FDR < 0.1 were considered significant. First, we investigated differences between samples of healthy rabbits and those that developed ERE, collected after ERE onset. Since samples of ERE rabbits, after they developed ERE, were collected on the day 11 after weaning, on average, samples from this same day in healthy rabbits were used for comparison (results shown in Figure 1). Subsequently we studied if significant changes in bacterial taxa or OTUs between both groups of rabbits could be detected before ERE onset (Additional file 5). Accordingly, both groups were compared applying the Wilcoxon test to the relative abundance of bacterial taxa or OTUs of samples detected the day of weaning and days 1, 3, 5 and 7 after weaning. In the comparison of day 7 after weaning, samples from two ERE rabbits were not included in the analysis, since these samples corresponded to the day of ERE onset in these two particular rabbits. Besides, as a third comparison, samples obtained the morning (8–9 a.m.) of the day of ERE onset (during that particular day a drop in food consumption was detected) were compared to samples obtained from age-matched healthy littermates using the same statistical analysis (Figure 3). On average, samples collected on the day of onset in ERE rabbits correspond to the day 9 after weaning. For this reason, we compared ERE onset samples with samples from healthy rabbits obtained on day 9 after weaning.
For determining if the evolution of the rabbits’ microbiota since the day of weaning until day 7 post-weaning is different in healthy rabbits compared to those that will develop ERE, linear model regression analysis (Figure 2A) was performed using Graphpad.