Sampling the anal gland and perianal region of domestic cats. Companion cats (N = 23) that presented to the Veterinary Medical Teaching Hospital at UC Davis between December 2021 and March 2022 for elective procedures were evaluated for their suitability for enrollment in the study. The elective procedures included dental cleaning, abdominal ultrasounds, radiographs, or oral examinations that required sedation or general anesthesia. Sedation or general anesthesia of companion cats at the UCD Veterinary Medical Teaching Hospital is carried by anesthesiologists and clinicians using established protocols (e.g. Dexmedetomidine at 4mcg/kg IV and Methadone at 0.2mg/kg IV for sedation, or Maropitant at 1 mg/kg IV and Butorphanol at 0.2–0.4 mg/kg IM for anesthesia). After obtaining informed written consent from cat owners, a board-certified internist and gastroenterologist (S.L.M) manually expressed the anal glands of cats while they were sedated by inserting a lubricated gloved index finger into the cat’s anus and digitally squeezing the anal gland between the index finger and thumb into a sterile 2” x 2” gauze sponge. Collected anal gland material was immediately transferred onto three sterile Puritan cotton swabs (one for microbiome analysis, one for microbial culturing, and one for metabolomics). In addition, a swab from the cat’s perianal region was collected for comparison to the microbiome found in the anal gland. Finally, an unused Puritan cotton swab (blank negative control) and a swab of the internist’s examination glove prior to the procedure were collected as well to assess background contamination.
The swabs were placed in 2 mL screw cap tubes (for microbiome or culturing analysis) or 20 mL borosilicate glass vials (for metabolomic analysis) and stored at -80°C until laboratory and chemical analysis. After sample collection, the cat’s perianal region was cleaned with warm water and sprayed with a feline deodorizer. The cat was closely monitored and observed during recovery. Data on the animal’s signalment (age, breed, sex), health status (including any medical diagnoses), lifestyle (living environment), body weight, body condition, and diet was also collected (Table S1). At the time of sample collection, none of the cats had diarrhea and although some had been previously diagnosed with periodontal disease, chronic enteropathy, renal disease, or intranasal disease, they were asymptomatic for these conditions and deemed overall healthy. With the exception of a single cat, individuals had not been treated with systemic antibiotics within 6-months of enrollment in the study. We kept this cat in our dataset as its microbiome was not compositionally anomalous in any obvious way compared to those from other surveyed cats.
The study was approved by the University of California, Davis, Institutional Animal Care and Use Committee (IACUC protocol # 22528). All methods were performed in accordance with the relevant guidelines and regulations, including the ARRIVE guidelines.
Bacterial culturing and isolation from the anal gland. Bacterial cultures of the anal gland were obtained from 11 of the 23 cats (Table S1); we did not culture microbes from the remaining 12 cats due to timing. Bacterial swabs from the anal gland were vortexed with 1 mL of PBS and two serial 1:10 dilutions were performed. For each mixture, 150 µL was pipetted into lysogeny broth (LB), brain heart infusion (BHI), and blood agar (BA) plates. The plates were placed in BD GasPak EZ Anaerobic System boxes (BD Biosciences, NJ, USA) with packets of CO2 generators to maintain an anaerobic environment. After growing plates at 37°C for 3–5 days, colonies with distinct morphologies, colors, and textures were picked from any or all of the three types of agar plates and plated until pure colonies were obtained. Single colonies were then added to 5 mL of sterile LB, BHI, or BA liquid broth, degassed with nitrogen to remove any oxygen, and grown for three days at 37°C in preparation for DNA extractions.
DNA extraction and Sanger sequencing of bacterial cultures. DNA was extracted from liquid bacterial broth cultures using the Wizard SV Genomic DNA Purification Kits (Promega, USA), according to their protocol for Gram-negative bacteria. Briefly, 1.5 mL of vortexed broth was centrifuged and subsequently exposed to nuclei lysis, RNAse digestion, incubation at 37°C, and protein precipitation on ice. The DNA was washed with isopropanol and ethanol, and the pellet was air-dried and then resuspended in nuclease-free water.
For each bacterial isolate, the 16S rRNA gene was amplified using the 27F (5′-AGAGTTTGATCMTGGCTCAG-3′) and 1391R (5′-GACGGGCGGTGTGTRCA-3′) bacterial-specific primers. The PCR conditions were as follows: an initial denaturation step at 95°C for 3 mins, followed by 30 cycles of 95°C for 45 seconds, 50°C for 60 seconds and 72°C for 90 seconds. A final extension occurred at 72°C for 10 minutes, and a final hold at 15°C. PCR products were purified with the NucleoSpin Gel and PCR Clean-Up kit (Takara Bio, CA, USA) and quantified with Qubit HS dsDNA assay (Thermo Scientific, MA, USA). PCR products were diluted with nuclease-free water until achieving a concentration of ~ 32 ng/µL in preparation for Sanger sequencing. Samples with initial concentrations between 10–35 ng/µL were not diluted, and those with concentrations < 10 ng/µL were not sequenced.
A total of 111 bacterial isolates were submitted for Sanger sequencing of the 16S rRNA gene (27F primer) at the UC Davis College of Biological Sciences DNA Sequencing Facility (Davis, CA, USA). Sanger chromatograms were uploaded to myRDP [54] for quality-trimming and base-calling. The trimmed sequences were searched against the bacterial NCBI RefSeq genomes database [55] using blastn for taxonomic identification setting default parameters. The top hit with the highest e-value and percent identity was selected as that organism’s taxon label (Table S18). Six of the 111 isolates did not meet sequence quality thresholds and could not be classified taxonomically.
DNA extraction and metagenomic sequencing of bacterial swabs. Genomic DNA was extracted directly from swabs of the anal gland (N = 23), perianal region (N = 6), and controls (a sterile swab and swab of glove; N = 2) using the QIAGEN DNeasy Powersoil Pro Kits (Qiagen, MD, USA) (Table S1). Swabs were incubated with the QIAGEN CD1 lysis buffer (800 mM guanidine hydrochloride; 30 mM Tris•Cl, pH 8.0; 30 mM EDTA, pH 8.0; 5% Tween 20; 0.5% Triton X-100) at 65°C for 10 minutes and then underwent bead-beating for 1.5 mins before resuming the manufacturer’s protocol for this kit. Genomic DNA from 31 samples was treated with RNAse A and sequenced on an Illumina NextSeq 500 (300 cycles; Illumina, CA, USA) at the UC Davis Genomics Core to generate PE x 150 bp reads.
The general protocol used by the UC Davis Genome Core for this work was as follows: briefly, barcode-indexed sequencing libraries were generated from genomic DNA samples sheared on an E220 Focused Ultrasonicator (Covaris, Woburn, MA). For each sample 1-130 ng sheared DNA was converted to sequencing libraries using a Kapa High Throughput Library Preparation Kit (Kapa Biosystems-Roche, Basel, Switzerland). The libraries were amplified with 5–15 PCR cycles and analyzed with a Bioanalyzer 2100 instrument (Agilent, Santa Clara, CA), quantified by fluorometry on a Qubit instrument (LifeTechnologies, Carlsbad, CA), and combined in one pool at equimolar ratios. The pools were quantified by qPCR with a Kapa Library Quant kit (Kapa Biosystems-Roche) and subsequently sequenced.
Sequence processing of metagenomic data. Sequenced samples had an average of 5,527,421 (± 997,985) metagenomic paired-end reads which were quality-filtered and trimmed using Trimmomatic v.0.38 setting default parameters [56]. On average, samples retained 95% (± 1.07%) of their sequences after quality-filtering (Table S2). Reads were filtered of host DNA by aligning them to two Felis catus reference genomes (GenBank accessions GCA_000181335.5 & GCA_013340865.1) using Bowtie2 [57]. Fifteen samples retained > 92% of their quality-filtered sequences after removal of host DNA, another four samples retained 62–75% of their sequences, and six samples retained < 50% of their sequences (Table S2). For samples from the perianal region, three retained 50–70% of reads after host DNA filtering, while the other three only retained 12–15% of their reads. Kraken 2 [58] assigned taxonomic classification to the host-filtered quality-trimmed reads for each sample and Bracken (Bayesian Reestimation of Abundance with Kraken) [59] estimated the taxon abundances at the Family, Genus, and species level. On average, 55% of reads were classified by Kraken2 for each sample (range: 26–86%) (Table S2). When Kraken was not able to assign a species label to a sequence, it used the genus or family label followed by spp. (e.g. Bacteroides spp.).
The resulting forward and reverse metagenomic reads were interleaved using a python script from the Ray assembler v.2.3.1 [60] and concatenated into a single file in preparation for metagenome assembly. Sequences were assembled into contigs with (meta)SPAdes v.3.14.1 [61] and the quality of assembly was evaluated with QUAST v.5.0.0 [62] (Table S3). After removing contigs shorter than 300bp using BBMap v.38.87, contigs were uploaded into Anvi’o v.7.1 [63] for gene prediction and functional annotation. Gene prediction was accomplished with Prodigal v.2.6.3 [64], and functional annotation was done using the Clusters of orthologous groups (COGs) [65], and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases [66]. To determine the abundances of genes in each sample, host-filtered reads were mapped to predicted genes using Salmon (v1.8.0) [67]. Gene relative abundances were in units of Transcripts Per Million (TPM), which normalizes for both gene length and sample sequencing depth. On average, 70.35% of host-filtered reads mapped to putative ORFs (range: 46–79%) (Table S2). Lesser percentages were assigned KEGG or COG annotation (~ 46% of mapped reads).
Contigs with a minimum length of 600 bp were binned into metagenome-assembled genomes (MAGs) using MetaBat2 v.2.15 [68]. A total of 85 high-quality MAGs were recovered with completeness scores > 80% and contamination scores < 5% as assessed by CheckM (v1.1.3) [69] (Table S4). The Genome Taxonomy Database Toolkit (GTDB-Tk) (v1.5.0) [70] was used to assign taxonomic identity to MAGs using database release 202 [71]. The abundance of each MAG in a sample was estimated using CoverM (v0.6.1) (https://github.com/wwood/CoverM) by mapping interleaved host-filtered reads to each MAG (Table S5). On average, 60.37% of reads in each sample were able to be mapped to MAGs (range: 32–88%) (Table S2). A phylogeny of these MAGs was constructed with RAxML [72] using the multiple-sequence alignments generated by GTDB-Tk. We had no outgroup and instead rooted our tree to the only member of the phylum Synergistota (Fretibacterium).
Extraction and analysis of anal gland metabolites. Volatile compounds were extracted from anal swab samples using two techniques. First, solid phase microextraction (SPME) fibers (50/30 µm DVB/CAR/PDMS coating) were exposed to swabs in vials to extract VOCs from the headspace. Then, VOCs were extracted from swabs using liquid phase extraction with methanol, followed by a derivatization process before being directly analyzed.
For headspace technique, a 1 µL aliquot of 10 mL/L decane-d22 was added to the 20 mL borosilicate glass vials containing the swabs as an internal standard. Two previously conditioned SPME fibers were exposed to anal swabs for 24 hours at room temperature, then capped and placed in a -20°C freezer until spectrometric analysis. For the liquid phase extraction, these same swabs were placed into 20 mL of methanol for 24 hours at room temperature to extract VOCs into solution. A 2 mL aliquot of each extract was transferred into a new vial and completely dried under nitrogen. Dried extracts were subsequently derivatized by adding 50 µL MTBSFTA (N-tert-Butyldimethylsilyl-N-methyltrifluoroacetamide with 1% tert-Butyldimethylchlorosilane) and 50 µL acetonitrile. Reconstituted samples were left to react for 1 hour at 60°C and stored at -20°C until spectrometric analysis.
VOCs from SPME fibers and derivatized extracts were analyzed with gas chromatography-mass spectrometry. For SPME fibers, the fiber was inserted into the inlet of an Agilent 6890N gas chromatograph (Agilent Technologies Inc.) set to 260°C. VOCs were desorbed from the fiber for 5 min in splitless mode while the GC oven was held at 40°C. For liquid extracts, 1 µL of each sample was injected into the GC inlet held at 260°C. For both samples, the oven was ramped to 120°C at 5°C/min, and then ramped to 280°C at 15°C/min, holding for 10 minutes. VOCs were separated on a DB-5ms column (30 m × 250 µm × 0.25 µm, Agilent Technologies Inc.) with a 1 mL/min constant flow of helium. Compounds were eluted through a 300°C transfer line into an Agilent 5795C mass spectrometer, which scanned 50 to 500 m/z with its source set to 230°C and quadrupole to 150°C.
Samples and blanks/controls were injected in a random order to produce reliable data. A standard Grob mixture was injected in triplicate to monitor instrument performance, and a standard mix of C8-C30 alkanes was analyzed to calculate the Kovats Retention Indices of each VOC. Raw data were first checked for qualitative reasons using Agilent’s Mass Hunter Qualitative Analysis B.06.00 software. GC-MS data were then deconvoluted and aligned using the recursive feature extraction on Profinder (Version B.08.00, Agilent Technologies Inc.) and Mass Profiler Professional (MPP, V13.0). An initial table was then obtained containing all peak intensities (rows) from each sample (column). Peaks from contaminants like siloxanes (base peaks 207, 221 and 281 m/z). Features that appeared in blanks with a signal more than five times the signal from samples (peak sample/blank ratio) were removed. These blanks were composed of: system blanks (instrumental blank without injection), Twister® blanks (injection of clean twisters), blank vials and blank cotton swabs; the latter two were treated as if they were biological samples and underwent SPME extraction and liquid extraction as described above.
Compounds were tentatively identified by matching the mass spectra with structures available in the National Institute of Standards and Technology (NIST) 2020 Library and by matching calculated retention times with those reported in the literature.
Statistical analysis of microbiome and metabolome data. Unless otherwise stated, all sequence data was analyzed and visualized using the R statistical software program (v3.6.2) [73, 74]. Prior to any statistical analysis, we used R decontam to identify and remove contaminant bacterial taxa from the Kraken2 dataset based on their prevalence in control samples compared to samples from the anal gland and perianal region. A total of 58 bacterial species (Table S6) were deemed contaminants by decontam (had scores below the specified threshold of 0.5) and were removed. Read counts assigned to “Homo” were also removed.
The composition of bacterial communities from the felid anal gland (or the perianal region) was visualized using data generated by Kraken2/Bracken (Table S7, Table S8). For this, bar plots showing the relative abundances of bacterial families and genera were constructed using ggplot2 (v3.3.6) [75]. Next, we examined whether five host factors of interest: age (in yrs), obesity (obese vs. not obese), living environment (indoor vs. indoor & outdoor), diet (dry food only vs. other diets), and a medical condition (moderate to severe periodontal disease vs. no disease) could account for any of the variance in anal gland microbiome beta-diversity. For analyses, cats with body condition scores (BCS) 8–9 were considered “obese” and all other cats were classified as “not obese” [76]. This classification was selected because of its clinical relevance and the multitude of obesity-associated comorbidities.
Genus-level abundance data generated by Kraken2/Bracken was converted to presence/absence (for Jaccard dissimilarity), proportions (for Bray-Curtis dissimilarity), or applied a Center-Log-Ratio (clr) transformation (for Aitchinson dissimilarity). Permutational multivariate analyses of variance (PERMANOVAs) tested whether microbiome beta-diversity varied with the five host factors. Tests evaluated all factors simultaneously, in a way where the order of terms did not influence statistical output (e.g. they were marginal PERMANOVA tests). PERMANOVAs employed 999 permutations, set alpha to 0.05, and were conducted with the vegan package (v2.6-2) [77].
We investigated whether microbiome profiles from the anal gland were significantly correlated with metabolome profiles (from solid-phase and liquid-derivatization extractions). For this, the matrix of metabolite absolute abundances (Table S9 for solid extraction, Table S10 for liquid extraction) were normalized by converting to relative abundances (e.g. proportions), log10-transformed to minimize the influence of heteroscedasticity, and scaled with pareto scaling, which is advised for metabolite data [78]. Metabolite Jaccard and Euclidean distances were estimated using the phyloseq package. Mantel tests correlated microbiome matrices (Jaccard, Bray-Curtis and Aitchison) to metabolite Jaccard or Euclidean distances using 999 permutations. Because a statistically significant relationship existed between the two datasets, we plotted this relationship; specifically, we extracted the first two principal coordinates from metabolite and microbiome dissimilarity matrices using the cmdscale function from the stats package [73, 74]. Scatterplots showcased the relationship visually.
We also tested whether the relative abundances of bacterial species were significantly associated with the relative abundances of specific metabolites using Spearman correlations (R stats package). P-values were adjusted for multiple comparisons using FDR. Only bacterial species found in at least 90% of anal gland samples and at a mean relative abundance > 0.2% were considered. For metabolites, all 37 putatively identified metabolites were considered, along with 187 unidentified metabolites that had normalized and scaled mean relative abundances > 0.
Analyses on COG and KEGG ortholog and pathway abundances were also conducted (Table S11-S14). Gene abundances which were estimated in TPM were normalized using total sum scaling. PERMANOVAs examined whether microbiome functions were significantly predicted by host factors, using the methods described above. Mantel tests evaluated whether microbiome functional profiles were correlated with metabolite abundance data using Mantel tests as described for the microbiome-metabolite correlations.
All of the methods described previously in this section were also used to analyze metagenome-assembled genomes (MAG) abundances. PERMANOVAs tested whether MAG relative abundances (Table S5) or their presence / absence were significantly associated with host age, obesity, living environment, diet, and a medical diagnosis of periodontitis. Mantel tests evaluated whether MAG profiles were associated with metabolite profiles. Spearman correlations were used to correlate the abundances of all MAGs classified to species level (47 MAGs) to the abundances of 37 putatively identified metabolites and 193 unidentified metabolites.
Lastly, marginal PERMANOVA tests were also employed to test whether the taxonomic and functional compositions of microbiomes in the anal gland were distinct from those in the perianal region.