How many microbial species are there in human tumor and normal tissues?

Introduction

Recent reviews have proposed that polymorphic microbial communities, also known as heterogeneous microbiomes, could represent an emerging cancer hallmark or enabling characteristic. Cancer hallmarks refer specifically to core biological capabilities that drive tumor development and progression. In contrast, enabling factors facilitate the acquisition of hallmark abilities without being hallmarks themselves. In his 2022 synthesis review, Hanahan (2022) suggested phenotypic flexibility and disrupted development as potential additional hallmark capabilities, adding to the currently well recognized eight hallmarks that he co-proposed about a decade ago. He also posited that non-mutational epigenetic changes and polymorphic microbiomes act as distinctive enabling characteristics that help tumors develop hallmark properties. Similarly, Lythgoe et al (2022) directly referred to microbes as an emerging hallmark of cancer in their 2022 review. Hanahan’s distinction of microbiomes as enabling factors rather than hallmarks reflects the slightly different level of the recognition of the critical role of microbes in cancer research. Regardless of the ultimate position on the role of microbiomes in cancer development and progression, their critical importance has been well recognized in the recent decade. For example, microbiome-immune cell interactions likely influence a tumor’s ability to evade immune destruction, one of the currently established key cancer hallmarks. Heterogeneous microbiomes can profoundly impact cancer phenotypes by differentially affecting hallmark processes. Obviously, the microbiome also contributes heterogeneity between patients (Sepich-Poore et al. 2021).

Sepich-Poore et al. (2021) discussed the role of microbes in cancer from a historical and modern perspective. Early studies dating back 4000 years linked cancer to microbes (Sepich-Poore et al. 2021). One of the first clinical studies in 1868 observed tumor regressions in streptococcus-infected patients, providing an early demonstration of immunotherapy. However, these claims faced reproducibility issues and toxicity concerns over the next century. While the viral theory of Rous sarcoma virus gained traction in 1911, decades of searching failed to find viruses causing human cancers. Instead, somatic mutations are now primarily linked to many cancers. Recent studies reconsider the importance of bacteria and fungi in cancer and immunotherapy through immune-mediated mechanisms (Sepich-Poore et al. 2021). The human microbiome project enabled large-scale metagenomic sequencing, fueling renewed interest in microbe-cancer relationships over the past decade. Landmark studies in 2020 (Nejman et al. 2020, Poore et al. 2020) provided experimental and computational evidence, signifying a breakthrough in understanding the cancer microbiome’s role and potential for improving immunotherapy.

Of the estimated 10¹² microbial species on Earth, only 11 have been identified as human carcinogens (Sepich-Poore et al. 2021). Although only a small number of microbes directly cause cancer, many appear to promote tumor development and progression through immune-mediated interactions, a process known as the immuno-oncology-microbiome axis. Key questions in this field include the roles of microbes - whether they are causal, complicit, or merely passive bystanders, and the current understanding of intertumoral microbes (Sepich-Poore et al. 2021). This article seeks to tackle one outstanding issue, part of the open problem regarding intertumoral or tissue microbiomes, seemingly straightforward as indicated by the title: the diversity of microbial species found in human tumor and normal tissues. Our primary focus is on the overall diversity of microbes in human tissues, a topic that, to our knowledge, has not been previously explored. While comprehensive studies have been conducted on local microbial diversity, known as alpha-diversity, and even beta-diversity, which evaluates regional or inter-locality diversity (including our own studies, Ma 2024a, b), the global diversity of microbes in human tissues remains an uncharted territory.

The challenge of estimating total or global biodiversity was first posed in the 19th century (Watson 1835), during the era of Charles Darwin when naturalists were passionately cataloguing the flora and fauna of our planet. A straightforward approach might be to aggregate the data collected by naturalists worldwide. However, two significant issues arise with this approach. Firstly, the vastness of the Earth makes it impossible for naturalists to reach every location where organisms exist. Secondly, there may be overlapping data in the catalogues compiled by different naturalists, whether from the same or different regions. A potential solution to the first issue could be implementing effective sampling schemes, as it is neither feasible nor necessary to conduct exhaustive counts of species numbers, especially considering the ongoing processes of speciation and extinction. The second issue could be addressed through an automated algorithm designed to eliminate overlaps. The complexity of these issues necessitates the use of statistical or stochastic algorithms to handle the challenges of sampling, stochasticity, and overlap (e.g., Connor and McCoy 1979, Chao et al. 2014).

