To collect RNA samples for metatranscriptomic analysis, oral swabs were taken from two adult domestic racing pigeons (Columba livia), from the Northeast of England. Birds sampled during this study were routinely treated with antibiotics including ronidazole (antitrichomonal) within 3 months previous to sampling, and showed no signs of disease. All birds were allowed to fly outside daily. A FLOQswab (Copan) was moistened with nuclease-free water and rolled around one half of the inside of the bird’s mouth before suspending in 0.7 mL RNAlater (ThermoFisher Scientific). Swabs were transported on dry ice and stored at −80°C until RNA extraction (within 2 months). Each bird was also screened for the presence of T. gallinae. An additional sterile cotton-tipped swab was moistened with sterile distilled water and rolled around the opposite side of the bird’s mouth before suspending in a tube with 13 mL TYM medium (Clark and Diamond, 2002) modified with heat-inactivated horse serum instead of bovine serum, pH adjusted to 7.2, supplemented with 87 μM ferric ammonium citrate, 1 unit/ml penicillin and 100 μg/mL each of streptomycin, nystatin and kanamycin. Tubes were sealed in parafilm and incubated at 37°C. An additional swab was moistened in water and exposed to the air for a few seconds before suspending in identical media as a negative control. Cultures were examined daily under a microscope for the presence of motile trophozoites.

Material stored in RNAlater (ThermoFisher Scientific) was thawed on ice, diluted with 0.7 mL nuclease-free PBS and pelleted by centrifugation at 6 k × g for 5 min at 4°C. RNA was extracted from the resulting pellet using TRIzol (ThermoFisher Scientific) according to the manufacturer’s instructions. RNase-free glycogen (10 μg) was used as a carrier during the precipitation phase to improve RNA yield. RNA yield was measured using the Qubit RNA high sensitivity kit (ThermoFisher Scientific), and sample purity and integrity was verified by TapeStation (Agilent) electrophoresis using RNA ScreenTape. To process samples for metantranscriptomic sequencing, total RNA was depleted of rRNA using the QIAseq FastSelect 5S/16S/23S kit (beta version, Qiagen), sequencing libraries were generated using the SMART-Seq Stranded for total RNA-seq kit (Takara) and 75 base pair, paired end reads were generated using an Illumina MiSeq v3 platform.

Metagenomics sample collection, processing, sequencing, and data quality control were performed in the lab of Professor Zheng Wang (Beijing Academy of Agriculture and Forestry Sciences). Oral samples were collected from domestic farm pigeons from Beijing, China. A sterile cotton-tipped swab was rolled across the inside of the crop before suspension in 0.5 mL 0.9% saline. Three swabs were collected from each bird. Swabs were used for diagnosis of T. gallinae infection by wet mount microscopy; smears were examined by light microscopy for motile parasites. To screen samples for the bacteria of interest by PCR, saline samples from T. gallinae-infected birds were centrifuged at 18,630 × g for 5 min, resuspended in 50 μL TE (Tris-EDTA) buffer, heated at 100°C for 5 min and stored at −20°C until required. PCR was performed in 25 μL volumes with 1 × Premix Taq^™ DNA Polymerase (TaKaRa), 400 nM forward and reverse primer and 2 μL Template DNA. Primers targeted the 23S rRNA gene sequence of the target putative Mycoplasma sp. (F: 5’-TAGGACCCGACTAACCCAGA-3′, R: 5’-TTCTGCGCCGAAGATTCAAC-3′). Thermal cycling conditions were performed using a SCILOGEX TC1000-G cyler as follows: initial denaturation at 95°C for 7 min, followed by 35 cycles of denaturing at 94°C for 1 min, annealing at 55°C for 1 min and extension at 72°C for 1 min, with a final extension at 72°C for 10 min. Products were analysed by electrophoresis on a 1.5% agarose gel.