The first attempt to estimate biodiversity was made by British plant biogeographer and evolutionist, Hewett Cottrell Watson (1804–1881), who introduced the concept known as the species-area relationship (SAR) (Watson 1835). The SAR correlates the accumulation of species numbers (S) with the size of the area (A) where the species are found, using a power function (S = cA^z). Intuitively, a larger area should be inhabited by more species, or in other words, the larger the area a naturalist surveys or samples, the more species they should collect. For this reason, SAR is also referred to as the collector’s (or naturalist’s) accumulation curve. The SAR can be transformed into a linear relationship on a logarithmic scale, i.e., ln(S) = In(c) + zln(A), which makes the increasing correlation between species number and area size even more apparent. Since Watson’s pioneering work, the SAR has been extensively and intensively studied in community ecology and biogeography (e.g., Preston 1962, Connor and McCoy 1979, Rosenzweig 1995, Plotkin et al. 2000, Ulrich & Buszko 2003, Tjørve & Tjørve 2008, Triantis et al. 2012). For instance, the SAR served as a foundational model for MacArthur and Wilson’s (1967) island biogeography theory, which significantly influenced community ecology in the 1960s and 1970s. In practical terms, the SAR has arguably become the most critical model in conservation biology, particularly for protecting biodiversity and endangered species, influencing decisions such as the appropriate size of conservation zones for endangered species (Rosenzweig 1995).

Despite its widespread applications, the classic SAR model presents two issues. The first pertains to the ever-increasing or decreasing nature of the power function in the SAR model, which lacks “saturation” points or extreme (maximum or minimum) values. Given the Earth’s finite size, it’s reasonable to argue that the number of species should also be finite. To address this issue, Plotkin et al. (2000) and Ulrich & Buszko (2003) introduced the power law model with exponential cutoff (PLEC) and power law with inverse exponential cutoff (PLIEC), both of which incorporate saturation points (extreme values). With PLEC or PLIEC, the SAR curve can reach saturation or maximums, corresponding to the maximum number of species on Earth, or alternatively, to the maximum number of species in a specific region, such as human tissues.

The second issue associated with the classic SAR relates to the definition of biodiversity, which, in its simplest form, is the number of species in a region (area), known as species richness (S or R). The problem with this simplified definition is that it fails to account for the fact that not all species are equal: some are abundant (such as ants and many insects), while others are not only rare (e.g., pandas and tigers), but also potentially more valuable. In the context of this study, as previously mentioned, only 11 microbial species have been identified as human carcinogens or oncomicrobes according to the International Association for Cancer Registries (IACR) (cited in Sepich-Poore et al. 2021). It is clear that different microbial species in human tissues have varying oncological significances, and our interest should lie in estimating the numbers of microbial species at various levels of significance.

Indeed, there are numerous diversity metrics (indexes) beyond the simplest measure of species richness. So many, in fact, that practitioners often find themselves overwhelmed by the multitude of choices, lacking a solid standard to guide their proper applications (Magurran 2013, Henderson 2021). Among these many diversity metrics, Shannon’s entropy and Simpson’s index are two of the most commonly used. It might seem that a simple substitution of species richness in the classic SAR with these metrics would solve the problem, but the solution is not that straightforward. Diversity metrics such as Shannon entropy and Simpson’s index are not only incompatible with each other, but they also do not scale in terms of simple mathematical functions like power-law models. This complexity may explain why a breakthrough in this area has been elusive for a long time.

One diversity metric that possesses such appropriate properties is the so-called Hill numbers, first introduced as biodiversity metrics by Hill (1973) from economics. However, it did not garner the attention it deserved among ecologists until its rediscovery by Chao et al. (2014), possibly due to the somewhat abstruse interpretations Hill (1973) used to explain its central concept of “numbers equivalent of elements” in economics. The accomplishment of Hill numbers is actually similar to linking the US dollar to gold at the rate of $35 per ounce under the Bretton Woods system, as per Ma & Li (2024). Hill numbers are now considered the most suitable system for biodiversity metrics, unifying Shannon, Simpson, and other diversity indexes. Against this backdrop, Ma (2018a, 2019) extended the classic SAR into the Diversity-Area Relationship (DAR) using Hill numbers, and also to the Diversity-Time Relationship and Diversity-Time-Area Relationship (DTAR). The extensions also incorporated the adoption of PLEC and PLIEC as DAR models and derivations of maximal accrual diversity (MAD). The MAD or D_max essentially represents the potential diversity or ‘dark’ diversity, accounting for species that may be locally absent but exist in the regional species pool (and may therefore immigrate at a certain time). In the context of this study, it means that MAD can account for microbial species that are absent in some individuals while present in others.

With the DAR-PL and DAR-PLEC models, and leveraging a recent breakthrough in computational approaches to distill large datasets of tumor tissue microbiomes, we can address the question raised in the article title. This has been made possible thanks to the revolutionary AI-machine learning approach by Poore et al. (2020), who produced a substantial dataset of tissue microbiomes from the TCGA (The Cancer Genome Atlas) database. To enhance the robustness of our estimations, we employ extensive permutation tests to manage the inherent stochasticity (uncertainty). Figure 1 and Table 1 provide a roadmap of our study.

• Download figure
• Open in new tab

Fig 1.

A diagram for illustrating the study design and its relationships with previous foundations works (also see Table 1 for the supplementary interpretations on the design). Legends for the equations: ^qD is the diversity in Hill numbers at diversity order q=0, 1, 2, 3; A is the ‘area’ size (e.g., the number of samples or individuals); ^qD_max is the MAD (maximal accrual diversity) of diversity order q, the so-termed potential or ‘dark’ diversity, which includes both local species and species that are absent locally but present in regional species pool.

Material and Methods

Results

Pan-tumor Microbiome Diversity-Scaling Analysis

Given the nearly universal invariance in major DAR scaling parameters, especially z, it is justified to combine all cancer disease types, and even tissue types to analyze diversity scaling across cancer/tissue types. Table S7 exhibits the results of fitting the DAR-PL and DAR-PLEC models with all cancer types and/or tissue types combined, based on the designs of Scheme-III and IV in Table 1. Table 2 below excepts the key DAR parameters for convenience of illustration, with the combined tissue types and cancer types. We further compared the DAR parameters with permutation tests based on 1,000 repetitions of re-sampling, and the test results were exhibited in Table S8. It turned out that no statistically significant differences were detected in any of the comparisons (P<0.05) between different tissue types in Table S8. That is, there are no statistically significant differences between tumor tissue types (PT, SN, or B) in any of the DAR parameters for the models built based on Scheme-III and IV. In other words, on a pan-tumor basis, diversity scaling makes no difference between tissue types. This suggests that tumors do not significantly influence the total microbiomes on a pan-tumor basis or across cancer types.

View this table:

• View inline
• View popup
• Download powerpoint

Table 2. The DAR modeling parameters from the combined datasets of all tissue types and cancer types, for each of the three taxa (Archaea, Bacteria, and Viruses) under each sequencing protocol (RNA-Seq or WGS), excerpted from Table S7.

From Table S7 and Table 2, some very interesting numbers that directly answer the question raised in the article title are worthy of emphasis here. Since no statistically significant differences were detected between tissue types, there was no need to distinguish between PT, SN and B samples to estimate the numbers of microbial species (microbial species richness) or MAD (maximum accrual diversity) in terms of the Hill numbers or species equivalents. Due to the nature of MAD, which is estimated from the saturation point of the DAR-PLEC curve (i.e., the diversity accumulation in terms of Hill numbers), the title question can be converted into how many species there are in human tissue microbiomes. We can depend on the total tissues (PT, SN, B) samples combined to estimate MAD. As to the sequencing protocols, due to potential differences in processing samples and computational pipelines, we do not combine their samples (results). Instead, we treat their estimation results as the range of MAD estimation.