Metagenomics sequencing was performed on samples positive for the bacteria of interest. DNA was extracted from saline-suspended swabs using the E.Z.N.A.^® Soil DNA Kit (Omega Bio-tek) according to the manufacturer’s instructions. DNA concentration was measured using the QuantiFluor dsDNA System (Promega) measuring fluorescence on a TBS-800 fluortimeter. DNA purity and integrity were assessed, respectively, using a NanoDrop2000 spectrophotometer and by AGE using a 1% gel. DNA was fragmented to approximately 400 bp using a Covaris^® M220 ultrasonicator (Gene Company Limited) and used to construct a paired end sequencing library using the NEXTflex^TM Rapid DNA-Seq kit (Bioo Scientific). Approximately 40 million 150 bp paired end reads were generated using a NovaSeq sequencer (Illumina). Fastp (Chen et al., 2018) was used to trim adapater sequences and low quality reads (any reads shorter than 50 bp, with average Phred quality score lower than 20, and with any undetermined bases).

PCR was used to close gaps in the resulting metagenomic scaffold and confirm their circularity as described above. The primer pair F: 5’-ATGATTTGAATTGCCCTTCCAG-3′ and R: 5’-AGTTGTTTCGCGTGGTTATG-3′ was used to close a single gap in the target scaffold at an annealing temperature of 56°C. The primer pair F: 5’-GGGATATTGTTATTCCGGCTAAGG-3′ and R: 5’-TCCAAGAGAAGTTACACTATGAGG-3′ was used to confirm scaffold circularity at an annealing temperature of 58°C. Products were sent for GENEWIZ^® sequencing by Azenta Biotechnology Co. Ltd. Sequencing was performed in both directions.

We re-sequenced the target bacterial genome using a pigeon oral sample collected for metatranscriptomics analysis (derived from a pigeon from Northeast England). Total DNA was isolated from the remaining material from swab by the Chelex method. Material from the swab tip was excised using a scalpel and was rehydrated by incubating at RT in 300 μL sterile distilled water for 30 min. Swab material was pelleted at 16.1 k × g for 3 min, all but approximately 50 μL of the supernatant was removed, and 150 μL 10% Chelex 100 (Bio-Rad) was added. Samples were vortexed for 15 s, centrifuged at 16.1 k x g for 10 s, incubated at 95°C for 20 min, vortexed again for 15 s, and centrifuged at 16.1 k × g for 3 min. The gDNA-containing supernatant was collected and stored at −20°C until required. Whole genome amplification was performed using a REPLI-g kit (Qiagen), before sequencing library preparation using the KAPA hyper kit (Roche). Sequencing was performed on an Illumina NovaSeq 6,000 sequencer.

Quality of raw sequencing data was assessed with fastQC (Andrews, 2018) and adapter sequences were clipped using the fastx-toolkit (Hannon, 2010). Read taxonomy was assigned using Kraken2 (Wood et al., 2019) with default parameters and the NCBI non-redundant nucleotide database as a reference. A de novo transcriptome assembly was generated with rnaSPAdes v3.11.1 (Bushmanova et al., 2019) under default parameters, using only reads classified as “Candidatus Malacoplasma girerdii.” Assembly quality was assessed by mapping “Candidatus Malacoplasma girerdii” reads back to the transcriptome using Bowtie2 (Langmead and Salzberg, 2012), reporting up to 20 alignments per read, and summary statistics were generated using Samtools (Li et al., 2009). Contigs used for downstream analysis were visually inspected as part of an alignment with a range of homologous sequences to check for obvious errors.

Sequence data from 16S rRNA profiling of the pigeon digestive tract was downloaded from the NCBI’s SRA database using accessions SRR9191928-SRR9192035 (Ji et al., 2020). OTUs (operational taxonomic units), clustered at 97% similarity were generated and quantified by read mapping as described by (Ji et al., 2020). The Ribosomal Database Project Classifier (Wang et al., 2007) was used to assign taxonomy to OTUs at the lowest available taxonomic rank with greater than 70% confidence score. BLASTn search (Altschul et al., 1990) was used to identify OTUs similar to query sequences of interest. Assessment of presence or absence of OTUs of interest between samples was defined by requiring at least 0.1% of total reads to be assigned to the OTU. As indicated in the figure legends, differential abundance of OTUs between conditions was assessed using the ANCOM-BC test (Lin and Peddada, 2020) and the Mann–Whitney U test, the latter using abundances normalised as a percentage of the total library size. ANCOM-BC-derived p-values were adjusted to correct for the multiple comparison problem using the Bonferroni method. All statistical analyses were performed in R.