Figure 2 illustrates the potential microbial diversity (Dmax) of pan-tumor microbiome for each taxon (archaea, bacteria, or viruses) and each tissue type (primary tumor [PT], solid normal [SN], or blood [B]), as well as the ‘Total’ of combined tissue types and taxa, at different diversity orders (q=0, 1, 2, 3), under each sequencing protocol (WGS or RNA-Seq). The pan-tumor microbiome concept means that the DAR model is built across cancer types—the microbiomes of all cancer types are combined as the microbial landscape. Therefore, Fig 2 actually shows the potential diversity at the largest pan-tumor scale of tissue microbiomes across all cancer types, including the ‘Total’ that also combined the tissue types besides cancer types. Given the lack of statistically significant differences between tissue types and the saturation nature of D_max estimation, Figure 2 also illustrates the potential diversity of human tissue microbiomes (the last bar in purple).

• Download figure
• Open in new tab

Fig 2.

The potential microbial diversity (D_max) of pan-tumor microbiome for each taxon (archaea, bacteria, or viruses) and each tissue type (PT, SN, or B), as well as the ‘Total’ of combined tissue types and taxa, at different diversity order (q=0, 1, 2, 3), under each sequencing protocol (WGS or RNA-Seq). The pan-tumor microbiome concept means that the DAR model is built across cancer types—the microbiomes of all cancer types are combined as the microbial landscape and the DAR models were built on the landscape (pan-tumor scale). Given the lack of statistically significant differences between tissue types and the saturation nature of D_max estimation, the last bar (in purple) represents the potential diversity of the human tissue microbiomes.

After the previous rational simplification, from Table S7 we conclude that the species richness (q=0) or the total number of tissue microbial species is between 1946 (estimated from RNA samples) to 1948 (estimated from WGS samples), which are rather close to each other and demonstrate the robustness of the Poore et al. (2020) TCGA datasets used in this study. The categorical breakdown of species richness for taxa are archaea between 102 (WGS) and 103 (RNA), bacteria between 1528 (RNA-Seq) and 1536 (WGS), and viruses between 330 (RNA-Seq) and 368 (WGS). Except for viruses, the numbers from both sequencing protocols are rather close and their differences are negligible.

For diversity order q=1, or the number of species equivalents of typical abundances, the numbers are between 21 (WGS) to 36 (RNA-Seq) or 26 on average across protocols and tissue types. The taxa breakdowns of typical species are archaea between 2.7 (WGS) and 6.5 (RNA-Seq), bacteria 21 (WGS) and 37 (RNA-Seq), and viruses between 20 (RNA-Seq) and 26 (WGS).

For diversity order q=2, or the number of species equivalents of dominant abundances, the numbers are between 21 (WGS) and 22 (RNA-Seq). The taxa breakdowns of dominant species are archaea between 2 (WGS) and 4 (RNA-Seq), bacteria between 9 (WGS) and 22 (RNA-Seq), and viruses between 9 (RNA-Seq) and 11 (WGS).

Conclusions and Discussion

Regarding the diversity-scaling analysis of the tumor microbiomes, the previous findings can be summarized in the following main conclusions. Except for D_max, the diversity scaling parameters are generally not significantly different between tissue types and among cancer types (except for q=0 or species richness with RNA-Seq). For D_max, on average, the differences among tissue types are around 5%, and around 25% (16%-34%) among cancer types. Additionally, the differences in D_max among cancer types also vary in terms of sequencing protocol (RNA-Seq < WGS) and tissue types (PT > SN).

This level of diversity analysis performed for tumor microbiomes summarized here mirrors the population-level diversity-disease relationship (DDR) previously reported by Li and Ma (2021), in which they found, based on the analysis of 23 microbiome-associated diseases, that the population-level DDR or p-DDR was only significant in approximately 5% of cases for the parameter D_max, and was insignificant for other DAR parameters. A significant relationship means that the parameters are different between the disease and health states, or between PT and SN or between PT and B. On this point, both the diversity-cancer relationship (DCR) and general DDR at the host population level show the same pattern—virtually all major DAR scaling parameters are invariant (especially z) except for D_max, which is variable in approximately 5% of cases. Note that the p-DDR or p-DCR are different from DDR/DCR at the individual host level, or microbial alpha-diversity level, where Ma et al. (2019) found that the DDR relationships were significant in approximately 1/3 of the cases they studied. Both studies analyzed the same datasets of 23 microbiome-associated diseases but generated different results. Li and Ma (2021) postulated that it should be the mutual cancellations of individual-level DDR differences (ups and downs) that generated the ‘flat’ or insignificant scaling parameters at the host population level. Here, we believe that the findings on p-DCR in this study simply cast supportive evidence on that previous hypothesis—that mutual cancellations of the ups and downs at the individual host level differences are responsible for the general lack of differences at the host population level.