Pigeon oral metagenome data was utilised to generate the genome sequence of a newly discovered Mycoplasma sp. associated with T. gallinae. Quality of sequencing data was confirmed using fastQC (Andrews, 2018). Host-derived reads were filtered from the data by alignment to the C. livia reference genome (NCBI genome ID: 10719) using Bowtie2 (Langmead and Salzberg, 2012) with default parameters. A de novo assembly was generated from the remaining microbial reads using metaSPAdes (Nurk et al., 2017) with default parameters. To assess circularity and coverage of the resulting scaffolds, paired end microbial reads were aligned to the assembly using Bowtie2 (Langmead and Salzberg, 2012) with default parameters. Sambamba (Tarasov et al., 2015) was used to calculate read coverage across scaffolds using a window size of 2000 bp. Artemis genome browser (Carver et al., 2012) was used to examine read alignment to the assembly and assess scaffold circularity. The Seqinr R package (Charif and Lobry, 2007) was used to calculate GC and tetranucleotide frequency across scaffolds. BUSCO (Bacterial Universal Single Copy Orthologs; Manni et al., 2021) was used to assess genome completeness using the Mollicutes reference dataset. Whole genome alignment was performed using Mauve (Darling et al., 2004) with default parameters.

A taxonomic sample of 16S and 23S rRNA genes sequences of interest was obtained by searching the literature and BLASTn search against the NCBI non-redundant nucleotide database (O’leary et al., 2016). Sequences were aligned using Clustal Omega (Sievers et al., 2011) and trimmed using trimAl with a gap threshold of 0.9 (Capella-Gutiérrez et al., 2009). All alignments were visually inspected in SeaView (Gouy et al., 2010). IQ-TREE v2.2.0 (Minh et al., 2020) was used to generate maximum likelihood phylogenies using automatic model selection based on the best fit between the model and data according to the Akaike information criterion (AIC) and Bayesian information criterion (BIC). The selected model is listed in the figure legend for each phylogeny. Branch confidence was assessed by computing 1,000 standard bootstrap replicates. iTOL (Letunic and Bork, 2019) was used to generate annotated figures.

The amino acids of 55 protein coding gene (single copy genes) datasets, with 83 taxa each, were aligned with muscle v3.8.31 and sites with poor alignment were identified with BMGE using default settings (version 1.12) and removed (Supplementary Table 2). Prior the phylogenetic analyses any duplicate sequences were then removed. Phylogenetic analyses of these alignments used IQ-TREE, using ModelFinder to choose the best fitting model, and using 1,000 ultrafast bootstrap pseudoreplicates. Each best fitting model is listed in Supplementary Table 2. Taxa that were removed as duplicates in the previous step were then restored to the trees as multifurcations. Compatibility of the gene trees was assessed before concatenating the alignments for those genes. To check for compatibility, we first made a 95% majority-rule consensus tree. Then for each split in that consensus tree we examined each of the 55 gene trees for incompatible splits. In order to be flagged as incompatible, we required that the split in the gene tree have ultrafast bootstrap support of 95% or more, which flagged 1 of the 55 gene trees as being incompatible (Family #2279at2085, Xylose isomerase-like superfamily; Supplementary Table 2). That gene tree was removed, and a 95% majority-rule consensus tree was made with the remaining 54 gene trees, which identified no incompatible trees. The 54 compatible alignments were concatenated for further analysis, making up an alignment of 13,845 sites (Supplementary Data File 2). Phylogenetic analysis of the multi-gene alignment was performed using IQ-TREE v2.2.0. A guide tree was constructed using the LG4X model (Le et al., 2012) with 1,000 ultrafast bootstrap pseudoreplicates, and the final phylogeny was generated using the LG + C20 + F + G model and the PMSF strategy (Wang et al., 2018) with 1,000 ultrafast bootstrap pseudoreplicates.

Network analysis was used to predict homologous gene families among the annotated proteins from the genomes of a range of Mollicutes sp. using EGN (Halary et al., 2013). All vs. all BLASTp was used as a similarity search, and genes were clustered into networks with edges linking gene nodes according to BLAST hits. Thresholds for significant BLAST hits were an E value of less than 10–5, percentage identity of greater than 20% and an alignment length of at least 25 AAs and 20% of the length of the shortest sequence and 90% of both sequences. Genes with no BLAST hits above this threshold were considered to be singletons. Networks were assessed using the igraph R package version 1.2.11 (Csardi and Nepusz, 2006). Figures were generated using Cytoscape version 3.6.1 (Shannon et al., 2003).

A Malacoplasma sp. related to the T. vaginalis symbiont “Ca. M. girerdii” was serendipitously discovered during screening of domestic pigeon oral samples in Northeast England. Two pigeon samples were examined for the presence of T. gallinae using the culture-based method, one of which was positive for the parasite. One sample for each condition in which T. gallinae was present and absent were subjected to metatranscriptomics analysis by Illumina sequencing. Results and metadata are summarised in Table 1. Classification of reads by k-mer alignment showed that the T. gallinae-positive sample contained a substantial number of reads with similarity to “Ca. M. girerdii,” whereas similar reads in the T. gallinae negative sample were almost absent. Reads assigned as M. hominis were also virtually absent from both samples, consistent with the absence of this organism. For this data, library rRNA depletion was unsuccessful due to the use of a beta-version of the QIAseq FastSelect 5S/16S/23S rRNA depletion kit (Qiagen), precluding accurate quantification of relative taxonomic abundances in these samples.

To further investigate the identity of the putative bacterial species, we generated a de novo assembly using the “Ca. M. girerdii”-like reads. A total of 15 contigs were generated (Supplementary Table 1; Supplementary Data File 1 with the sequences), and a full-length 23S rRNA sequence was identified from the largest contig by alignment with a range of Mycoplasma 23S rRNA sequences. Phylogenetic analysis of the 23S rRNA sequence (Figure 1) revealed that the putative bacteria is a distinct but closely related sister taxon to “Ca. M. girerdii,” with maximum bootstrap support.

To allow sequence comparison with published 16S rRNA sequence profiling we generated a partial 16S rRNA sequence for the bird-associated mycoplasma species. We designed a reverse primer based on the metatranscriptomics 23S rRNA contig used in combination with the generic bacterial 16S rRNA forward primer 338F (Yu et al., 2005; Supplementary Figure 1). We generated a PCR product of the expected size (~1.8 kb), using the metagenomic gDNA template from pigeon 017. The sequences corresponding to the 16S rRNA V3-V4 hypervariable region from 3 independent clones were 100% identical. Phylogenetic analysis demonstrated the pigeon-derived 16S rRNA sequence was closely related to “Ca. M. girerdii,” with 99% bootstrap support for the lineage (Figure 2). Herein we describe this putative mycoplasma species as “Candidatus Malacoplasma trichamica.”

Our preliminary metatransciptomics results were consistent with a putative association between T. gallinae and “Ca. M. trichamica.” To investigate this further, we searched for “Ca. M. trichamica”¬-like 16S sequences among the operational taxonomic units (OTUs) identified in a published 16S ribosomal profiling dataset derived from the crop and gut of 36 farmed pigeons in Beijing, China (Ji et al., 2020). Otu13, showed 100% identity to “Ca. M. trichamica” and 97.2% identity to the “Ca. M. girerdii” 16S rRNA sequence. In contrast, there were no close hits among the OTUs using M. hominis as a query. Consistent with these observations, Otu13 was assigned as Malacoplasma (historically called Mycoplasma) by the Ribosomal Database Project (RDP) classifier (Wang et al., 2007). A phylogenetic analysis confirmed Otu13 to be closely related to the pigeon-derived mycoplasma 16S rRNA sequence from bird 017 with maximum bootstrap support and “Ca. M. girerdii” with 99% support (Figure 2).

Among the cohort of 36 pigeons, 12 each originated from groups heavily infected (high grade; HG), lightly infected (low grade; LG), and uninfected with T. gallinae (CG). For each bird, samples were collected from the crop, ilium and rectum (Ji et al., 2020). We assessed the distribution Otu13 among the cohort, to investigate any potential association between T. gallinae and the putative Mycoplasma sp. (Table 2). Among samples from the crop, Otu13 was found exclusively in birds infected with T. gallinae, and was present in the crop of 37.5% of the T. gallinae-infected birds. Otu13 was found in only three gut-derived samples, two of which were from birds infected with T. gallinae.

Differential abundance (DA) analysis also supported the association (Figure 3). Out of 954 OTUs, 261 were identified to be significantly differentially abundant in the crop of T. gallinae-infected vs. control birds (adjusted Q-value less than 0.05). Among all OTUs, Otu13 showed the 4th greatest significant difference in abundance between the crop of T. gallinae-infected birds and uninfected controls, although there was no significant difference in abundance between the HG and LG conditions. Interestingly, the abundance of Otu13 was relatively high in several samples (the highest sample abundance was 46% of 16SrRNA reads; Figure 4). Otu13 also showed a somewhat mutually exclusive pattern of abundance with some of the other OTUs most strongly associated with T. gallinae, particularly Otu34 (Figure 4). We found no significant difference in abundance of Otu13 between squabs of different ages, however there was tentative evidence that Otu13 was more abundant in the crop than the gut in T. gallinae-infected birds, consistent with its co-localisation with T. gallinae (Figure 2). Intriguingly, among the 5 most significantly increased OTUs in the T. gallinae-infected crop, Otu77 (Figure 2) shares 97% sequence similarity with Pelistega europaea, a species isolated from pigeons with respiratory disease (Vandamme et al., 1998). In addition, our results suggested a positive correlation between T. gallinae presence and two OTUs assigned to the family Prevotellaceae (Out34 and Otu92) and a negative correlation between T. gallinae presence and two OTUs assigned to the Lactobacillus genus (Otu177 and Otu200).

We attempted to generate a whole genome sequence for “Ca. M. trichamica.” A single domestic pigeon was confirmed to be naturally infected with T. gallinae by wet mount microscopy. Infection with the “Ca. M. trichamica” was confirmed by PCR using primers designed based on the metatranscriptome-derived ““Ca. M. trichamica” 23S rRNA sequence. Metagenome sequencing was performed on a single sample derived from the crop of this bird. After de novo assembly of the sequence data, the largest scaffold (scaffold 1) showed characteristics suggesting a near-complete genome. Numerous paired sequence reads spanned between the beginning and end of the scaffold, indicating circularity. A complete metagenomics assembled genome (MAG) sequence was finalised via targeted PCR and Sanger sequencing to confirm scaffold circularity and fill a single ~100 bp sequence gap (Figure 5). The scaffold was ~657 kb, with a GC content of 28.8%, close to the 629 kb, 28.6% GC genome sequence of “Ca. M. girerdii” isolate UC_B3 (NCBI assembly accession ASM221542v1). GC content and tetranucleotide frequency were relatively uniform across the MAG (Figures 6A,B), suggesting that the assembly was not chimeric.

We performed a BUSCO (Bacterial Universal Single Copy Orthologs) analysis of the genome completeness, using a reference dataset of Mollicute genes. The results suggested that 147 out of 151 expected single copy genes were present, with an additional putative gene fragment representing one of the remaining 4 genes (annotated as 50S ribosomal protein L15). For comparison, BUSCO analysis of the “Ca. M. girerdii” UC_B3 genome revealed 148 out of 151 expected single copy genes were present. The list of expected genes which were entirely absent was identical between scaffold 1 and “Ca. M. girerdii.” However, the 50S ribosomal protein L15 existed as a complete gene in “Ca. M. girerdii.” We used the genomic location of the dnaA gene as a position indicator for the replication origin, as the two are typically adjacent (Mackiewicz et al., 2004). The dnaA gene was located by tBLASTn (Altschul et al., 1990) search using the DnaA protein sequence from “Ca. M. girerdii” UC_B3 (GenBank accession: ASJ88910.1). Coverage across the MAG (proportional to DNA abundance) decreased uniformly from ~5,000X at the putative replication origin to approximately half this value at the opposite region of the presumed circular scaffold (Figure 6C). This suggests that the genome was undergoing replication, indicative of active growth. Consistent with this, 45.4% of non-host metagenomic sequence reads aligned to the MAG, suggesting high abundance of “Ca. M. trichamica” in the sampled pigeon.

We re-sequenced an additional isolate of “Ca. M. trichamica” via metagenomics derived from the UK sample pigeon 017, from which the initial metatranscriptomics data was generated (Figure 1). Overall, 0.04% of the UK-derived metagenomics reads aligned to the complete MAG, indicating a low overall abundance. A consensus sequence derived from the UK data aligned with the complete MAG with 95.2% coverage and 91.9% overall percentage identity (Figure 5), indicating a high level of similarity between the two genomes as expected.

To complement the rRNA based phylogenies (Figures 1, 3) we performed a phylogenomic analysis exploiting the “Ca. M. trichamica” genome sequence using a set of 54 single copy orthologous proteins shared among a diverse sample of available Mollicutes genomes (Figure 7). Consistent with phylogenies based on the rRNA sequence, “Ca. M. trichamica” was most closely related to “Ca. M. girerdii” with maximum bootstrap support. We performed a whole genome alignment of the “Ca. M. trichamica” MAG alongside the two available “Ca. M. girerdii” genomes at GenBank (Figure 8). The two “Ca. M. girerdii” genomes shared extensive homologous regions and were largely syntenic, although our results supported a large-scale inversion at the ~170 kb position previously described (Costello et al., 2017). The “Ca. M. trichamica” MAG shared several large homologous regions with “Ca. M. girerdii,” although fewer than between the “Ca. M. girerdii” genomes. There appeared to be some syntenic regions between the “Ca. M. girerdii” and “Ca. M. trichamica” MAGs, but results also suggested extensive genomic rearrangement.

We performed gene family analysis in order to identify potential distinguishing gene features of the “Ca. M. girerdii”-like lineage of Trichomonas-symbiotic Mollicutes. We included the annotated coding sequences from the “Ca. M. trichamica” MAG alongside a lineage of related Mollicutes sequences as defined by the phylogenomic analysis (Figure 8). We also included M. hominis due to its relevance as a symbiont of T. vaginalis (Dessì et al., 2019). The analysis included 24 Mollicute genomes. Overall, there was a low CDS content; among all genomes analysed, “Ca. M. trichamica,” “Ca. M. girerdii” and M. hominis had the 4th, 6th and 7th fewest total CDS, respectively (Table 3). However, this did not correspond to a reduced protein family repertoire as “Ca. M. trichamica,” “Ca. M. girerdii” and M. hominis ranked midrange for combined gene families and singletons (17th, 14th and 13th).

A general feature specific to “Ca. M. trichamica” and “Ca. M. girerdii” compared with related Mollicutes appeared to be an adaptation towards an anaerobic lifestyle. Related gene families unique to “Ca. M. trichamica” and “Ca. M. girerdii” within our analysis included Pyruvate:ferredoxin oxidoreductase, rubrerythrin, desulfoferrodoxin, ferrous iron transport protein B and FeoA domain-containing protein (Chabrière et al., 1999; Lombard et al., 2000; Sztukowska et al., 2002; Lau et al., 2016). Additionally, amino acid catabolism was shared feature, indicated by alanine dehydrogenase and L-serine ammonia-lyase gene families also unique to “Ca. M. trichamica” and “Ca. M. girerdii” within our analysis. Notable other gene families apparently unique to “Ca. M. girerdii” and “Ca. M. trichamica” included 2′,3′-cyclic-nucleotide 2′-phosphodiesterase, and anaerobic ribonucleoside-triphosphate reductase (ARTR). A single contig (Node 11) assembled from the pigeon metatranscriptomics data shared 96.8% sequence identity at the nucleotide level with the putative ARTR encoded by the “Ca. M. trichamica” MAG across its full length. Apparently unique to “Ca. M. girerdii” in our analysis was a putative class IIb bacteriocin (Genbank accession: ASJ89026.1).

BspAs, a diverse group of cell surface proteins mediating cell–cell interactions (Sharma, 2010), are of interest due to their annotation in the genome of both “Ca. M. girerdii” and T. vaginalis. In our analysis, BspAs formed a large protein family with 159 members distributed across 11 Mycoplasma spp. genomes. These genomes appear to originate from three separate linages as highlighted in Figure 7. Due to the repetitive nature of BspA proteins (characterised by leucine-rich repeats) phylogenetic analyses are problematic and so we employed gene network analysis to investigate their relationship. Graphical network analysis of the BspA family based on sequence similarity suggested two major distinct groups, one of which was mostly confined to “Ca. M. girerdii,” “Ca. M. trichamica” and M. hominis but also included a single member encoded by the genome of Mycoplasma pirum (Figure 9). BspA sequences from “Ca. M. girerdii,” M. hominis and “Ca. M. trichamica” were absent from the other group, which was associated with the first group via a single BLASTp hit (33.9% sequence identity) with M. pirum. Clustering of putative M. hominis BspAs was apparent from the network, but those from “Ca. M. girerdii” and “Ca. M. trichamica” were interspersed.