In this study, our p-DDR level analysis goes beyond single disease or single cancer type; instead, our diversity-scaling analyses were performed at the individual cancer type level (for comparing different tissue types such as PT vs. SN or B) and at the pan-tumor level (across cancer types for comparing different cancer types such as lung vs. breast cancer), respectively. This pan-disease analysis was not performed in previous studies by Li & Ma (2021) and Ma et al. (2019) because those microbiome datasets were not from single tissue types—instead they included rather heterogeneous samples from gut, oral, skin, and vaginal microbiomes. In this study, all microbiomes were from human tissues or blood. Although we consider the tissue microbiome samples to be more homogeneous than those analyzed in the previous studies (Ma et al. 2019, Li & Ma 2021), there should still be a certain level of heterogeneity among different tissue types. We postulate that the relatively larger D_max differences (around 25%) at the pan-tumor level among cancer types, compared to approximately 5% difference at the tumor level, should be attributed to the heterogeneity among different cancer types including different tissues such as lung versus breast tissues.

Given the general lack of differences in major scaling parameters, especially z, at the pan-tumor level analysis, we further built pan-tumor DAR models (Table S7) by combining all disease types. Further tests (Table S8) of the pan-tumor population DAR models (pp-DAR models) revealed no significant differences between tissue types (PT, SN, or B), which prompted us to further build pp-DAR models by combining tissue types (Tables S7 & S8). In Table S7 and Table 2 (summary version), we show that the total species richness, typical species equivalents, and dominant species equivalents of human tissues including blood are approximately 1948, 36, and 22, respectively.

Among the total numbers, the proportions of archaea, bacteria, and viruses are about 5%, 78%, and 17%, respectively, in terms of RNA-Seq. The proportions are slightly higher (about 1%) for bacteria and viruses in terms of the WGS protocol, and that for archaea would be approximately 3% only.

To the best of our knowledge, there is no existing study that has estimated the total (gamma) microbial diversity at the host population level for tumor tissue microbiomes, or actually for any human tissue microbiomes. Since no significant differences were detected between tumor tissues and solid normal tissues (or even blood-derived normal), our results with combined cancer types and tissue types can be considered an estimate of gamma diversity of the human tissue microbiome. Although no similar estimates exist for human tissues, there were similar estimates for non-tissue human microbiomes. In fact, the development of the DAR modeling method extending the classic SAR (species-area relationship) was demonstrated with Human Microbiome Project (HMP) datasets (HMP Consortium 2012; Ma 2018a, 2018b, 2019). We choose two human microbiome sites, skin and gut (stool), to compare with the tissue microbiome. For species richness (q=0), the total species in the tissue microbiome is approximately 1⁄4 (25%) of the gut microbiome (7706) and 1/8 (12.5%) of the skin microbiome (16206). For the number of typical species (q=1), the number of typical species equivalents in tissue microbiomes is 3.5% of gut microbiomes (1020), and 6% of skin microbiomes (605). For the number of dominant species (q=2), the number of dominant species equivalents is 8% of the gut microbiome (272), and 40% of the skin microbiome (55). The higher similarity between tissue and skin microbiomes than between tissue and gut microbiome in terms of dominant species equivalents is somewhat puzzling and deserves further investigation, given that one is ‘internal’ and another is ‘external’ and therefore should be dissimilar. In contrast, for the total species number, the tissue microbiome is indeed more similar to the gut microbiome than to the skin microbiome. Note that our comparisons here are purely based on numbers of species, rather than on compositional comparisons. It is likely that the comparisons are not very informative and instead simply offer relative numbers of species equivalents.