Machine learning models reveal microbial signatures in healthy human tissues, challenging the sterility of human organs

Machine learning models reveal microbial signatures in healthy human tissues, challenging the sterility of human organs | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Machine learning models reveal microbial signatures in healthy human tissues, challenging the sterility of human organs ANARGYROS SKOULAKIS, GIORGOS SKOUFOS, ARMEN OVSEPIAN, ARTEMIS G. HATZIGEORGIOU This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4610181/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 25 Jan, 2025 Read the published version in Frontiers in Microbiology → Version 1 posted You are reading this latest preprint version Abstract Background The presence of microbiomes within healthy human internal organs still remains under question. Our study endeavors to discern distinctive microbial signatures within normal human internal tissues from the Genotype-Tissue Expression (GTEx) consortium through the application of machine learning (ML) models. Machine learning models were developed to classify each tissue type against all others using solely the microbial profile of the tissues. The identification of tissue-specific microbial signatures that enable machine learning models to accurately predict different tissue types would suggest the presence of distinct microbial communities within tissues. To identify these microbial signatures, we analyzed 13,871 normal RNA-seq samples from 28 tissues obtained from the GTEx database. Results Tissue-specific microbial signatures were identified in 11 out of the 28 tissues while the signatures for 8 tissues (Muscle, Heart, Stomach, Colon tissue, Testis, Blood, Liver, and Bladder tissue) demonstrated resilience to in silico contamination. As the GTEx samples were derived from post-mortem biopsies, healthy living tissue samples were also analyzed to further validate the results. The models for Heart, Colon tissue, and Liver displayed high discriminatory performance in the living dataset, suggesting the presence of a tissue-specific microbiome for these tissues even in a living state. Notably, the most crucial features were the fungus Sporisorium graminicola for the heart, the gram-positive bacterium Flavonifractor plautii for the colon tissue, and the gram-negative bacterium Bartonella machadoae for the liver. The presence of contamination-resilient tissue-specific microbial signatures implies that microbial colonization is not a random occurrence. Conclusions The presence of tissue-specific microbial signatures in certain tissues suggests that these organs are not devoid of microorganisms even in healthy conditions and probably they harbor low-biomass microbial communities unique to each tissue. The discoveries presented here confront the enduring dogma positing the sterility of internal tissues, yet further validation through controlled laboratory experiments is imperative to substantiate this hypothesis. Exploring the microbiome of internal tissues holds promise for elucidating the pathophysiology underlying both health and a spectrum of diseases, including sepsis, inflammation, and cancer. human tissues human organs tissues microbiome microbial communities microorganisms microbial signatures machine learning Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Background Microbes exhibit a remarkable spectrum of functions and capabilities that enable them to colonize diverse and extreme habitats. The human body serves as a host for a vast array of microbes with varied functionalities 1,2 . Our understanding regarding interactions between humans and microbes has evolved through the progression of sequencing technologies. Pioneering this transformative trajectory was the NIH Human Microbiome Project (HMP), which revealed the complex interplay between human and microbial cells 3,4,5 . Recent advancements have prompted a reevaluation of long-standing beliefs, showcased by the dismissal of the traditional notion that healthy human lung and vagina, are sterile environments. 6,7 The predominant focus of microbiome studies has been on characterizing microbial communities in easily accessible sampling sites, including the human intestinal tract, the oral cavity, and the skin. Internal tissues, less accessible by conventional sampling methods, have historically been presumed to be devoid of microorganisms due to protective layers of epithelial and endothelial tissues. However, this assumption has been challenged, particularly in settings such as cancer, where various studies have revealed the presence of tumor-specific microbial sequences within human tumors, offering a potential new avenue for cancer diagnosis 8–14 . As the exploration of microbes inside internal tissues has primarily been conducted in the context of pathological conditions like cancer, fundamental questions regarding healthy status persist. Are there microbial residents within human tissues engaged in a commensal relationship with human cells? What is the role of their presence? Can these microbial signatures trigger an immunological response or serve as diagnostic markers for organ health? Gita Mahmoudabadi et al. observed a notable presence of bacterial genera in tumors that is also detectable and similar in adjacent tumor-free tissues 15 . This suggests that tumor microbiomes may partially originate from neighboring normal tissues, indicating that, under healthy conditions, internal tissues host a microbiome. Hieken et al also detected microbial communities in aseptically collected human breast tissues in benign and malignant conditions 16 . Moreover, investigations in germ-free mice have unveiled microbiomes within various organs, including the brain, muscle, adipose tissue, liver, and heart. This challenges the conventional notion that the internal organs of mammals are devoid of microbial presence. 17 To elucidate the microbiome inside healthy human tissues, we conducted a re-analysis of RNA-seq sequencing data derived from Genotype-Tissue Expression (GTEx) consortium 18 . Originally designed to explore variations in gene expression within healthy human tissues, GTEx's extensive repository contains samples exclusively derived from healthy specimens. Acknowledging certain limitations within the GTEx analysis pipeline, such as the employment of a poly-A selection protocol and the reliance on post-mortem biopsy samples, it remains noteworthy that the GTEx consortium provides the most extensive and analytically robust dataset about RNA expression within healthy human tissues. To date, GTEx sequencing data have not been explored from a microbiome perspective. Here, we present what, to our knowledge, stands as the inaugural and comprehensive healthy human tissue microbiome dataset. Leveraging ML models, we identified microbial signatures capable of discriminating among various tissue types and investigated the potential influence of phenotypic traits (e.g., age, sex, BMI) on these signatures. The presence of tissue-specific microbial signatures in certain tissues would suggest that these organs are not devoid of microorganisms even in healthy conditions, thus challenging the long-standing dogma of the sterility of internal tissues. Results Low-biomass microbial RNA is detected in GTEx samples A total of 13,871 RNA-seq sequencing samples originating from 28 distinct tissues (e.g., Lung, Liver, Pituitary, Blood Vessel, Thyroid, Skin, Salivary Gland, Esophagus, Heart, Muscle, Pancreas, Adipose Tissue, Vagina, Blood, Ovary, Spleen, Prostate, Adrenal Gland, Nerve, Stomach tissue, Colon tissue, Testis, Brain, Breast, Uterus, Small Intestine tissue, Kidney, Bladder tissue) within the GTEx consortium have been utilized. All samples were designated as non-diseased and free of pathology, as the specimens were reviewed by a panel of 2–3 pathologists and any specimen found with an incidental finding had been systematically excluded. The workflow of the study is shown in Fig. 1 . Following the removal of reads of human origin, the remaining reads were aligned against an extensive database comprising 4,034 bacterial, 489 archaeal, 11,259 viral, and 81 fungal entities. The microbial composition of each sample was determined utilizing AGAMEMNON 19 . AGAMEMNON provides precise genus, species, and strain abundances through an efficient indexing scheme for rapid pattern matching, facilitating analysis of extensive datasets using common computational resources. The abundance estimation employs expectation maximization algorithm and targets maximizing the likelihood of the observed reads by gradually altering the abundance value associated to different taxa. On average, out of ~ 93×10 6 raw reads per sample, ~ 2×10 6 reads (2.35% of the raw reads) did not align to the human genome. After quality control (QC) and filtering of the non-human reads, ~ 3.5×10 5 reads (0.38% of the raw reads) were preserved. Employing Pufferfish 21 , ~ 4.3×10 4 (0.046% of the raw reads) were attributed to bacteria, archaea, fungi, or viruses. To eliminate any residual human-associated sequences, a realignment of microbial reads to the most recent human genome assembly (T2T-CHM13v2.0) was executed, resulting in ~ 4.2×10 4 reads (0.045% of raw reads) remaining. Finally, using AGAMEMNON, all the remaining reads (4.2×10 4 reads, 0.045% of raw reads) were successfully classified into microbial species (Fig. 2 a). The reported read counts are the mean values across all samples at each step of the analysis. The bladder exhibited the most pronounced abundance of microbial reads, whereas the adrenal gland, brain, heart, muscle, salivary gland, and testis demonstrated the lowest microbial read counts, as illustrated in Supplementary Fig. 1 and in Supplementary File 1 . In each of the 28 tissues examined, a discernible presence of low-biomass microbial RNA was identified. The microbial communities detected within these tissues may either genuinely inhabit the respective tissues or be a result of potential contamination. As in the case of microbial reads, the bladder exhibited the highest microbial species richness in its microbial community, whereas the brain, heart, muscle, salivary gland, and testis displayed the lowest richness, as depicted in Fig. 2 b. Notably, the Shannon index of the bladder was markedly lower compared to the rest of the samples (Fig. 2 c ) , which arises from the fact that bladder samples demonstrated elevated counts for a few selected species, and comparatively lower counts for the remaining species.. To identify species that are consistently present in each tissue and are not sporadic opportunistic pathogens, we retained species present in at least 10% of samples from each tissue. The retained species are hereafter referred to as the core microbiome of the tissue. Species meeting this criterion for at least one tissue were included in subsequent analyses, resulting in a total of 1,708 different species across all four kingdoms. As expected, the tissue with the smallest core microbiome was the brain, consisting of 720 species, while the bladder exhibited the most diverse core microbiome, comprising 1,344 species (Fig. 2 d and Supplementary File 2 for details). Additionally we noticed that across all tissues, bacteria constituted the predominant component, making up over 95% of the core microbiome for each tissue. The core microbiome of each tissue for all four kingdoms, namely bacteria, fungi, viruses and archaea, is documented in Supplementary Files 3 to 6 respectively. Notably, a total of 529 species, accounting for ~ 31% of the overall core microbiome, were shared across all tissues, showing that tissues share in some extent a common microbiome (Fig. 2 e). Among fungi, the majority (24 out of 37 species) were present in all tissues, as indicated in Supplementary Fig. 1 and Supplementary File 7 . However, for bacteria, only 504 out of 1,297 species (38.86%) and for viruses, only Geobacillus virus E2 out of 13 viruses were present in all tissues, as illustrated in Supplementary Fig. 2, 3 and Supplementary Files 8, 9 , respectively. Archaea were exclusively identified in specific tissues, including the bladder and kidney (species Methanocaldococcus jannaschii ), the colon and small intestine tissues (species Methanosarcina sp. WH1 ), and the spleen (species Methanocaldococcus jannaschii ) ( Supplementary File 6) . 8 tissues harbor specific microbial signatures at species level While solely the identification of microbiomes within tissues doesn't warrant the existence of microbial communities within these tissues, the presence of discernible biological footprints within these communities would imply a non-coincidental occurrence of microbes inside these tissues. To explore this, Machine Learning models were employed to detect potential microbial footprints for each tissue. By using Cum-Sum Scaling method 21 to normalize microbial profiles across samples, stochastic Gradient Boosting Machine Learning models (GBM models) 22 were trained to distinguish individual tissues from all other tissues. Notably, 11 out of the 28 models (models for Brain, Small Intestine tissue, Liver, Bladder, Muscle, Heart, Salivary Gland, Stomach tissue, Colon tissue, Testis, and Blood) exhibited robust performance (mean AUROC > = 0.70 and relative mean AUPR > = 1.4) discriminating each specific tissue type from the collective representation of all other tissues, as illustrated in Fig. 3 a and in Supplementary File 10 ( Methods section ML models 1vs27AllOtherTissues). There was no significant correlation between the sample size and mean AUROC performance (rho = 0.207, p-value = 0.28) and the normalized mean AUPR performance (rho=-0.152, p-value = 0.43). The identified microbial signatures based on species abundances within these 11 tissues were further analyzed to investigate if the microbial presence could be attributed to sample contamination. Given that the primary objective of the GTEx initiative was to investigate gene expression, the samples were not processed under sterile conditions. Due to the nature of GTEx experiment protocol, it was inefficient to employ an in-silico decontamination approach and to overcome this inefficiency, an alternative strategy was implemented. A randomized in silico contamination approach was used to evaluate the potential impact of systematic contamination on the performance of the models. To scrutinize whether possible contamination could drive the performance of ML models exhibiting high performance, 12 pseudo-contaminants were strategically introduced to the data. These pseudo-contaminants aimed to simulate diverse scenarios of contamination throughout the analytical process, encompassing potential contamination scenarios from the initial sampling to the sequencing procedure ( Methods section “In silico Contamination”). We retrained the ML models of the 11 tissues presenting high discriminatory power (mean AUROC > = 0.70 and relative mean AUPR > = 1.4) and compared their performances with and without in silico contamination ( Supplementary file 11 ). Only the models of Salivary Gland, Brain, and Small Intestine had significantly higher performance (both mean AUROC values and mean AUPR values) in the contaminated iterations, with in silico contaminants playing an important role in the models’ performance. On the other hand, for the rest of the models in the contaminated iterations, models did not rely on the inserted pseudo-contaminants (Fig. 3 b), Therefore, the models for the Salivary Gland, Brain, and Small Intestine tissues were deemed unreliable as their high performance could be the result of contamination. The contamination type, that could potentially influence the performance of these 3 models, is a low volume contamination (contamination type 8) which imitates a contaminant that is systematically present in the samples of a collection site center ( Methods section “In silico contamination”). Nevertheless, in the rest 8 tissues (Blood, Testis, Colon tissue, Stomach tissue, Muscle, Bladder, Liver, and Heart), the tissue-specific microbial signatures exhibited resilience against potential contamination. The ability of GBM models to discriminate these tissues solely based on their microbial profiles independently of any potential contaminations, underscores the presence of distinctive microbial signatures specific to each of these tissues. To investigate microbial differences among these 8 tissues, GBM models were recreated using data exclusively from these contamination-resilient tissues, distinguishing each tissue from the combined representation of all the other 7 tissues ( Methods section “ML models 1vs7Tissues”). Remarkably, all 8 tissues demonstrated robust performance in terms of AUROC and relative AUPR (mean AUROC > = 0.70 and relative mean AUPR > = 1.4), as illustrated in Fig. 3 c and in Supplementary File 12 . Collectively, these results suggest that each of these 8 tissues harbors a distinct and unique microbiome. To assess the biological relevance of these microbial signatures, a detailed examination of the most important features from the 8 models was conducted. Boxplots of the abundance of the most important features across the different tissues are depicted in Fig. 3 D- 3 K. Supplementary File 13 contains the full catalog of features importance scores for each tissue. In the colon tissue model, the 2 most predominant features were the bacterium Enterocloster bolteae , formerly known as Clostridium bolteae , a recognized constituent of human feces and Flavonifractor plautii (formerly Eubacterium plautii ), which has been isolated from human feces, blood, intra-abdominal pus, and infected soft tissues in humans 24 (Figs. 3 D, 3 E). In the case of Heart, the most crucial feature was the fungus Sporisorium graminicola and for the Blood model, the 2 most important features were Schizosaccharomyces pombe and Sporisorium graminicola (Fig. 3 F, 3 G respectively). These two fungal species, discussed in the study by Hidayat et al. , are capable of synthesizing mannosylerythritol lipids (MELs) 25 . MELs belong to the glycolipid class of biosurfactants and are known for their outstanding interfacial and biochemical characteristics, as highlighted by Morita et al. 26 . For the Muscle model, the overall importance of features was relatively subdued, and many features contributed equally lightly in the model’s discriminatory ability. In the Stomach model, the pivotal feature was the virus Human betaherpesvirus 7 , previously identified as an inhabitant of gastric mucosa 27 (Fig. 3 H). For the Testis model, the critical feature was the Human endogenous retrovirus K (HERV-K) (Fig. 3 I). The human genome harbors numerous copies of HERV-K, many of which retain intact open reading frames (ORFs). These ORFs are capable of being transcribed and translated, particularly during early embryonic development and in cancerous conditions. 28 For the Bladder model, the foremost feature was the bacterium Caldimonas thermodepolymerans , an underexplored microorganism phylogenetically proximate to the Comamonadaceae group that was systematically less abundant in the samples of bladder (Fig. 3 J). Lastly, in the Liver model, the primary feature was Bartonella machadoae (Fig. 3 K). Bartonella machadoea is a bacteria belonging to Proteobacteria genus , and it was recently reported that the liver is inhabited mainly by proteobacteria by a gut-liver-specific axis 29 . For additional validation, the samples from the aforementioned eight tissues underwent analysis using Kaiju, a different computational method for microbial taxonomic profiling. Kaiju assigns each sequencing read to a taxon in the NCBI taxonomy by comparing it to a reference protein database. The used reference database was the microbial subset of the NCBI BLAST non-redundant protein database, encompassing bacteria, archaea, virus, fungi and microbial eukaryotes. Subsequently, GBM models for these 8 tissues were reconstructed using the taxonomic profiles generated by Kaiju. Remarkably, the performance of these models exhibited a high performance (mean AUROC > = 0.70 and relative mean AUPR > = 1.4) as depicted in Fig. 3 L ( Supplementary file 14 ), indicating that regardless of the bioinformatics method analysis, these tissues present a distinctive microbial signature. The most important features of the microbial signatures found by Kaiju ( Supplementary file 15 ) were different compared to the microbial profiles derived by AGAMEMNON. Employing different approaches for read classification, the taxonomic profiles and respectively the most important features of the models are differing, mainly due to the different reference databases that each tool uses for the classification of the reads and to the algorithm used in assigning the reads to each organism. Heart, Colon tissue, and Liver preserve the tissue-specific microbial signatures in the living state As the specimens sourced from the GTEx consortium originate from post-mortem biopsies, an investigation was conducted to examine whether the identified signatures for the eight tissues resulted from post-mortem microbial colonization or were present during the subjects' lifetime. To address this, the eight models, comparing one tissue against the seven others, were subjected to testing using data derived from an entirely distinct project (NCBI ID PRJEB4337). This dataset sourced from living tissues of subjects within the Swedish Biobank 30 . This project did not contain samples from muscle tissue and blood, so only 6 out of the 8 models (heart, colon, stomach, liver, bladder, testis models) were tested with data from living tissues. The data were processed utilizing the same pipeline as applied to GTEx data, and normalization was separately carried out on this specific distinct dataset using the CSS normalization method. For each tissue, all 100 iterations of the 1vs7Tissues models (that were created before) were tested again using the living dataset. The obtained results were then compared with the performance of models, generated with randomly assigned tissue labels to ascertain the significant contribution of tissue labeling to model performance ( Supplementary file 16 ). Remarkably, 5 out of 6 tissues (only the model for bladder tissue did not) exhibit a statistically significant difference in AUROC and AUPR in the living dataset (p-value < 0.05) between the true tissue labeled model and the random tissue labeled model. Among the six models, the models for Heart, Colon tissue, and Liver demonstrated robust AUROC and relative AUPR in the living dataset (AUROC > = 0.7 and relative AUPR > = 1.4) (Fig. 4 a, Supplementary file 17 ). This suggests that the microbial signatures identified in these tissues (Heart, Colon tissue, and Liver) also exist in the healthy living state. The 5 most abundant species in the living dataset for the 3 tissues are depicted in Fig. 4 B-D (refer for frequencies of all microbiomes inside living tissues of the project PRJEB4337 in Supplemetary file 18 ). The lack of high performance for the remaining three tissues may be attributed to the likelihood that the microbial signature of these tissues undergoes significant alterations during the event of death. The high performance of the Heart, Colon tissue, and Liver models was not observed in the models created with the taxonomic profiles derived from Kaiju ( Supplementary File 19 ). Lifestyle factors and medical history do not influence tissues’ microbiome Additionally, The putative impact of lifestyle factors on the composition of tissue-specific microbiomes across the eight human tissues was examined. The following variables were studied : Age, Body Mass Index (BMI), Sex, Alcohol consumption, Smoking status, Ancestry). To elucidate potential associations, GBM models were constructed for each tissue to predict the specific trait. However, despite rigorous computational analyses, the findings reveal a lack of robust predictive capability of the tissue microbiomes for these traits across all examined tissues. Detailed performance metrics for each trait within each tissue are provided in Supplementary File 20 . In addition, each tissue underwent a thorough examination to discern whether its microbial profile harbored predictive potential for various medical conditions, encompassing Alzheimer's or Dementia, Asthma, Cerebrovascular Disease, COPD-CLRD, Depression, Diabetes Type I and II, Dialysis Treatment (Renal Failure), Rheumatoid Arthritis, Hypertension, Ischemic Heart Disease, Liver disease, Schizophrenia, Seizures, and Arthritis (as a comprehensive category inclusive of various types of arthritis). Only in the case of dialysis treatment (renal failure) could discernment be achieved between heart tissues of subjects with renal failure and heart tissues of subject without renal failure. However, for the rest of the diseases, discrimination of subjects’ diagnosis across the tissues was not feasible. Detailed performance metrics for each disease in each of the eight tissues are presented in the Supplementary file 21 . Overall, the analysis of tissue microbiomes suggests a lack of inherent microbial signatures indicative of lifestyle factors or systemic disease status. The utmost significance in unraveling tissue-specific microbial signatures lies within the species-level analysis To elucidate the optimal taxonomic level for investigating microbial signatures across the eight vital tissues, we constructed Gradient Boosting Machine (GBM) models utilizing microbial profiles at the genus level. AGAMEMNON supports the aggregation of microbial counts at the taxonomic level of choice. Employing microbial profiles at the genus level, we developed GBM models for the aforementioned eight tissues. GBM models based on genera exhibited slightly inferior performance compared to those derived from species-level data (refer to Supplementary File 22 ). Although the distinctions are subtle, the majority of the eight models demonstrated enhanced performance when analyzed at the species level; only, the bladder model exhibited superior mean AUROC and mean AUPR at the genera level. To investigate if any other of the twenty tissues, that did not present microbial footprints at species level analysis, manifested distinct microbial signatures at the genera level, we generated GBM models utilizing taxonomic profiles at the genera level for each tissue. Only the models from the same tissues as species level models presented high performance (AUROC > = 0.7 and relative AUPR > = 1.4), most of the times slightly diminished when compared to their species-level counterparts. Only the brain tissue models did not present high performance at the genera level, but as previously shown, the species-level brain model is possibly a result of contamination (refer to Supplementary file 23 ). In light of the known redundancy in the functional repertoire (i.e. genes) across diverse microorganisms, we endeavored to investigate whether tissue-specific microbial signatures predominantly arise from a tissue-specific microbial gene repertoire rather than a specific microbial taxonomy. To delineate the functional profile of samples, microbial reads identified by AGAMEMNON were subjected to analysis using Humann3 31 , specifically employing translational searches against the Uniref90 database. Subsequently, GBM models for the aforementioned eight tissues were reconstructed utilizing the CSS-normalized functional profiles of the samples (refer to Supplementary File 24 ). Remarkably, only the Blood, Colon, and Bladder tissues exhibited high performance, with a mean AUROC > = 0.70 and a relative mean AUPR > = 1.4. In an attempt to integrate both functional and species-level information, we reconstructed GBM models using both the functional and the taxonomic profiles of the samples. However, even with this comprehensive approach, a very slight enhanced performance was attained and not in all models (see Supplementary File 25 ), indicating that in our scientific setup the information of functional repertoire contributes very slightly in discriminating tissues. It is crucial to note that the low performance of the functional models, could most probably be attributed to the technical characteristics of the GTEx RNA-seq data. The data, generated using a poly-A selection sequencing protocol, inherently filtered out a significant proportion of microbial RNAs. Discussion In each of the 28 examined tissues, a modest yet substantiated presence of transcriptionally active microbial communities was observed. Contrary to traditional notions of sterility of internal tissues, this analysis, alongside analogous investigations in other mammalian species, suggests that tissues potentially harbor a low-biomass microbiome not only in the context of disease but also in normal state. Established knowledge in the field of cancer underscores the presence of a tumor microbiome, with mounting evidence indicating its pivotal role in tumor development and influencing treatment outcomes 32 . Extending beyond pathological conditions, our analysis suggests that certain human internal tissues consistently maintain a low-biomass microbiome. This phenomenon prompts speculation that the microbiome may serve as a functional reservoir contributing to tissue well-being or act as a regulator of the immune system. The robust discriminatory performance demonstrated by the models for eight tissues, relying solely on microbial profiles of samples, strongly suggest that the observed microbiome is not a result of random chance. Despite the impracticality of conducting in silico decontamination due to the uniform utilization of the same RNA quantity for each sample in the processing of GTEx samples and the absence of negative controls, a reverse analysis —implementing in silico contamination— illustrated the robustness of microbial signatures in eight tissues against various potential contaminations. It is imperative to clarify that the in silico contamination approach does not seek to substitute a laboratory analysis conducted under the most stringent sterile conditions, nor is it exhaustive in considering all conceivable contaminations. Rather, it serves as a validation step to evaluate whether contamination occurred at different stages of the analysis and to assess its potential impact on models performance. In total, eight tissues demonstrated a tissue-specific microbiome that exhibited resilience against various potential contaminations, thereby indicating the presence of a distinct microbiome in each tissue. Interestingly, this tissue-specific microbiome primarily diverges not in terms of species but predominantly in the composition of these species. The microbial signatures for heart, colon and liver tissues seem to be distinctive also in an independent dataset from healthy living tissues, supporting the presence of microbiomes inside these tissues also in a living healthy state. In colon tissue, the most pivotal component of its microbial signature is Flavonifractor plautii , a prevalent bacterium in the human gastrointestinal tract recognized for its notable butyrate production 33 . In liver tissue, the primary microbial entity of significance is Bartonella machadoae . Bartonella species have been associated with compromised liver function 34 . Lastly, the presence of the phytopathogenic fungus Sporisorium graminicola in cardiac tissue presents an intriguing finding. S. graminicola has also been isolated from human fecal samples 35 . Despite the apparent peculiarity of a phytopathogenic fungus inhabiting cardiac tissue, its consistent identification as a significant feature in both the GTEx dataset containing samples from USA and the independent test dataset from Swedish Biobank renders the possibility of mere contamination or chance occurrence highly unlikely. The ability of machine learning models to discriminate distinct microbial signatures across these three tissues in separate datasets suggests that these microorganisms are likely indigenous inhabitants of their respective tissues. In the majority of tissues (20 out of 28 tissues), a distinct tissue-specific microbial signature was not detected. However, the absence of such a signature does not necessarily imply tissue sterility. These tissues are more likely to harbor distinct microbial signatures in species that may not have been successfully detected, possibly due to limitations imposed by the poly-A protocol or the lack of reference genomes. Indeed, the methodologies employed by GTEx protocols unequivocally revealed only a fraction of the microbial load present inside these tissues. Another plausible explanation for the absence of tissues' microbial signatures is that the models for these tissues may not have performed optimally, as these tissues may exhibit a more versatile microbiome without a clearly defined tissue-specific signature. Lastly, for the three tissues that did not perform well in living tissues, it is conceivable that these tissues undergo a significant shift in their microbiome after death, or there may be notable differences attributable to geographical variations, as the living subjects were from a different continent than GTEx subjects. Significantly, beyond the considerations of the poly-A sequencing protocol, metagenomics and metatranscriptomics necessitates tailored laboratory procedures for efficient DNA or RNA isolation from both gram-positive and gram-negative microorganisms. As these considerations were not initially integrated into the GTEx pipeline, it is evident that the analytical procedure failed to unveil the full spectrum of microorganisms present in the tissues. Moreover, aside from the critical role of laboratory protocols, the bioinformatics analysis protocols also play a pivotal role in discriminating and accurately describing tissue microbiomes. Two metatranscriptomics tools, Kaiju and AGAMEMNON, each employing distinct approaches for read classification, yielded different taxonomic profiles and exhibited variations in performance on the models. To advance microbiome research, it is paramount to develop specific guidelines for both laboratory and in silico analysis of microbiome data. This step is crucial for ensuring comprehensive and accurate insights into the diversity and composition of microbial communities within tissues. It is crucial to consider that, in the majority of the identified microbiomes, the reads classified as belonging to these microbiomes were relatively scarce. This implies that RNA transcripts of these microbiomes were present within the tissues; however, these transcripts did not adequately cover a substantial portion of the genetic material of these microbiomes. This limitation is likely attributable to constraints imposed by the poly-A protocol. The limitations of the poly-A protocol highlight the necessity for additional investigations to validate the existence of microbiomes within tissues, rather than merely isolated microbial RNA transcripts. It is pertinent to acknowledge that bioinformatic analyses in scenarios such as the present study are subject to probabilistic constraints and, on their own, are insufficient to definitively prove the existence of microbiomes. However, they can serve as valuable tools to guide subsequent investigations, indicating tissues and organisms that warrant further scrutiny. Conclusion Collectively, a growing body of evidence substantiates the existence of microbial interactions with human cells in anatomical sites traditionally considered sterile. The identification of encoded signatures specific to certain tissues within normal human tissues suggests for the first time that these tissues harbor a low-biomass microbiome. This exploration signifies a paradigmatic shift, heralding an era in which we acknowledge that our bodies are not solitary entities but rather collaborative ecosystems housing diverse microbial species. We anticipate that our study will serve as a foundational resource, providing crucial guidance for future investigations and facilitating targeted laboratory validations aimed at confirming the presence of microbial communities within internal tissues. Methods GTEx data accession The data used for the analyses described in this manuscript were obtained from dbGaP accession number phs000424.v9.p2 on 04/01/2023. The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. All GTEx RNA-seq data and matched samples protected metadata were accessed via NHGRI Analysis Visualizaion and Informatics Lab-space (AnVIL )( https://gen3.theanvil.io/login ). Details of how these data were downloaded are comprehensively described in the AnVIL site ( https://anvilproject.org/learn/reference/gtex-v8-free-egress-instructions#step-4---download-object-files ). For bulk data acquisition, custom Snakemake files were employed. Due to constraints in storage capacity, an approach was implemented wherein individual sample BAM files were downloaded, and only the unmapped reads (as was originally mapped by GTEx pipeline using STAR (version 2.5.3a) 36 https://gtexportal.org/home/methods ) were retained using samtools (version 1.10) 37 . GTEx sample and subject QC The original GTEx dataset initially comprised 17,350 RNA-seq samples derived from 948 subjects, all of which had successfully met the GTEx original quality control (QC) criteria. Subsequently, 52 samples were excluded from the analysis due to missing data in critical variables, namely, Genotype or Expression Batch ID, Date of genotype or expression batch, and Total Ischemic time. The dataset was further refined to encompass only 14,478 samples originating from tissues preserved in the PAXgene tissue fixative solution, specifically categorized under "PAXgene" and "Whole Blood:PAXgene" in the "current_material_type" variable. Furthermore, an additional refinement process involved excluding 32 samples originating from tissue sites with fewer than 20 samples, specifically samples from the following tissues : Kidney Medulla, Fallopian Tube, Cervix Extocervix, and Cervix Endocervix were removed. Additionally, 3 samples possessing an RNA Integrity Number (RIN) less than 5 were excluded, along with 139 samples with Total Ischemic time exceeding 1440 minutes (24 hours). Following these rigorous sample quality control measures, a total of 14,304 samples from 942 subjects and spanning 28 different tissues remained in the dataset. Following, quality control of the subjects was performed. 9 subjects were initially excluded due to their ineligibility based on GTEx original criteria. 7 subjects were then removed from the dataset as they had a current cancer diagnosis, and 16 subjects were excluded due to a history of cancer diagnosis within the past 5 years. All the samples from the excluded subjects were removed. Following this comprehensive sample and subject quality control process, the resulting dataset comprised 13,871 RNA-seq samples originating from 28 distinct tissues (refer to S upplementary file 25 ) and derived from a cohort of 910 subjects. The quality control procedures were implemented utilizing custom Python scripts, available on the corresponding GitHub repository. Isolation of unmapped reads and quality control For each sample, sequencing reads that failed to align with the human reference genome, as indicated by the mapping information in the raw BAM files obtained from GTEx, were selectively retained. To isolate the unmapped reads where both paired reads were unaligned, and to eliminate reads classified as non-primary alignments, bioinformatic tool Samtools was employed using the arguments “-f 12 -F 256 “. The unmapped reads underwent a comprehensive quality control process, involving the exclusion of reads with a length shorter than 35 nucleotides (nt), given that the GTEx normal sequencing length was 76 nt. Additionally, steps were taken to remove adapters and perform quality trimming with a threshold at Phred quality score of 15. These quality control procedures were executed using Atropos (version 1.1.31) 38 . To streamline and automate this sequence of steps, custom Snakemake scripts were implemented, available on the corresponding GitHub repository. Taxonomic Assignment using AGAMEMNON To construct the reference database utilized for the taxonomic algorithm AGAMEMNON (version 0.1.0), a custom Bash script was employed on 02/14/2023 to download microbial genomes from RefSeq. The dataset comprised all bacterial representative or reference genomes with complete genome assembly level (4,034 bacterial genomes), all archaeal genomes with complete genome assembly level (489 archaea genomes), and all viral genomes with complete genome assembly level (11,259 viral genomes). Furthermore, fungal transcriptomes from representative or reference genomes with complete genome or chromosome assembly levels (81 fungal transcriptomes) were also included in the reference database. In the case of fungi, the transcriptome data was specifically employed to account for the intricacies of the splicing process. This approach ensures a more comprehensive representation of fungal genomic information, taking into consideration the variations introduced during the splicing of transcripts. The quality-controlled non-human reads of each sample underwent mapping against the custom reference database described above, using the Puffaligner algorithm 39 . For Puffaligner, the flag “--noOrphans” was used in order to discard the orphans reads. The term “orphan” refers to one end of paired-end read that is confidently aligned to some genomic position, but for which the other read end is not jointly aligned nearby (and paired). Puffaligner aligns the reads to the compiled microbial genomes, enabling the identification of reads of microbial origin within the samples. Puffaligner is a fast, sensitive and accurate aligner based on a compacted sequence graph and is meticulously crafted to embody a dual emphasis on high sensitivity in alignment tasks and efficient computational performance. Its design capitalizes on the utilization of a colored compacted de Bruijn graph to efficiently identify and factor out recurring sub-sequences within the reference. As highlighted by Gitawi et al. 40 , in order to be sure that no human read had succeeded to infiltrate in our analysis, we remapped the reads that were classified by Puffaligner as reads of microbial origin to the most recent human genome assembly (T2T-CHM13v2.0 from T2T Consortium) using bowtie2 (version 2.2.3) 41, 42 with the "preset" parameters of “--very-sensitive” and kept all the reads that were not mapped concordantly. Custom snakemake scrips were used to automate the process. Then by using only the remaining reads, we analysed them with AGAMEMNON (a changed version, deposited in the github repo of the present study) to find the microbial profile of each sample. AGAMEMNON represents a metagenomics and metatranscriptomics algorithm, notable for its integration of a time and space-efficient indexing scheme. This feature facilitates rapid pattern matching, allowing for the efficient indexing and analysis of extensive datasets using commonly available computational resources. In the abundance estimation step, the primary approach relies on the expectation maximization (EM) algorithm. The goal is to maximize the likelihood of observed reads by iteratively adjusting the abundance values linked to various taxa. At last, after running AGAMEMNON, the taxonomic profile of each sample was generated. Diversity metrics and Core microbiome per tissue The taxonomic profiles generated by AGAMEMNON at the species level were employed to calculate the microbiome richness of each sample. Microbiome richness was defined as the count of species with non-zero abundance in each sample. To compute the Shannon diversity index, the 'diversity' function from the R package vegan (version '2.6.4') was utilized. The core microbiome for each tissue was defined by considering all species present in at least 10% of the respective tissue’s samples. The 10% threshold was chosen to exclude species identified in a limited number of samples, that may be indicative of opportunistic infections, and to retain only those species consistently present in tissue samples. Given the inherent diversity and dynamic nature of the microbiome, a relatively low threshold (10%) was applied to avoid excluding too many species. For comparisons of the core microbiome across tissues, UpSet plots were generated using the R library UpSet (version 1.4.0). Normalization of taxonomic profiles To address biases, particularly due to differences in sequencing depth among samples, we employed Cumulative Sum Scaling (CSS) normalization in the taxonomic profiles. CSS, functioning as a median-like quantile normalization method, corrects for variations in sampling depth or library size. Unlike standard relative abundance normalization, which rescales all samples to a uniform total sum ( e.g. , 100%), CSS retains variability in total counts across samples. This normalization method adjusts samples based on a subset (quartile) of lower-abundance taxa that remain relatively constant and independent, mitigating the impact of high-abundance taxa that may dominate a study. For CSS normalization, the metagenomeSeq library (version 1.40.0) in R was employed. This normalization approach was also applied to normalize taxonomic profiles at the genus level and functional profiles generated by HumanN3. ML models 1vsAllOther27tissues To discern potential biological significance within the taxonomic profiles of samples, we endeavored to construct ML models for each tissue to discriminate it from all the other 27 tissues. Utilizing Gradient Boosting ML models known for their efficacy in classification tasks and resilience to imbalanced datasets, we trained stochastic gradient boosting machine (GBM) learning models. The models were implemented and fine-tuned using the R libraries GBM (version 2.1.8.1), Caret (version 6.0.94), and PRROC (version 1.3.1) (for calculating AUROC and AUPR values). The training and testing phases occurred on distinct, randomly selected, stratified sampling splits of 70% and 30% of the data, respectively, with a fixed random seed for reproducibility. Data CSS normalization was performed separately in each split. After normalisation of each split, only the total core species of the 28 tissues were retained for training the model as we were interested in the signatures present in the consistent microbiome within tissues. A two-fold cross-validation strategy was employed to create subsets of the training data and we conducted a basic grid search optimization of GBM parameters, with interaction depth (1, 2, or 3) and the number of trees (50, 100, or 150), with the aim of maximizing the area under the receiver operating characteristic (AUROC) of the final, fully trained model. The learning rate (shrinkage) remained constant at 0.1, and the minimum observations per node was fixed at 3. Final model performance metrics, including AUROC and area under the precision-recall curve (AUPR) values, were derived by applying the final model to the unseen 30% holdout test set. For calculating AUROC and AUPR the library PRROC was used. Given the stochastic nature of GBM models, 100 iterations of GBM models were created for each tissue model to calculate the mean AUROC and mean AUPR values, along with 95% confidence intervals. After that as the minority class for the different tissues were quite different across samples, we calculated a "normalized performance”, the relative AUPR, which was given by the following formula: Relative AUPR is equal to the ratio of AUPR(model)​ to AUPR(random). To expedite processes, the script for ML model creation utilized 10 cores, leveraging the R libraries “parallel” (version 4.2.2) and “doMC” (version 1.3.8). For the correlation of mean AUROC and relative mean AUPR of the different models with the size of the samples of each tissue, Spearman's rank correlation coefficient was applied using the cor.test from the R library stats (version 4.2.2). The same approach was used for creating the GBM models for the genera level. As features of the models were used the total core genera microbiome of the 28 tissues, which consisted of 738 genera. In silico Contamination For the in silico contamination approach, we strategically introduced 12 distinct in silico contaminants into our dataset, simulating various potential types of contamination that may occur within the GTEx framework. The contaminants were classified into two primary categories: high-volume contaminants, characterized by their high volume (high number of reads attributed to them) on a small subset of samples, positing that their detection would be possible if they affected a larger sample pool; and low-volume contaminants, which, although affecting a greater number of samples, did so at lower volumes (small number of reads attributed to them), thereby evading their detection. Consequently, we incorporated 6 high-volume and 6 low-volume contaminants into our dataset. The diverse types of contaminants is depicted in the Fig. 5 . We conducted 100 iterations, randomly selecting contaminated samples in each iteration using a custom Python script. The in silico contaminations were introduced post CSS normalization. Subsequently, Gradient Boosting Machine (GBM) models were constructed for the 11 tissues exhibiting a tissue-specific microbiome. These models aimed to distinguish each tissue from the other 10 tissues using similar methodology as described previously. For these models, the entire core microbiome of the aforementioned 11 tissues, comprising 1612 distinct microbial species, was retained for analysis. Additionally, 100 iterations of GBM models without any type of in silico contamination were generated to discriminate one tissue from the other 10 tissues, enabling comparison of results. Features’ importance scores of non-zero model features were estimated using the GBM and Caret packages. The percentage of contribution of a specific feature to the model's predictions was computed by dividing its importance score by the sum of all features’ importance scores for the given model. To assess the performance of contaminated versus uncontaminated models, we employed Wilcoxon's test on the AUROC and AURP values derived from the 100 iterations from contaminated and non-contaminated models across the 11 tissues, utilizing the R library stats (version 4.2.2). ML models 1vs7 tissues To focus exclusively on the signatures identified within the 8 contamination-resilient tissues (Blood, Testis, Colon tissue, Stomach tissue, Muscle, Bladder, Liver, and Heart) and analyze their distinctions, we reconstructed GBM models utilizing only samples derived from these tissues. These models utilized the concatenation of the core microbiome of these 8 tissues, consisting of 1597 species for the species models. The configurations of these models are similar to the ones described in the section ML models 1vsAllOther27tissues . To compute the mean AUROC and AUPR, along with their associated margins of error, 100 iterations of each model were executed. Features’ importance scores for the model features were estimated as previously described, employing the first iteration of the models. The identified signatures from these models underscored the differences in the microbial compositions inside these 8 tissues. The same methodology was applied to construct GBM models at the genus level. The features for these models comprised the total core microbiome (at genus level) present within the 8 tissues, encompassing 695 genera. ML models with Kaiju To validate the robustness of findings for the 8 tissue-specific signatures, we employed a different taxonomic tool to identify the taxonomic profiles of samples from the 8 contamination-resilient tissues. Utilizing the QCed reads (post-Atropos processed reads), we employed Kaiju (version 1.9.2) for taxonomic identification with 10 parallel threads allocated to each sample. Kaiju is a rapid and sensitive program designed for taxonomic classification of high-throughput sequencing reads. Each sequencing read is assigned to a taxon within the NCBI taxonomy by comparing it against a reference protein database. We utilized the microbial subset of the NCBI BLAST non-redundant protein database (nr), incorporating fungi and microbial eukaryotes as reference database. The database was downloaded on 03/30/2023 from https://bioinformatics-centre.github.io/kaiju/downloads.html . Reads were translated into amino acid sequences, which were then queried against the database using a modified backward search on a memory-efficient implementation of the Burrows-Wheeler transform, identifying maximum exact matches (MEMs) and optionally allowing mismatches. The tool kaiju2table was employed to convert Kaiju's output files into a summary table for the taxonomic rank of "species." Using these taxonomic profiles, GBM models were reconstructed for each tissue to discriminate it from all other 7 tissues, following a similar methodology as described previously. In these GBM models, the features used comprised the concatenation of the core microbiome of the 8 tissues based on the Kaiju tool, encompassing 1864 microbial species. To determine the mean AUROC and AUPR, along with their respective margins of error, 100 iterations of each model were executed. Feature importance scores for the model features were estimated in one iteration of the models, as described earlier. Functional Assignment using HumanN3 To characterize the functional profile of samples from the GTEx consortium, reads that were identified as microbial by AGAMEMNON were analyzed using HUMAnN3 (version 3.6.1) for profiling the abundance of microbial genes. HUMAnN3 is a method designed to efficiently and accurately conduct functional profiling in metagenomic or metatranscriptomic sequencing data. Due to constraints in computational resources and time, HUMAnN3 was executed with the option "--bypass-nucleotide-search" to skip all alignment steps before the translated search. For this translated search, the full UniRef90 database (version 201901b) served as the reference. Subsequently, the functional profiles of all samples were concatenated and normalized using CSS normalization. GBM models were then created following the approach described previously, utilizing the 1-vs-7 tissues strategy and incorporating only the core gene repertoire. This repertoire consisted of genes present in at least 10% of each tissue’s samples, amounting to 1,831 different microbial genes across the 8 tissues. For models integrating both functional and taxonomic profiles, the core gene repertoire (1,831 genes) and the core species microbiome (1,597 species) were utilized in the creation of GBM models. Each component was separately normalized using CSS normalization, and then the normalized gene and species taxonomy profiles were merged for the training and testing datasets. The GBM models were constructed using the same parameters as described earlier. To calculate the mean AUROC and AUPR, along with their margins of error, 100 iterations of each model were generated. Factors associated with tissue microbiome To investigate the influence of various traits (age, Body Mass Index - BMI, smoking status, drinking status, ancestry), as well as disease history (hypertension history, ischemic heart history, diabetes II history, diabetes I history, arthritis history, seizures history, schizophrenia history, rheumatoid arthritis history, liver disease history, dialysis treatment, depression history, COPD or CLRD history, cerebrovascular disease history, asthma history, alzheimer or dementia history), on the core microbiome of the 8 tissues with contamination-resilient tissue-specific microbiomes, GBM models were created for each trait and tissue combination with aim to discriminate in each tissue the samples with the specific trait from the samples without this trait. Tissues with fewer than 20 samples within each different group of specific trait/disease were excluded due to insufficient data for model creation. For continuous traits such as age and BMI, GBM models were constructed using 5-fold cross-validation and Root Mean Squared Error (RMSE) was used to select the optimal model with the smallest RMSE value. Mean Absolute Error (MAE) and R-squared were calculated using the "postResample" function from the Caret package. The mean RMSE, MAE, and R-squared were computed across 100 iterations of each model. For categorical traits (smoking, drinking, sex, and ancestry), as well as models for disease history, GBM models were developed similarly to tissue models but with 4-fold cross-validation to mitigate overfitting. Mean AUROC and mean AUPR were computed across 100 iterations of each model. For the ancestry trait, samples with ancestries other than "White" or "Black or African American" were removed due to limited representation. Validating significant ML models with living samples To validate the previously created 1-vs-7 tissues models in living tissues, we analyzed data of the 8 aforementioned tissues from project PRJEB4337, comprising RNA-seq data from living samples of the Swedish Biobank. From this project, only the samples from the 8 tissues were kept for testing the models, unfortunately this project did not contain samples from muscle tissue and blood. The RNA-seq data from this project were generated using a poly-A selection protocol and processed similarly to GTEx samples. Firstly, the fastq files were aligned to the same human genome as GTEx using STAR (version 2.7.10b). Then, the unmapped reads were isolated and quality controlled using Atropos with the same configurations as described in the GTEx analysis. Subsequently, the QCed reads were mapped to the microbial database using Puffaligner, and the microbial reads were remapped to the latest human genome using Bowtie2. Lastly, the taxonomical microbial profile of the 38 samples was assessed using Agamemnon. The identical pipeline was employed for both the GTEx and PRJEB4337 samples to ensure consistency. Next, the previously created 1-vs-7 tissues models (100 iterations per tissue) from the GTEx dataset were tested using the living samples from PRJEB4337. The AUROC and AUPR were calculated using the PRROC library (version 1.3.1). To verify that the observed high performance in the living tissues was attributed to the specific microbial signatures within these tissues, 100 iterations of GBM models were created using the GTEx data as described in the 1-vs-7 tissues models section, but the tissue labels of the samples in the training dataset were randomly shuffled to create random models without any specific tissue information. These random models were then tested with the living dataset to assess their performance. The performances in the living dataset of the 100 iterations of the original 1-vs-7 tissues models were compared with the performances of the 100 iterations of the random 1-vs-7 tissues models using a Wilcoxon test. This comparison was crucial for determining whether the performance observed in the living tissues was attributable to the specific microbial signatures within those tissues or was merely a result of chance. Lastly, samples from the PRJEB4337 project were also analyzed using the Kaiju tool to create their microbial profiles, following the same pipeline employed for GTEx samples. Subsequently, the 100 iterations of GBM models, created using Kaiju taxonomic profiles of GTEx data as described in the 'ML Models with Kaiju' section, were tested on the samples from the PRJEB4337 project. System configurations The analyses described in this study were conducted on a computational platform equipped with an Intel® Xeon(R) Gold 6226R CPU @ 2.90GHz and 376 GB of RAM. The operating system used was Ubuntu Linux version 20.04. Data processing and modeling were implemented using Python version 3.9.15 and R version 4.2.2. Declarations Availability of data and materials The taxonomic profiles of GTEx samples, generated using AGAMEMNON and Kaiju tools, along with their functional profiles generated via Humann3, are accessible in Zenodo (link https://zenodo.org/uploads/10980664 ). Additionally, taxonomic and functional profiles for the distinctive living dataset (NCBI ID PRJEB4337) are also available on the same link, along with its metadata R object used in the analysis. Although the full metadata for GTEx samples is restricted, a subset is available at https://www.gtexportal.org/home/downloads/adult-gtex/metadata; access to the complete metadata requires approval through a dbGaP application at https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000424.v8.p2. All programming scripts used to download and analyse the data of the GTEx and NCBI project PRJEB4337 as well as the pipelines for normalization, in silico contamination, ML models can be found at our GitHub repository https://github.com/dianalabgr/GTEx_microbiome_analysis . These scripts are designed to directly process the summarized count data available on Zenodo. However, analysis of microbial signatures concerning phenotypic traits and medical history necessitates access to the complete metadata. For the other analytical pursuits, the open access metadata suffices; nonetheless, it is essential to accommodate variations in column names between the open access and protected access metadata files. Acknowledgements The authors acknowledge the members of DIANA-lab for the very useful comments and ideas, and specifically Spyros Tatsoglou, PhD for his help with the procedures of accessing the data from GTEx Consortium. Authors’ contributions Anargyros Skoulakis: Conceptualization, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. Giorgos Skoufos: Conceptualization, Methodology, Writing – review & editing. Armen Ovsepian: Methodology, Writing – review & editing. Artemis G. Hatzigeorgiou: Conceptualization, Supervision, Funding acquisition, Writing – review & editing Funding This research has been co‐financed by the European Regional Development Fund of the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship and Innovation, and Greece 2.0, under the call “RESEARCH – CREATE – INNOVATE” (ID 16971), with project id: TAEDK-06179. Also, specifically for Anargyros Skoulakis, the research work was supported also by the Hellenic Foundation for Research and Innovation (HFRI) under the 5th Call for HFRI PhD Fellowships (Fellowship Number: 20480) Ethics approval and consent to participate Not applicable Consent for publication Not applicable Conflict of Interest The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. References Merino, N., Aronson, H. S., Bojanova, D. P., Feyhl-Buska, J., Wong, M. L., Zhang, S., & Giovannelli, D. (2019). Living at the extremes: extremophiles and the limits of life in a planetary context. Frontiers in microbiology, 10, 780. Ursell, L. K., Metcalf, J. L., Parfrey, L. W., & Knight, R. (2012). Defining the human microbiome. Nutrition reviews, 70(suppl_1), S38-S44. The Human Microbiome Project Consortium. A framework for human microbiome research. Nature 486, 215–221 (2012). https://doi.org/10.1038/nature11209 The Human Microbiome Project Consortium. Structure, function and diversity of the healthy human microbiome. Nature 486, 207–214 (2012). https://doi.org/10.1038/nature11234 The Integrative HMP (iHMP) Research Network Consortium. The Integrative Human Microbiome Project. Nature 569, 641–648 (2019). https://doi.org/10.1038/s41586-019-1238-8 Natalini, J.G., Singh, S. & Segal, L.N. The dynamic lung microbiome in health and disease. Nat Rev Microbiol 21, 222–235 (2023). https://doi.org/10.1038/s41579-022-00821-x Chen, X., Lu, Y., Chen, T., & Li, R. (2021). The female vaginal microbiome in health and bacterial vaginosis. Frontiers in cellular and infection microbiology, 11, 631972. Poore, G. D., Kopylova, E., Zhu, Q., Carpenter, C., Fraraccio, S., Wandro, S., … Knight, R. (2020). Microbiome analyses of blood and tissues suggest cancer diagnostic approach. Nature, 579(7800), 567–574. Dohlman, A. B., Mendoza, D. A., Ding, S., Gao, M., Dressman, H., Iliev, I. D., … Shen, X. (2021). The cancer microbiome atlas: a pan-cancer comparative analysis to distinguish tissue-resident microbiota from contaminants. Cell host & microbe, 29(2), 281–298. Narunsky-Haziza, L., Sepich-Poore, G. D., Livyatan, I., Asraf, O., Martino, C., Nejman, D., … Straussman, R. (2022). Pan-cancer analyses reveal cancer-type-specific fungal ecologies and bacteriome interactions. Cell, 185(20), 3789–3806. Nejman, D., Livyatan, I., Fuks, G., Gavert, N., Zwang, Y., Geller, L. T., … Straussman, R. (2020). The human tumor microbiome is composed of tumor type–specific intracellular bacteria. Science, 368(6494), 973–980. Riquelme, E., Zhang, Y., Zhang, L., Montiel, M., Zoltan, M., Dong, W., … McAllister, F. (2019). Tumor microbiome diversity and composition influence pancreatic cancer outcomes. Cell, 178(4), 795–806. Ghaddar, B., Biswas, A., Harris, C., Omary, M. B., Carpizo, D. R., Blaser, M. J., & De, S. (2022). Tumor microbiome links cellular programs and immunity in pancreatic cancer. Cancer Cell, 40(10), 1240–1253. Aykut, B., Pushalkar, S., Chen, R., Li, Q., Abengozar, R., Kim, J. I., … Miller, G. (2019). The fungal mycobiome promotes pancreatic oncogenesis via activation of MBL. Nature, 574(7777), 264–267. Mahmoudabadi, G., Consortium, T. S., & Quake, S. R. (2022). Single cell transcriptomics reveals the hidden microbiomes of human tissues. bioRxiv, 2022–10. Hieken, T. J., Chen, J., Hoskin, T. L., Walther-Antonio, M., Johnson, S., Ramaker, S., … Degnim, A. C. (2016). The microbiome of aseptically collected human breast tissue in benign and malignant disease. Scientific reports, 6(1), 30751. Lluch, J., Servant, F., Païssé, S., Valle, C., Valiere, S., Kuchly, C., … Lelouvier, B. (2015). The characterization of novel tissue microbiota using an optimized 16S metagenomic sequencing pipeline. PloS one, 10(11), e0142334. Lonsdale, J., Thomas, J., Salvatore, M., Phillips, R., Lo, E., Shad, S., … Moore, H. F. (2013). The genotype-tissue expression (GTEx) project. Nature genetics, 45(6), 580–585. Skoufos, G., Almodaresi, F., Zakeri, M., Paulson, J. N., Patro, R., Hatzigeorgiou, A. G., & Vlachos, I. S. (2022). AGAMEMNON: an Accurate metaGenomics And MEtatranscriptoMics quaNtificatiON analysis suite. Genome biology, 23(1), 39. Menzel, P., Ng, K. L., & Krogh, A. (2016). Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nature communications, 7(1), 11257. Paulson, J. N., Stine, O. C., Bravo, H. C., & Pop, M. (2013). Differential abundance analysis for microbial marker-gene surveys. Nature methods, 10(12), 1200–1202. Friedman, J. H. (2001). Greedy function approximation: a gradient boosting machine. Annals of statistics, 1189–1232. Almodaresi, F., Sarkar, H., Srivastava, A., & Patro, R. (2018). A space and time-efficient index for the compacted colored de Bruijn graph. Bioinformatics , 34 (13), i169-i177. Carlier, J. P., Bedora-Faure, M., K'ouas, G., Alauzet, C., & Mory, F. (2010). Proposal to unify Clostridium orbiscindens Winter et al. 1991 and Eubacterium plautii (Séguin 1928) Hofstad and Aasjord 1982, with description of Flavonifractor plautii gen. nov., comb. nov., and reassignment of Bacteroides capillosus to Pseudoflavonifractor capillosus gen. nov., comb. nov. International journal of systematic and evolutionary microbiology, 60(3), 585–590. Morita, T., Koike, H., Hagiwara, H., Ito, E., Machida, M., Sato, S., … Kitamoto, D.(2014). Genome and transcriptome analysis of the basidiomycetous yeast Pseudozyma antarctica producing extracellular glycolipids, mannosylerythritol lipids. PLoS One , 9 (2), e86490. Morita, T., Fukuoka, T., Imura, T., & Kitamoto, D. (2015). Mannosylerythritol lipids: production and applications. Journal of Oleo Science, 64(2), 133–141. Gonelli, A., Boccia, S., Boni, M., Pozzoli, A., Rizzo, C., Querzoli, P., … Di Luca,D. (2001). Human herpesvirus 7 is latent in gastric mucosa. Journal of medical virology,63(4), 277–283. Garcia-Montojo, M., Doucet-O’Hare, T., Henderson, L., & Nath, A. (2018). Human endogenous retrovirus-K (HML-2): a comprehensive review. Critical reviews in microbiology, 44(6), 715–738. Broderick, N. A., & Nagy, L. (2022). Bacteria may be in the liver, but the jury is still out. The Journal of Clinical Investigation, 132(8). Fagerberg, L., Hallström, B. M., Oksvold, P., Kampf, C., Djureinovic, D., Odeberg,J., … Uhlén, M. (2014). Analysis of the human tissue-specific expression by genome-wide integration of transcriptomics and antibody-based proteomics. Molecular & cellular proteomics, 13(2), 397–406. Beghini, F., McIver, L. J., Blanco-Míguez, A., Dubois, L., Asnicar, F., Maharjan,S., … Segata, N. (2021). Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3. elife, 10, e65088. Ciernikova, S., Sevcikova, A., Stevurkova, V., & Mego, M. (2022). Tumor microbiome–an integral part of the tumor microenvironment. Frontiers in Oncology, 12, 1063100. Rajilić-Stojanović, M., & De Vos, W. M. (2014). The first 1000 cultured species of the human gastrointestinal microbiota. FEMS microbiology reviews, 38(5), 996–1047. VanderHeyden, T. R., Yong, S. L., Breitschwerdt, E. B., Maggi, R. G., Mihalik, A. R., Parada, J. P., & Fimmel, C. J. (2012). Granulomatous hepatitis due to Bartonella henselae infection in an immunocompetent patient. BMC Infectious Diseases, 12, 1–7. Natalia, S., Luiza, H., Anna, S. C., Katarzyna, T., Schmidt, M., Łukasz, P., … Philips,A. (2023). Host Factors Associated with Gut Mycobiome Structure. MSystems, 8(2). Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., … Gingeras,T. R. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29(1), 15–21. Danecek, P., Bonfield, J. K., Liddle, J., Marshall, J., Ohan, V., Pollard, M. O.,… Li, H. (2021). Twelve years of SAMtools and BCFtools. Gigascience, 10(2), giab008. Didion, J. P., Martin, M., & Collins, F. S. (2017). Atropos: specific, sensitive, and speedy trimming of sequencing reads. PeerJ, 5, e3720. Almodaresi, F., Zakeri, M., & Patro, R. (2021). PuffAligner: a fast, efficient and accurate aligner based on the Pufferfish index. Bioinformatics, 37(22), 4048–4055. Gihawi, A., Ge, Y., Lu, J., Puiu, D., Xu, A., Cooper, C. S., … Salzberg, S. L. (2023).Major data analysis errors invalidate cancer microbiome findings. Mbio, 14(5), e01607-23. Langmead, B., Wilks, C., Antonescu, V., & Charles, R. (2019). Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics, 35(3), 421–432. Langmead, B., & Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nature methods, 9(4), 357–359. Additional Declarations No competing interests reported. Supplementary Files Supplementary.zip SupplementaryFiles.docx Cite Share Download PDF Status: Published Journal Publication published 25 Jan, 2025 Read the published version in Frontiers in Microbiology → Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4610181","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":321232259,"identity":"a0db06c0-b2ae-414d-b9fc-161b7d6d5e20","order_by":0,"name":"ANARGYROS SKOULAKIS","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA1UlEQVRIie2PMQrCQBBFx35lOxlR2CuMBMRCyFU2eAFLK/EAK7YRj2ATCaReEawSbFdstFdImcJC18LSrJ3gvu7Df8wfAI/nF2HQ0KARuA1y7KbAS2nPrELuCgBpm1wUoXLSVT4Ig8M6Ks8Egrf0Z4UKRVtlMMrMLcXnsN5yJWsUzkizEmXfFIlVJJ1qFLFgtL2XGAZxkVZOCjyH7ZjBRsLnmdsVyvfjXTfHKDbNbCAJ638RarS5XPfTkC+K9FhNhoJ36oa9wVcTXesWrr9pezwezz/xAGpASYN+hEIAAAAAAElFTkSuQmCC","orcid":"","institution":"DIANA-Lab, Department of Computer Science and Biomedical Informatics, University of Thessaly","correspondingAuthor":true,"prefix":"","firstName":"ANARGYROS","middleName":"","lastName":"SKOULAKIS","suffix":""},{"id":321232260,"identity":"d4008a3a-d1c4-4276-a4dc-f8adf00457c4","order_by":1,"name":"GIORGOS SKOUFOS","email":"","orcid":"","institution":"DIANA-Lab, Department of Computer Science and Biomedical Informatics, University of Thessaly","correspondingAuthor":false,"prefix":"","firstName":"GIORGOS","middleName":"","lastName":"SKOUFOS","suffix":""},{"id":321232261,"identity":"bf0c9229-35a5-477d-a061-0b373f491474","order_by":2,"name":"ARMEN OVSEPIAN","email":"","orcid":"","institution":"DIANA-Lab, Department of Computer Science and Biomedical Informatics, University of Thessaly","correspondingAuthor":false,"prefix":"","firstName":"ARMEN","middleName":"","lastName":"OVSEPIAN","suffix":""},{"id":321232262,"identity":"9de7de5e-97f4-43db-aecc-c34b22739de0","order_by":3,"name":"ARTEMIS G. HATZIGEORGIOU","email":"","orcid":"","institution":"DIANA-Lab, Department of Computer Science and Biomedical Informatics, University of Thessaly","correspondingAuthor":false,"prefix":"","firstName":"ARTEMIS","middleName":"G.","lastName":"HATZIGEORGIOU","suffix":""}],"badges":[],"createdAt":"2024-06-20 08:10:24","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-4610181/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4610181/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.3389/fmicb.2024.1512304","type":"published","date":"2025-01-26T00:00:00+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":60392516,"identity":"62c0697d-6652-4d19-840b-f13990f54ef8","added_by":"auto","created_at":"2024-07-16 09:17:12","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":169921,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFlowchart illustrating the employed methodology. \u003c/strong\u003eRNA-seq samples were quality controlled and pre-processed. Subsequently, microbial reads were assigned to specific microbial species utilizing AGAMEMNON\u003csup\u003e19\u003c/sup\u003e, a very accurate metagenomics and metatranscriptomics quantification analysis suite. The microbial taxonomic profiles of samples originating from diverse tissues were then utilized to discriminate the different tissues using Gradient Boosting ML models. 11 tissues presented tissue-specific microbial signatures, capable of discriminating these tissues against all the others. Following an\u003cem\u003e in-silico\u003c/em\u003e contamination analysis, eight tissues presented contamination-resilient microbial signatures, underscoring the potential existence of tissue-specific microbiomes in these tissues. Additionally, the models of these 8 tissues, trained in GTEx data, were further tested on living samples (NCBI project PRJEB4337), and 3 living tissues (heart, colon and liver) were discriminated based on their microbial signatures. Furthermore, these 8 tissues were further analysed to find the most important features of these signatures. An evaluation of the potential influence of various lifestyle factors and medical history was conducted on the identified signatures. Lastly, the presence of tissue-specific models was further validated using Kaiju\u003csup\u003e20\u003c/sup\u003e, a different metatranscriptomics taxonomic tool.\u003c/p\u003e","description":"","filename":"floatimage123.png","url":"https://assets-eu.researchsquare.com/files/rs-4610181/v1/d54a1e1249ea3f25e8b3355a.png"},{"id":60393295,"identity":"16cf8305-2544-4dc0-a1a2-bb51298c0dbf","added_by":"auto","created_at":"2024-07-16 09:25:12","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":614363,"visible":true,"origin":"","legend":"\u003cp\u003eA) The mean of reads (log\u003csub\u003e10 \u003c/sub\u003espace) at each step of the analysis, including raw reads, unmapped to human genome reads, QCed reads, microbial reads, non-human microbial reads and lastly classified microbial reads, B) Boxplots illustrating the microbial richness of the samples across 28 tissues (the darker shades of blue correspond to lower mean values, while brighter shades indicate higher mean values.), C) Boxplots depicting the microbial Shannon index of samples in each tissue, D) The distribution of the core microbiome per tissue in the four kingdoms of life, E) An UpSet plot illustrating the shared species present in the core microbiome across all tissues, comprising a total of 1,708 species.\u003c/p\u003e","description":"","filename":"floatimage219.png","url":"https://assets-eu.researchsquare.com/files/rs-4610181/v1/0671ccac0248442c01d7626f.png"},{"id":60392520,"identity":"98e1a29f-2962-4ed4-9619-0c4cdd55b201","added_by":"auto","created_at":"2024-07-16 09:17:12","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":566494,"visible":true,"origin":"","legend":"\u003cp\u003eA) The performance of the models (AUROC and relative AUPR) for all tissues in the models 1TissuevsAll27OtherTissues. It is evident that eleven tissues presented high AUROC and AUPR values. B) A bubble plot showing the mean importance score of each contamination type in the 11 tissues models. The contamination 8 plays an important role in in small intestine, brain and salivary gland models. C) The performance of the models (AUROC and relative AUPR) for the tissues that presented tissue-specific contamination-resilient microbiome. All eight tissues presented high AUROC and AUPR values in the 1vs7OtherTissues models. D-K) Boxplots of the normalized abundances of the most important features of each model, in the headline of each boxplot is written the species name and in parenthesis the tissue for which this species is the most important feature. L) The performance of the models (AUROC and relative AUPR) for the tissues that presented tissue-specific contamination-resilient microbiome but with the taxonomic profiles produced by Kaiju. All eight tissues presented high AUROC and AUPR values in the 1vs7OtherTissues models even with the employment of a total different taxonomic tool.\u003c/p\u003e","description":"","filename":"floatimage314.png","url":"https://assets-eu.researchsquare.com/files/rs-4610181/v1/400e49694ba891e68fbd3596.png"},{"id":60392523,"identity":"8798aa4f-b7a0-40fd-b485-0932125f76c1","added_by":"auto","created_at":"2024-07-16 09:17:12","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":364345,"visible":true,"origin":"","legend":"\u003cp\u003eA) The performance of the 8 models created from GTEx samples tested on the samples from living tissues. Colon tissue, heart and Liver presented a high AUROC and relative AUPR in the living samples, showing that the microbial signatures of these tissues are present during the lifetime of individuals too. B-D) The 5 most abundant microbial species on the living samples of the colon tissue, heart, and liver. The species that are common among tissues are depicted with same color.\u003c/p\u003e","description":"","filename":"floatimage410.png","url":"https://assets-eu.researchsquare.com/files/rs-4610181/v1/d4daefb52b36643982d2c52e.png"},{"id":60392517,"identity":"96a90a0b-e55d-45ec-8b0a-6acdcb9e9493","added_by":"auto","created_at":"2024-07-16 09:17:12","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":265664,"visible":true,"origin":"","legend":"\u003cp\u003eThe different types of in silico contaminants. In silico contaminants can be classified into two main categories: high abundance contamination and low abundance contamination. The primary routes of contamination in the GTEx original analysis protocol include the collection site, the DNA isolation step (DNA/RNA isolation procedure), and the expression analysis step (RNA sequencing procedure). At each step, we introduced different types of in silico contaminants to mitigate possible contamination scenarios.\u003c/p\u003e","description":"","filename":"floatimage511.png","url":"https://assets-eu.researchsquare.com/files/rs-4610181/v1/696565ebf97cf0bb4f85ecfe.png"},{"id":75444029,"identity":"0a62c52a-0544-44dd-ba05-dbe966b7f221","added_by":"auto","created_at":"2025-02-04 16:11:37","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":3361635,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4610181/v1/d6d2a759-464b-4be9-9339-3bb4447ad255.pdf"},{"id":60392519,"identity":"2c38de03-0fef-4192-915c-fa269d75e206","added_by":"auto","created_at":"2024-07-16 09:17:12","extension":"zip","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":2373465,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementary.zip","url":"https://assets-eu.researchsquare.com/files/rs-4610181/v1/d77ee2b15c34ebd6cdee426d.zip"},{"id":60392522,"identity":"b2cde811-a373-42d9-ba6e-4974fa44ba41","added_by":"auto","created_at":"2024-07-16 09:17:12","extension":"docx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":17886,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryFiles.docx","url":"https://assets-eu.researchsquare.com/files/rs-4610181/v1/88d3732a16ef5abea5b409aa.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Machine learning models reveal microbial signatures in healthy human tissues, challenging the sterility of human organs","fulltext":[{"header":"Background","content":"\u003cp\u003eMicrobes exhibit a remarkable spectrum of functions and capabilities that enable them to colonize diverse and extreme habitats. The human body serves as a host for a vast array of microbes with varied functionalities\u003csup\u003e1,2\u003c/sup\u003e. Our understanding regarding interactions between humans and microbes has evolved through the progression of sequencing technologies. Pioneering this transformative trajectory was the NIH Human Microbiome Project (HMP), which revealed the complex interplay between human and microbial cells\u003csup\u003e3,4,5\u003c/sup\u003e. Recent advancements have prompted a reevaluation of long-standing beliefs, showcased by the dismissal of the traditional notion that healthy human lung and vagina, are sterile environments.\u003csup\u003e6,7\u003c/sup\u003e\u003c/p\u003e \u003cp\u003eThe predominant focus of microbiome studies has been on characterizing microbial communities in easily accessible sampling sites, including the human intestinal tract, the oral cavity, and the skin. Internal tissues, less accessible by conventional sampling methods, have historically been presumed to be devoid of microorganisms due to protective layers of epithelial and endothelial tissues. However, this assumption has been challenged, particularly in settings such as cancer, where various studies have revealed the presence of tumor-specific microbial sequences within human tumors, offering a potential new avenue for cancer diagnosis\u003csup\u003e8\u0026ndash;14\u003c/sup\u003e. As the exploration of microbes inside internal tissues has primarily been conducted in the context of pathological conditions like cancer, fundamental questions regarding healthy status persist. Are there microbial residents within human tissues engaged in a commensal relationship with human cells? What is the role of their presence? Can these microbial signatures trigger an immunological response or serve as diagnostic markers for organ health?\u003c/p\u003e \u003cp\u003eGita Mahmoudabadi \u003cem\u003eet al.\u003c/em\u003e observed a notable presence of bacterial genera in tumors that is also detectable and similar in adjacent tumor-free tissues\u003csup\u003e15\u003c/sup\u003e. This suggests that tumor microbiomes may partially originate from neighboring normal tissues, indicating that, under healthy conditions, internal tissues host a microbiome. Hieken \u003cem\u003eet al\u003c/em\u003e also detected microbial communities in aseptically collected human breast tissues in benign and malignant conditions\u003csup\u003e16\u003c/sup\u003e. Moreover, investigations in germ-free mice have unveiled microbiomes within various organs, including the brain, muscle, adipose tissue, liver, and heart. This challenges the conventional notion that the internal organs of mammals are devoid of microbial presence.\u003csup\u003e17\u003c/sup\u003e\u003c/p\u003e \u003cp\u003eTo elucidate the microbiome inside healthy human tissues, we conducted a re-analysis of RNA-seq sequencing data derived from Genotype-Tissue Expression (GTEx) consortium\u003csup\u003e18\u003c/sup\u003e. Originally designed to explore variations in gene expression within healthy human tissues, GTEx's extensive repository contains samples exclusively derived from healthy specimens. Acknowledging certain limitations within the GTEx analysis pipeline, such as the employment of a poly-A selection protocol and the reliance on post-mortem biopsy samples, it remains noteworthy that the GTEx consortium provides the most extensive and analytically robust dataset about RNA expression within healthy human tissues. To date, GTEx sequencing data have not been explored from a microbiome perspective. Here, we present what, to our knowledge, stands as the inaugural and comprehensive healthy human tissue microbiome dataset. Leveraging ML models, we identified microbial signatures capable of discriminating among various tissue types and investigated the potential influence of phenotypic traits (e.g., age, sex, BMI) on these signatures. The presence of tissue-specific microbial signatures in certain tissues would suggest that these organs are not devoid of microorganisms even in healthy conditions, thus challenging the long-standing dogma of the sterility of internal tissues.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eLow-biomass microbial RNA is detected in GTEx samples\u003c/h2\u003e \u003cp\u003eA total of 13,871 RNA-seq sequencing samples originating from 28 distinct tissues (e.g., Lung, Liver, Pituitary, Blood Vessel, Thyroid, Skin, Salivary Gland, Esophagus, Heart, Muscle, Pancreas, Adipose Tissue, Vagina, Blood, Ovary, Spleen, Prostate, Adrenal Gland, Nerve, Stomach tissue, Colon tissue, Testis, Brain, Breast, Uterus, Small Intestine tissue, Kidney, Bladder tissue) within the GTEx consortium have been utilized. All samples were designated as non-diseased and free of pathology, as the specimens were reviewed by a panel of 2\u0026ndash;3 pathologists and any specimen found with an incidental finding had been systematically excluded.\u003c/p\u003e \u003cp\u003eThe workflow of the study is shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e. Following the removal of reads of human origin, the remaining reads were aligned against an extensive database comprising 4,034 bacterial, 489 archaeal, 11,259 viral, and 81 fungal entities. The microbial composition of each sample was determined utilizing AGAMEMNON\u003csup\u003e19\u003c/sup\u003e. AGAMEMNON provides precise genus, species, and strain abundances through an efficient indexing scheme for rapid pattern matching, facilitating analysis of extensive datasets using common computational resources. The abundance estimation employs expectation maximization algorithm and targets maximizing the likelihood of the observed reads by gradually altering the abundance value associated to different taxa. On average, out of ~\u0026thinsp;93\u0026times;10\u003csup\u003e6\u003c/sup\u003e raw reads per sample, ~\u0026thinsp;2\u0026times;10\u003csup\u003e6\u003c/sup\u003e reads (2.35% of the raw reads) did not align to the human genome. After quality control (QC) and filtering of the non-human reads, ~\u0026thinsp;3.5\u0026times;10\u003csup\u003e5\u003c/sup\u003e reads (0.38% of the raw reads) were preserved. Employing Pufferfish\u003csup\u003e21\u003c/sup\u003e, ~\u0026thinsp;4.3\u0026times;10\u003csup\u003e4\u003c/sup\u003e (0.046% of the raw reads) were attributed to bacteria, archaea, fungi, or viruses. To eliminate any residual human-associated sequences, a realignment of microbial reads to the most recent human genome assembly (T2T-CHM13v2.0) was executed, resulting in ~\u0026thinsp;4.2\u0026times;10\u003csup\u003e4\u003c/sup\u003e reads (0.045% of raw reads) remaining. Finally, using AGAMEMNON, all the remaining reads (4.2\u0026times;10\u003csup\u003e4\u003c/sup\u003e reads, 0.045% of raw reads) were successfully classified into microbial species (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea). The reported read counts are the mean values across all samples at each step of the analysis.\u003c/p\u003e \u003cp\u003eThe bladder exhibited the most pronounced abundance of microbial reads, whereas the adrenal gland, brain, heart, muscle, salivary gland, and testis demonstrated the lowest microbial read counts, as illustrated in \u003cb\u003eSupplementary Fig.\u0026nbsp;1\u003c/b\u003e and in \u003cb\u003eSupplementary File 1\u003c/b\u003e. In each of the 28 tissues examined, a discernible presence of low-biomass microbial RNA was identified. The microbial communities detected within these tissues may either genuinely inhabit the respective tissues or be a result of potential contamination. As in the case of microbial reads, the bladder exhibited the highest microbial species richness in its microbial community, whereas the brain, heart, muscle, salivary gland, and testis displayed the lowest richness, as depicted in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb. Notably, the Shannon index of the bladder was markedly lower compared to the rest of the samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec\u003cb\u003e)\u003c/b\u003e, which arises from the fact that bladder samples demonstrated elevated counts for a few selected species, and comparatively lower counts for the remaining species..\u003c/p\u003e \u003cp\u003eTo identify species that are consistently present in each tissue and are not sporadic opportunistic pathogens, we retained species present in at least 10% of samples from each tissue. The retained species are hereafter referred to as the core microbiome of the tissue. Species meeting this criterion for at least one tissue were included in subsequent analyses, resulting in a total of 1,708 different species across all four kingdoms. As expected, the tissue with the smallest core microbiome was the brain, consisting of 720 species, while the bladder exhibited the most diverse core microbiome, comprising 1,344 species (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed and \u003cb\u003eSupplementary File 2\u003c/b\u003e for details). Additionally we noticed that across all tissues, bacteria constituted the predominant component, making up over 95% of the core microbiome for each tissue. The core microbiome of each tissue for all four kingdoms, namely bacteria, fungi, viruses and archaea, is documented in \u003cb\u003eSupplementary Files 3 to 6\u003c/b\u003e respectively.\u003c/p\u003e \u003cp\u003eNotably, a total of 529 species, accounting for ~\u0026thinsp;31% of the overall core microbiome, were shared across all tissues, showing that tissues share in some extent a common microbiome (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ee). Among fungi, the majority (24 out of 37 species) were present in all tissues, as indicated in \u003cb\u003eSupplementary Fig.\u0026nbsp;1\u003c/b\u003e and \u003cb\u003eSupplementary File 7\u003c/b\u003e. However, for bacteria, only 504 out of 1,297 species (38.86%) and for viruses, only \u003cem\u003eGeobacillus\u003c/em\u003e virus E2 out of 13 viruses were present in all tissues, as illustrated in \u003cb\u003eSupplementary Fig.\u0026nbsp;2, 3\u003c/b\u003e and \u003cb\u003eSupplementary Files 8, 9\u003c/b\u003e, respectively. Archaea were exclusively identified in specific tissues, including the bladder and kidney (species \u003cem\u003eMethanocaldococcus jannaschii\u003c/em\u003e), the colon and small intestine tissues (species \u003cem\u003eMethanosarcina sp. WH1\u003c/em\u003e), and the spleen (species \u003cem\u003eMethanocaldococcus jannaschii\u003c/em\u003e) (\u003cb\u003eSupplementary File 6)\u003c/b\u003e.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003e8 tissues harbor specific microbial signatures at species level\u003c/h2\u003e \u003cp\u003eWhile solely the identification of microbiomes within tissues doesn't warrant the existence of microbial communities within these tissues, the presence of discernible biological footprints within these communities would imply a non-coincidental occurrence of microbes inside these tissues. To explore this, Machine Learning models were employed to detect potential microbial footprints for each tissue. By using Cum-Sum Scaling method\u003csup\u003e21\u003c/sup\u003e to normalize microbial profiles across samples, stochastic Gradient Boosting Machine Learning models (GBM models)\u003csup\u003e22\u003c/sup\u003e were trained to distinguish individual tissues from all other tissues. Notably, 11 out of the 28 models (models for Brain, Small Intestine tissue, Liver, Bladder, Muscle, Heart, Salivary Gland, Stomach tissue, Colon tissue, Testis, and Blood) exhibited robust performance (mean AUROC\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;0.70 and relative mean AUPR\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;1.4) discriminating each specific tissue type from the collective representation of all other tissues, as illustrated in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ea and in \u003cb\u003eSupplementary File 10\u003c/b\u003e (\u003cspan refid=\"Sec10\" class=\"InternalRef\"\u003eMethods\u003c/span\u003e section ML models 1vs27AllOtherTissues). There was no significant correlation between the sample size and mean AUROC performance (rho\u0026thinsp;=\u0026thinsp;0.207, p-value\u0026thinsp;=\u0026thinsp;0.28) and the normalized mean AUPR performance (rho=-0.152, p-value\u0026thinsp;=\u0026thinsp;0.43).\u003c/p\u003e \u003cp\u003eThe identified microbial signatures based on species abundances within these 11 tissues were further analyzed to investigate if the microbial presence could be attributed to sample contamination. Given that the primary objective of the GTEx initiative was to investigate gene expression, the samples were not processed under sterile conditions. Due to the nature of GTEx experiment protocol, it was inefficient to employ an \u003cem\u003ein-silico\u003c/em\u003e decontamination approach and to overcome this inefficiency, an alternative strategy was implemented. A randomized \u003cem\u003ein silico\u003c/em\u003e contamination approach was used to evaluate the potential impact of systematic contamination on the performance of the models. To scrutinize whether possible contamination could drive the performance of ML models exhibiting high performance, 12 pseudo-contaminants were strategically introduced to the data. These pseudo-contaminants aimed to simulate diverse scenarios of contamination throughout the analytical process, encompassing potential contamination scenarios from the initial sampling to the sequencing procedure (\u003cspan refid=\"Sec10\" class=\"InternalRef\"\u003eMethods\u003c/span\u003e section \u0026ldquo;In silico Contamination\u0026rdquo;).\u003c/p\u003e \u003cp\u003eWe retrained the ML models of the 11 tissues presenting high discriminatory power (mean AUROC\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;0.70 and relative mean AUPR\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;1.4) and compared their performances with and without \u003cem\u003ein silico\u003c/em\u003e contamination (\u003cb\u003eSupplementary file 11\u003c/b\u003e). Only the models of Salivary Gland, Brain, and Small Intestine had significantly higher performance (both mean AUROC values and mean AUPR values) in the contaminated iterations, with in silico contaminants playing an important role in the models\u0026rsquo; performance. On the other hand, for the rest of the models in the contaminated iterations, models did not rely on the inserted pseudo-contaminants (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eb), Therefore, the models for the Salivary Gland, Brain, and Small Intestine tissues were deemed unreliable as their high performance could be the result of contamination. The contamination type, that could potentially influence the performance of these 3 models, is a low volume contamination (contamination type 8) which imitates a contaminant that is systematically present in the samples of a collection site center (\u003cspan refid=\"Sec10\" class=\"InternalRef\"\u003eMethods\u003c/span\u003e section \u0026ldquo;In silico contamination\u0026rdquo;). Nevertheless, in the rest 8 tissues (Blood, Testis, Colon tissue, Stomach tissue, Muscle, Bladder, Liver, and Heart), the tissue-specific microbial signatures exhibited resilience against potential contamination. The ability of GBM models to discriminate these tissues solely based on their microbial profiles independently of any potential contaminations, underscores the presence of distinctive microbial signatures specific to each of these tissues. To investigate microbial differences among these 8 tissues, GBM models were recreated using data exclusively from these contamination-resilient tissues, distinguishing each tissue from the combined representation of all the other 7 tissues (\u003cspan refid=\"Sec10\" class=\"InternalRef\"\u003eMethods\u003c/span\u003e section \u0026ldquo;ML models 1vs7Tissues\u0026rdquo;). Remarkably, all 8 tissues demonstrated robust performance in terms of AUROC and relative AUPR (mean AUROC\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;0.70 and relative mean AUPR\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;1.4), as illustrated in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003ec and in \u003cb\u003eSupplementary File 12\u003c/b\u003e.\u003c/p\u003e \u003cp\u003eCollectively, these results suggest that each of these 8 tissues harbors a distinct and unique microbiome. To assess the biological relevance of these microbial signatures, a detailed examination of the most important features from the 8 models was conducted. Boxplots of the abundance of the most important features across the different tissues are depicted in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD-\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eK. \u003cb\u003eSupplementary File 13\u003c/b\u003e contains the full catalog of features importance scores for each tissue. In the colon tissue model, the 2 most predominant features were the bacterium \u003cem\u003eEnterocloster bolteae\u003c/em\u003e, formerly known as \u003cem\u003eClostridium bolteae\u003c/em\u003e, a recognized constituent of human feces and \u003cem\u003eFlavonifractor plautii\u003c/em\u003e (formerly \u003cem\u003eEubacterium plautii\u003c/em\u003e), which has been isolated from human feces, blood, intra-abdominal pus, and infected soft tissues in humans\u003csup\u003e24\u003c/sup\u003e (Figs.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD, \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eE). In the case of Heart, the most crucial feature was the fungus \u003cem\u003eSporisorium graminicola\u003c/em\u003e and for the Blood model, the 2 most important features were \u003cem\u003eSchizosaccharomyces pombe\u003c/em\u003e and \u003cem\u003eSporisorium graminicola\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eF, \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eG respectively). These two fungal species, discussed in the study by Hidayat \u003cem\u003eet al.\u003c/em\u003e, are capable of synthesizing mannosylerythritol lipids (MELs)\u003csup\u003e25\u003c/sup\u003e. MELs belong to the glycolipid class of biosurfactants and are known for their outstanding interfacial and biochemical characteristics, as highlighted by Morita \u003cem\u003eet al.\u003c/em\u003e\u003csup\u003e26\u003c/sup\u003e. For the Muscle model, the overall importance of features was relatively subdued, and many features contributed equally lightly in the model\u0026rsquo;s discriminatory ability. In the Stomach model, the pivotal feature was the virus \u003cem\u003eHuman betaherpesvirus 7\u003c/em\u003e, previously identified as an inhabitant of gastric mucosa\u003csup\u003e27\u003c/sup\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eH). For the Testis model, the critical feature was the \u003cem\u003eHuman endogenous retrovirus K\u003c/em\u003e (HERV-K) (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eI). The human genome harbors numerous copies of HERV-K, many of which retain intact open reading frames (ORFs). These ORFs are capable of being transcribed and translated, particularly during early embryonic development and in cancerous conditions.\u003csup\u003e28\u003c/sup\u003e For the Bladder model, the foremost feature was the bacterium \u003cem\u003eCaldimonas thermodepolymerans\u003c/em\u003e, an underexplored microorganism phylogenetically proximate to the \u003cem\u003eComamonadaceae\u003c/em\u003e group that was systematically less abundant in the samples of bladder (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eJ). Lastly, in the Liver model, the primary feature was \u003cem\u003eBartonella machadoae\u003c/em\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eK). \u003cem\u003eBartonella machadoea\u003c/em\u003e is a bacteria belonging to \u003cem\u003eProteobacteria genus\u003c/em\u003e, and it was recently reported that the liver is inhabited mainly by proteobacteria by a gut-liver-specific axis\u003csup\u003e29\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eFor additional validation, the samples from the aforementioned eight tissues underwent analysis using Kaiju, a different computational method for microbial taxonomic profiling. Kaiju assigns each sequencing read to a taxon in the NCBI taxonomy by comparing it to a reference protein database. The used reference database was the microbial subset of the NCBI BLAST non-redundant protein database, encompassing bacteria, archaea, virus, fungi and microbial eukaryotes. Subsequently, GBM models for these 8 tissues were reconstructed using the taxonomic profiles generated by Kaiju. Remarkably, the performance of these models exhibited a high performance (mean AUROC\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;0.70 and relative mean AUPR\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;1.4) as depicted in Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eL (\u003cb\u003eSupplementary file 14\u003c/b\u003e), indicating that regardless of the bioinformatics method analysis, these tissues present a distinctive microbial signature. The most important features of the microbial signatures found by Kaiju (\u003cb\u003eSupplementary file 15\u003c/b\u003e) were different compared to the microbial profiles derived by AGAMEMNON. Employing different approaches for read classification, the taxonomic profiles and respectively the most important features of the models are differing, mainly due to the different reference databases that each tool uses for the classification of the reads and to the algorithm used in assigning the reads to each organism.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eHeart, Colon tissue, and Liver preserve the tissue-specific microbial signatures in the living state\u003c/h2\u003e \u003cp\u003eAs the specimens sourced from the GTEx consortium originate from post-mortem biopsies, an investigation was conducted to examine whether the identified signatures for the eight tissues resulted from post-mortem microbial colonization or were present during the subjects' lifetime. To address this, the eight models, comparing one tissue against the seven others, were subjected to testing using data derived from an entirely distinct project (NCBI ID PRJEB4337). This dataset sourced from living tissues of subjects within the Swedish Biobank\u003csup\u003e30\u003c/sup\u003e. This project did not contain samples from muscle tissue and blood, so only 6 out of the 8 models (heart, colon, stomach, liver, bladder, testis models) were tested with data from living tissues. The data were processed utilizing the same pipeline as applied to GTEx data, and normalization was separately carried out on this specific distinct dataset using the CSS normalization method.\u003c/p\u003e \u003cp\u003eFor each tissue, all 100 iterations of the 1vs7Tissues models (that were created before) were tested again using the living dataset. The obtained results were then compared with the performance of models, generated with randomly assigned tissue labels to ascertain the significant contribution of tissue labeling to model performance (\u003cb\u003eSupplementary file 16\u003c/b\u003e). Remarkably, 5 out of 6 tissues (only the model for bladder tissue did not) exhibit a statistically significant difference in AUROC and AUPR in the living dataset (p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05) between the true tissue labeled model and the random tissue labeled model.\u003c/p\u003e \u003cp\u003eAmong the six models, the models for Heart, Colon tissue, and Liver demonstrated robust AUROC and relative AUPR in the living dataset (AUROC\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;0.7 and relative AUPR\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;1.4) (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003ea, \u003cb\u003eSupplementary file 17\u003c/b\u003e). This suggests that the microbial signatures identified in these tissues (Heart, Colon tissue, and Liver) also exist in the healthy living state. The 5 most abundant species in the living dataset for the 3 tissues are depicted in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB-D (refer for frequencies of all microbiomes inside living tissues of the project PRJEB4337 in \u003cb\u003eSupplemetary file 18\u003c/b\u003e).\u003c/p\u003e \u003cp\u003eThe lack of high performance for the remaining three tissues may be attributed to the likelihood that the microbial signature of these tissues undergoes significant alterations during the event of death. The high performance of the Heart, Colon tissue, and Liver models was not observed in the models created with the taxonomic profiles derived from Kaiju (\u003cb\u003eSupplementary File 19\u003c/b\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003eLifestyle factors and medical history do not influence tissues\u0026rsquo; microbiome\u003c/h2\u003e \u003cp\u003eAdditionally, The putative impact of lifestyle factors on the composition of tissue-specific microbiomes across the eight human tissues was examined. The following variables were studied : Age, Body Mass Index (BMI), Sex, Alcohol consumption, Smoking status, Ancestry). To elucidate potential associations, GBM models were constructed for each tissue to predict the specific trait. However, despite rigorous computational analyses, the findings reveal a lack of robust predictive capability of the tissue microbiomes for these traits across all examined tissues. Detailed performance metrics for each trait within each tissue are provided in \u003cb\u003eSupplementary File 20\u003c/b\u003e.\u003c/p\u003e \u003cp\u003eIn addition, each tissue underwent a thorough examination to discern whether its microbial profile harbored predictive potential for various medical conditions, encompassing Alzheimer's or Dementia, Asthma, Cerebrovascular Disease, COPD-CLRD, Depression, Diabetes Type I and II, Dialysis Treatment (Renal Failure), Rheumatoid Arthritis, Hypertension, Ischemic Heart Disease, Liver disease, Schizophrenia, Seizures, and Arthritis (as a comprehensive category inclusive of various types of arthritis). Only in the case of dialysis treatment (renal failure) could discernment be achieved between heart tissues of subjects with renal failure and heart tissues of subject without renal failure. However, for the rest of the diseases, discrimination of subjects\u0026rsquo; diagnosis across the tissues was not feasible. Detailed performance metrics for each disease in each of the eight tissues are presented in the \u003cb\u003eSupplementary file 21\u003c/b\u003e. Overall, the analysis of tissue microbiomes suggests a lack of inherent microbial signatures indicative of lifestyle factors or systemic disease status.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec7\" class=\"Section2\"\u003e \u003ch2\u003eThe utmost significance in unraveling tissue-specific microbial signatures lies within the species-level analysis\u003c/h2\u003e \u003cp\u003eTo elucidate the optimal taxonomic level for investigating microbial signatures across the eight vital tissues, we constructed Gradient Boosting Machine (GBM) models utilizing microbial profiles at the genus level. AGAMEMNON supports the aggregation of microbial counts at the taxonomic level of choice. Employing microbial profiles at the genus level, we developed GBM models for the aforementioned eight tissues. GBM models based on genera exhibited slightly inferior performance compared to those derived from species-level data (refer to \u003cb\u003eSupplementary File 22\u003c/b\u003e). Although the distinctions are subtle, the majority of the eight models demonstrated enhanced performance when analyzed at the species level; only, the bladder model exhibited superior mean AUROC and mean AUPR at the genera level.\u003c/p\u003e \u003cp\u003eTo investigate if any other of the twenty tissues, that did not present microbial footprints at species level analysis, manifested distinct microbial signatures at the genera level, we generated GBM models utilizing taxonomic profiles at the genera level for each tissue. Only the models from the same tissues as species level models presented high performance (AUROC\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;0.7 and relative AUPR\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;1.4), most of the times slightly diminished when compared to their species-level counterparts. Only the brain tissue models did not present high performance at the genera level, but as previously shown, the species-level brain model is possibly a result of contamination (refer to \u003cb\u003eSupplementary file 23\u003c/b\u003e).\u003c/p\u003e \u003cp\u003eIn light of the known redundancy in the functional repertoire (i.e. genes) across diverse microorganisms, we endeavored to investigate whether tissue-specific microbial signatures predominantly arise from a tissue-specific microbial gene repertoire rather than a specific microbial taxonomy. To delineate the functional profile of samples, microbial reads identified by AGAMEMNON were subjected to analysis using Humann3\u003csup\u003e31\u003c/sup\u003e, specifically employing translational searches against the Uniref90 database. Subsequently, GBM models for the aforementioned eight tissues were reconstructed utilizing the CSS-normalized functional profiles of the samples (refer to \u003cb\u003eSupplementary File 24\u003c/b\u003e). Remarkably, only the Blood, Colon, and Bladder tissues exhibited high performance, with a mean AUROC\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;0.70 and a relative mean AUPR\u0026thinsp;\u0026gt;\u0026thinsp;=\u0026thinsp;1.4. In an attempt to integrate both functional and species-level information, we reconstructed GBM models using both the functional and the taxonomic profiles of the samples. However, even with this comprehensive approach, a very slight enhanced performance was attained and not in all models (see \u003cb\u003eSupplementary File 25\u003c/b\u003e), indicating that in our scientific setup the information of functional repertoire contributes very slightly in discriminating tissues. It is crucial to note that the low performance of the functional models, could most probably be attributed to the technical characteristics of the GTEx RNA-seq data. The data, generated using a poly-A selection sequencing protocol, inherently filtered out a significant proportion of microbial RNAs.\u003c/p\u003e \u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eIn each of the 28 examined tissues, a modest yet substantiated presence of transcriptionally active microbial communities was observed. Contrary to traditional notions of sterility of internal tissues, this analysis, alongside analogous investigations in other mammalian species, suggests that tissues potentially harbor a low-biomass microbiome not only in the context of disease but also in normal state. Established knowledge in the field of cancer underscores the presence of a tumor microbiome, with mounting evidence indicating its pivotal role in tumor development and influencing treatment outcomes\u003csup\u003e32\u003c/sup\u003e. Extending beyond pathological conditions, our analysis suggests that certain human internal tissues consistently maintain a low-biomass microbiome. This phenomenon prompts speculation that the microbiome may serve as a functional reservoir contributing to tissue well-being or act as a regulator of the immune system.\u003c/p\u003e \u003cp\u003eThe robust discriminatory performance demonstrated by the models for eight tissues, relying solely on microbial profiles of samples, strongly suggest that the observed microbiome is not a result of random chance. Despite the impracticality of conducting \u003cem\u003ein silico\u003c/em\u003e decontamination due to the uniform utilization of the same RNA quantity for each sample in the processing of GTEx samples and the absence of negative controls, a reverse analysis \u0026mdash;implementing \u003cem\u003ein silico\u003c/em\u003e contamination\u0026mdash; illustrated the robustness of microbial signatures in eight tissues against various potential contaminations. It is imperative to clarify that the \u003cem\u003ein silico\u003c/em\u003e contamination approach does not seek to substitute a laboratory analysis conducted under the most stringent sterile conditions, nor is it exhaustive in considering all conceivable contaminations. Rather, it serves as a validation step to evaluate whether contamination occurred at different stages of the analysis and to assess its potential impact on models performance. In total, eight tissues demonstrated a tissue-specific microbiome that exhibited resilience against various potential contaminations, thereby indicating the presence of a distinct microbiome in each tissue. Interestingly, this tissue-specific microbiome primarily diverges not in terms of species but predominantly in the composition of these species. The microbial signatures for heart, colon and liver tissues seem to be distinctive also in an independent dataset from healthy living tissues, supporting the presence of microbiomes inside these tissues also in a living healthy state.\u003c/p\u003e \u003cp\u003eIn colon tissue, the most pivotal component of its microbial signature is \u003cem\u003eFlavonifractor plautii\u003c/em\u003e, a prevalent bacterium in the human gastrointestinal tract recognized for its notable butyrate production\u003csup\u003e33\u003c/sup\u003e. In liver tissue, the primary microbial entity of significance is \u003cem\u003eBartonella machadoae\u003c/em\u003e. \u003cem\u003eBartonella\u003c/em\u003e species have been associated with compromised liver function\u003csup\u003e34\u003c/sup\u003e. Lastly, the presence of the phytopathogenic fungus \u003cem\u003eSporisorium graminicola\u003c/em\u003e in cardiac tissue presents an intriguing finding. \u003cem\u003eS. graminicola\u003c/em\u003e has also been isolated from human fecal samples\u003csup\u003e35\u003c/sup\u003e. Despite the apparent peculiarity of a phytopathogenic fungus inhabiting cardiac tissue, its consistent identification as a significant feature in both the GTEx dataset containing samples from USA and the independent test dataset from Swedish Biobank renders the possibility of mere contamination or chance occurrence highly unlikely. The ability of machine learning models to discriminate distinct microbial signatures across these three tissues in separate datasets suggests that these microorganisms are likely indigenous inhabitants of their respective tissues.\u003c/p\u003e \u003cp\u003eIn the majority of tissues (20 out of 28 tissues), a distinct tissue-specific microbial signature was not detected. However, the absence of such a signature does not necessarily imply tissue sterility. These tissues are more likely to harbor distinct microbial signatures in species that may not have been successfully detected, possibly due to limitations imposed by the poly-A protocol or the lack of reference genomes. Indeed, the methodologies employed by GTEx protocols unequivocally revealed only a fraction of the microbial load present inside these tissues. Another plausible explanation for the absence of tissues' microbial signatures is that the models for these tissues may not have performed optimally, as these tissues may exhibit a more versatile microbiome without a clearly defined tissue-specific signature. Lastly, for the three tissues that did not perform well in living tissues, it is conceivable that these tissues undergo a significant shift in their microbiome after death, or there may be notable differences attributable to geographical variations, as the living subjects were from a different continent than GTEx subjects.\u003c/p\u003e \u003cp\u003eSignificantly, beyond the considerations of the poly-A sequencing protocol, metagenomics and metatranscriptomics necessitates tailored laboratory procedures for efficient DNA or RNA isolation from both gram-positive and gram-negative microorganisms. As these considerations were not initially integrated into the GTEx pipeline, it is evident that the analytical procedure failed to unveil the full spectrum of microorganisms present in the tissues. Moreover, aside from the critical role of laboratory protocols, the bioinformatics analysis protocols also play a pivotal role in discriminating and accurately describing tissue microbiomes. Two metatranscriptomics tools, Kaiju and AGAMEMNON, each employing distinct approaches for read classification, yielded different taxonomic profiles and exhibited variations in performance on the models. To advance microbiome research, it is paramount to develop specific guidelines for both laboratory and \u003cem\u003ein silico\u003c/em\u003e analysis of microbiome data. This step is crucial for ensuring comprehensive and accurate insights into the diversity and composition of microbial communities within tissues.\u003c/p\u003e \u003cp\u003eIt is crucial to consider that, in the majority of the identified microbiomes, the reads classified as belonging to these microbiomes were relatively scarce. This implies that RNA transcripts of these microbiomes were present within the tissues; however, these transcripts did not adequately cover a substantial portion of the genetic material of these microbiomes. This limitation is likely attributable to constraints imposed by the poly-A protocol. The limitations of the poly-A protocol highlight the necessity for additional investigations to validate the existence of microbiomes within tissues, rather than merely isolated microbial RNA transcripts. It is pertinent to acknowledge that bioinformatic analyses in scenarios such as the present study are subject to probabilistic constraints and, on their own, are insufficient to definitively prove the existence of microbiomes. However, they can serve as valuable tools to guide subsequent investigations, indicating tissues and organisms that warrant further scrutiny.\u003c/p\u003e"},{"header":"Conclusion","content":"\u003cp\u003eCollectively, a growing body of evidence substantiates the existence of microbial interactions with human cells in anatomical sites traditionally considered sterile. The identification of encoded signatures specific to certain tissues within normal human tissues suggests for the first time that these tissues harbor a low-biomass microbiome. This exploration signifies a paradigmatic shift, heralding an era in which we acknowledge that our bodies are not solitary entities but rather collaborative ecosystems housing diverse microbial species. We anticipate that our study will serve as a foundational resource, providing crucial guidance for future investigations and facilitating targeted laboratory validations aimed at confirming the presence of microbial communities within internal tissues.\u003c/p\u003e"},{"header":"Methods","content":"\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eGTEx data accession\u003c/h2\u003e \u003cp\u003eThe data used for the analyses described in this manuscript were obtained from dbGaP accession number phs000424.v9.p2 on 04/01/2023. The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. All GTEx RNA-seq data and matched samples protected metadata were accessed via NHGRI Analysis Visualizaion and Informatics Lab-space (AnVIL )(\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://gen3.theanvil.io/login\u003c/span\u003e\u003cspan address=\"https://gen3.theanvil.io/login\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e).\u003c/span\u003e Details of how these data were downloaded are comprehensively described in the AnVIL site (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://anvilproject.org/learn/reference/gtex-v8-free-egress-instructions#step-4---download-object-files\u003c/span\u003e\u003cspan address=\"https://anvilproject.org/learn/reference/gtex-v8-free-egress-instructions#step-4---download-object-files\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). For bulk data acquisition, custom Snakemake files were employed. Due to constraints in storage capacity, an approach was implemented wherein individual sample BAM files were downloaded, and only the unmapped reads (as was originally mapped by GTEx pipeline using STAR (version 2.5.3a)\u003csup\u003e36\u003c/sup\u003e \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://gtexportal.org/home/methods\u003c/span\u003e\u003cspan address=\"https://gtexportal.org/home/methods\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e)\u003c/span\u003e were retained using samtools (version 1.10)\u003csup\u003e37\u003c/sup\u003e.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eGTEx sample and subject QC\u003c/h2\u003e \u003cp\u003eThe original GTEx dataset initially comprised 17,350 RNA-seq samples derived from 948 subjects, all of which had successfully met the GTEx original quality control (QC) criteria. Subsequently, 52 samples were excluded from the analysis due to missing data in critical variables, namely, Genotype or Expression Batch ID, Date of genotype or expression batch, and Total Ischemic time. The dataset was further refined to encompass only 14,478 samples originating from tissues preserved in the PAXgene tissue fixative solution, specifically categorized under \"PAXgene\" and \"Whole Blood:PAXgene\" in the \"current_material_type\" variable. Furthermore, an additional refinement process involved excluding 32 samples originating from tissue sites with fewer than 20 samples, specifically samples from the following tissues : Kidney Medulla, Fallopian Tube, Cervix Extocervix, and Cervix Endocervix were removed. Additionally, 3 samples possessing an RNA Integrity Number (RIN) less than 5 were excluded, along with 139 samples with Total Ischemic time exceeding 1440 minutes (24 hours). Following these rigorous sample quality control measures, a total of 14,304 samples from 942 subjects and spanning 28 different tissues remained in the dataset.\u003c/p\u003e \u003cp\u003eFollowing, quality control of the subjects was performed. 9 subjects were initially excluded due to their ineligibility based on GTEx original criteria. 7 subjects were then removed from the dataset as they had a current cancer diagnosis, and 16 subjects were excluded due to a history of cancer diagnosis within the past 5 years. All the samples from the excluded subjects were removed. Following this comprehensive sample and subject quality control process, the resulting dataset comprised 13,871 RNA-seq samples originating from 28 distinct tissues (refer to S\u003cb\u003eupplementary file 25\u003c/b\u003e) and derived from a cohort of 910 subjects. The quality control procedures were implemented utilizing custom Python scripts, available on the corresponding GitHub repository.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eIsolation of unmapped reads and quality control\u003c/h2\u003e \u003cp\u003eFor each sample, sequencing reads that failed to align with the human reference genome, as indicated by the mapping information in the raw BAM files obtained from GTEx, were selectively retained. To isolate the unmapped reads where both paired reads were unaligned, and to eliminate reads classified as non-primary alignments, bioinformatic tool Samtools was employed using the arguments \u003cem\u003e\u0026ldquo;-f 12 -F 256\u003c/em\u003e\u0026ldquo;.\u003c/p\u003e \u003cp\u003eThe unmapped reads underwent a comprehensive quality control process, involving the exclusion of reads with a length shorter than 35 nucleotides (nt), given that the GTEx normal sequencing length was 76 nt. Additionally, steps were taken to remove adapters and perform quality trimming with a threshold at Phred quality score of 15. These quality control procedures were executed using Atropos (version 1.1.31)\u003csup\u003e38\u003c/sup\u003e. To streamline and automate this sequence of steps, custom Snakemake scripts were implemented, available on the corresponding GitHub repository.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eTaxonomic Assignment using AGAMEMNON\u003c/h2\u003e \u003cp\u003eTo construct the reference database utilized for the taxonomic algorithm AGAMEMNON (version 0.1.0), a custom Bash script was employed on 02/14/2023 to download microbial genomes from RefSeq.\u0026nbsp;The dataset comprised all bacterial representative or reference genomes with complete genome assembly level (4,034 bacterial genomes), all archaeal genomes with complete genome assembly level (489 archaea genomes), and all viral genomes with complete genome assembly level (11,259 viral genomes). Furthermore, fungal transcriptomes from representative or reference genomes with complete genome or chromosome assembly levels (81 fungal transcriptomes) were also included in the reference database. In the case of fungi, the transcriptome data was specifically employed to account for the intricacies of the splicing process. This approach ensures a more comprehensive representation of fungal genomic information, taking into consideration the variations introduced during the splicing of transcripts.\u003c/p\u003e \u003cp\u003eThe quality-controlled non-human reads of each sample underwent mapping against the custom reference database described above, using the Puffaligner algorithm\u003csup\u003e39\u003c/sup\u003e. For Puffaligner, the flag \u0026ldquo;--noOrphans\u0026rdquo; was used in order to discard the orphans reads. The term \u0026ldquo;orphan\u0026rdquo; refers to one end of paired-end read that is confidently aligned to some genomic position, but for which the other read end is not jointly aligned nearby (and paired). Puffaligner aligns the reads to the compiled microbial genomes, enabling the identification of reads of microbial origin within the samples. Puffaligner is a fast, sensitive and accurate aligner based on a compacted sequence graph and is meticulously crafted to embody a dual emphasis on high sensitivity in alignment tasks and efficient computational performance. Its design capitalizes on the utilization of a colored compacted de Bruijn graph to efficiently identify and factor out recurring sub-sequences within the reference.\u003c/p\u003e \u003cp\u003eAs highlighted by Gitawi \u003cem\u003eet al.\u003c/em\u003e\u003csup\u003e40\u003c/sup\u003e, in order to be sure that no human read had succeeded to infiltrate in our analysis, we remapped the reads that were classified by Puffaligner as reads of microbial origin to the most recent human genome assembly (T2T-CHM13v2.0 from T2T Consortium) using bowtie2 (version 2.2.3)\u003csup\u003e41, 42\u003c/sup\u003e with the \"preset\" parameters of \u0026ldquo;--very-sensitive\u0026rdquo; and kept all the reads that were not mapped concordantly. Custom snakemake scrips were used to automate the process. Then by using only the remaining reads, we analysed them with AGAMEMNON (a changed version, deposited in the github repo of the present study) to find the microbial profile of each sample. AGAMEMNON represents a metagenomics and metatranscriptomics algorithm, notable for its integration of a time and space-efficient indexing scheme. This feature facilitates rapid pattern matching, allowing for the efficient indexing and analysis of extensive datasets using commonly available computational resources. In the abundance estimation step, the primary approach relies on the expectation maximization (EM) algorithm. The goal is to maximize the likelihood of observed reads by iteratively adjusting the abundance values linked to various taxa. At last, after running AGAMEMNON, the taxonomic profile of each sample was generated.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eDiversity metrics and Core microbiome per tissue\u003c/h2\u003e \u003cp\u003eThe taxonomic profiles generated by AGAMEMNON at the species level were employed to calculate the microbiome richness of each sample. Microbiome richness was defined as the count of species with non-zero abundance in each sample. To compute the Shannon diversity index, the 'diversity' function from the R package vegan (version '2.6.4') was utilized.\u003c/p\u003e \u003cp\u003eThe core microbiome for each tissue was defined by considering all species present in at least 10% of the respective tissue\u0026rsquo;s samples. The 10% threshold was chosen to exclude species identified in a limited number of samples, that may be indicative of opportunistic infections, and to retain only those species consistently present in tissue samples. Given the inherent diversity and dynamic nature of the microbiome, a relatively low threshold (10%) was applied to avoid excluding too many species. For comparisons of the core microbiome across tissues, UpSet plots were generated using the R library UpSet (version 1.4.0).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section2\"\u003e \u003ch2\u003eNormalization of taxonomic profiles\u003c/h2\u003e \u003cp\u003eTo address biases, particularly due to differences in sequencing depth among samples, we employed Cumulative Sum Scaling (CSS) normalization in the taxonomic profiles. CSS, functioning as a median-like quantile normalization method, corrects for variations in sampling depth or library size. Unlike standard relative abundance normalization, which rescales all samples to a uniform total sum (\u003cem\u003ee.g.\u003c/em\u003e, 100%), CSS retains variability in total counts across samples. This normalization method adjusts samples based on a subset (quartile) of lower-abundance taxa that remain relatively constant and independent, mitigating the impact of high-abundance taxa that may dominate a study. For CSS normalization, the metagenomeSeq library (version 1.40.0) in R was employed. This normalization approach was also applied to normalize taxonomic profiles at the genus level and functional profiles generated by HumanN3.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003eML models 1vsAllOther27tissues\u003c/h2\u003e \u003cp\u003eTo discern potential biological significance within the taxonomic profiles of samples, we endeavored to construct ML models for each tissue to discriminate it from all the other 27 tissues. Utilizing Gradient Boosting ML models known for their efficacy in classification tasks and resilience to imbalanced datasets, we trained stochastic gradient boosting machine (GBM) learning models. The models were implemented and fine-tuned using the R libraries GBM (version 2.1.8.1), Caret (version 6.0.94), and PRROC (version 1.3.1) (for calculating AUROC and AUPR values). The training and testing phases occurred on distinct, randomly selected, stratified sampling splits of 70% and 30% of the data, respectively, with a fixed random seed for reproducibility. Data CSS normalization was performed separately in each split. After normalisation of each split, only the total core species of the 28 tissues were retained for training the model as we were interested in the signatures present in the consistent microbiome within tissues. A two-fold cross-validation strategy was employed to create subsets of the training data and we conducted a basic grid search optimization of GBM parameters, with interaction depth (1, 2, or 3) and the number of trees (50, 100, or 150), with the aim of maximizing the area under the receiver operating characteristic (AUROC) of the final, fully trained model. The learning rate (shrinkage) remained constant at 0.1, and the minimum observations per node was fixed at 3. Final model performance metrics, including AUROC and area under the precision-recall curve (AUPR) values, were derived by applying the final model to the unseen 30% holdout test set. For calculating AUROC and AUPR the library PRROC was used. Given the stochastic nature of GBM models, 100 iterations of GBM models were created for each tissue model to calculate the mean AUROC and mean AUPR values, along with 95% confidence intervals. After that as the minority class for the different tissues were quite different across samples, we calculated a \"normalized performance\u0026rdquo;, the relative AUPR, which was given by the following formula: Relative AUPR is equal to the ratio of AUPR(model)​ to AUPR(random). To expedite processes, the script for ML model creation utilized 10 cores, leveraging the R libraries \u0026ldquo;parallel\u0026rdquo; (version 4.2.2) and \u0026ldquo;doMC\u0026rdquo; (version 1.3.8).\u003c/p\u003e \u003cp\u003eFor the correlation of mean AUROC and relative mean AUPR of the different models with the size of the samples of each tissue, Spearman's rank correlation coefficient was applied using the \u003cem\u003ecor.test\u003c/em\u003e from the R library stats (version 4.2.2).\u003c/p\u003e \u003cp\u003eThe same approach was used for creating the GBM models for the genera level. As features of the models were used the total core genera microbiome of the 28 tissues, which consisted of 738 genera.\u003c/p\u003e \u003cp\u003e \u003cb\u003eIn silico\u003c/b\u003e \u003cb\u003eContamination\u003c/b\u003e\u003c/p\u003e \u003cp\u003eFor the \u003cem\u003ein silico\u003c/em\u003e contamination approach, we strategically introduced 12 distinct \u003cem\u003ein silico\u003c/em\u003e contaminants into our dataset, simulating various potential types of contamination that may occur within the GTEx framework. The contaminants were classified into two primary categories: high-volume contaminants, characterized by their high volume (high number of reads attributed to them) on a small subset of samples, positing that their detection would be possible if they affected a larger sample pool; and low-volume contaminants, which, although affecting a greater number of samples, did so at lower volumes (small number of reads attributed to them), thereby evading their detection. Consequently, we incorporated 6 high-volume and 6 low-volume contaminants into our dataset. The diverse types of contaminants is depicted in the Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e. We conducted 100 iterations, randomly selecting contaminated samples in each iteration using a custom Python script. The \u003cem\u003ein silico\u003c/em\u003e contaminations were introduced post CSS normalization. Subsequently, Gradient Boosting Machine (GBM) models were constructed for the 11 tissues exhibiting a tissue-specific microbiome. These models aimed to distinguish each tissue from the other 10 tissues using similar methodology as described previously. For these models, the entire core microbiome of the aforementioned 11 tissues, comprising 1612 distinct microbial species, was retained for analysis. Additionally, 100 iterations of GBM models without any type of \u003cem\u003ein silico\u003c/em\u003e contamination were generated to discriminate one tissue from the other 10 tissues, enabling comparison of results. Features\u0026rsquo; importance scores of non-zero model features were estimated using the GBM and Caret packages. The percentage of contribution of a specific feature to the model's predictions was computed by dividing its importance score by the sum of all features\u0026rsquo; importance scores for the given model. To assess the performance of contaminated versus uncontaminated models, we employed Wilcoxon's test on the AUROC and AURP values derived from the 100 iterations from contaminated and non-contaminated models across the 11 tissues, utilizing the R library stats (version 4.2.2).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003eML models 1vs7 tissues\u003c/h2\u003e \u003cp\u003eTo focus exclusively on the signatures identified within the 8 contamination-resilient tissues (Blood, Testis, Colon tissue, Stomach tissue, Muscle, Bladder, Liver, and Heart) and analyze their distinctions, we reconstructed GBM models utilizing only samples derived from these tissues. These models utilized the concatenation of the core microbiome of these 8 tissues, consisting of 1597 species for the species models. The configurations of these models are similar to the ones described in the section \u003cspan refid=\"Sec17\" class=\"InternalRef\"\u003eML models 1vsAllOther27tissues\u003c/span\u003e. To compute the mean AUROC and AUPR, along with their associated margins of error, 100 iterations of each model were executed. Features\u0026rsquo; importance scores for the model features were estimated as previously described, employing the first iteration of the models. The identified signatures from these models underscored the differences in the microbial compositions inside these 8 tissues.\u003c/p\u003e \u003cp\u003eThe same methodology was applied to construct GBM models at the genus level. The features for these models comprised the total core microbiome (at genus level) present within the 8 tissues, encompassing 695 genera.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003eML models with Kaiju\u003c/h2\u003e \u003cp\u003eTo validate the robustness of findings for the 8 tissue-specific signatures, we employed a different taxonomic tool to identify the taxonomic profiles of samples from the 8 contamination-resilient tissues. Utilizing the QCed reads (post-Atropos processed reads), we employed Kaiju (version 1.9.2) for taxonomic identification with 10 parallel threads allocated to each sample. Kaiju is a rapid and sensitive program designed for taxonomic classification of high-throughput sequencing reads. Each sequencing read is assigned to a taxon within the NCBI taxonomy by comparing it against a reference protein database. We utilized the microbial subset of the NCBI BLAST non-redundant protein database (nr), incorporating fungi and microbial eukaryotes as reference database. The database was downloaded on 03/30/2023 from \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://bioinformatics-centre.github.io/kaiju/downloads.html\u003c/span\u003e\u003cspan address=\"https://bioinformatics-centre.github.io/kaiju/downloads.html\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e. Reads were translated into amino acid sequences, which were then queried against the database using a modified backward search on a memory-efficient implementation of the Burrows-Wheeler transform, identifying maximum exact matches (MEMs) and optionally allowing mismatches. The tool \u003cem\u003ekaiju2table\u003c/em\u003e was employed to convert Kaiju's output files into a summary table for the taxonomic rank of \"species.\"\u003c/p\u003e \u003cp\u003eUsing these taxonomic profiles, GBM models were reconstructed for each tissue to discriminate it from all other 7 tissues, following a similar methodology as described previously. In these GBM models, the features used comprised the concatenation of the core microbiome of the 8 tissues based on the Kaiju tool, encompassing 1864 microbial species. To determine the mean AUROC and AUPR, along with their respective margins of error, 100 iterations of each model were executed. Feature importance scores for the model features were estimated in one iteration of the models, as described earlier.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003eFunctional Assignment using HumanN3\u003c/h2\u003e \u003cp\u003eTo characterize the functional profile of samples from the GTEx consortium, reads that were identified as microbial by AGAMEMNON were analyzed using HUMAnN3 (version 3.6.1) for profiling the abundance of microbial genes. HUMAnN3 is a method designed to efficiently and accurately conduct functional profiling in metagenomic or metatranscriptomic sequencing data. Due to constraints in computational resources and time, HUMAnN3 was executed with the option \"--bypass-nucleotide-search\" to skip all alignment steps before the translated search. For this translated search, the full UniRef90 database (version 201901b) served as the reference.\u003c/p\u003e \u003cp\u003eSubsequently, the functional profiles of all samples were concatenated and normalized using CSS normalization. GBM models were then created following the approach described previously, utilizing the 1-vs-7 tissues strategy and incorporating only the core gene repertoire. This repertoire consisted of genes present in at least 10% of each tissue\u0026rsquo;s samples, amounting to 1,831 different microbial genes across the 8 tissues.\u003c/p\u003e \u003cp\u003eFor models integrating both functional and taxonomic profiles, the core gene repertoire (1,831 genes) and the core species microbiome (1,597 species) were utilized in the creation of GBM models. Each component was separately normalized using CSS normalization, and then the normalized gene and species taxonomy profiles were merged for the training and testing datasets. The GBM models were constructed using the same parameters as described earlier. To calculate the mean AUROC and AUPR, along with their margins of error, 100 iterations of each model were generated.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec21\" class=\"Section2\"\u003e \u003ch2\u003eFactors associated with tissue microbiome\u003c/h2\u003e \u003cp\u003eTo investigate the influence of various traits (age, Body Mass Index - BMI, smoking status, drinking status, ancestry), as well as disease history (hypertension history, ischemic heart history, diabetes II history, diabetes I history, arthritis history, seizures history, schizophrenia history, rheumatoid arthritis history, liver disease history, dialysis treatment, depression history, COPD or CLRD history, cerebrovascular disease history, asthma history, alzheimer or dementia history), on the core microbiome of the 8 tissues with contamination-resilient tissue-specific microbiomes, GBM models were created for each trait and tissue combination with aim to discriminate in each tissue the samples with the specific trait from the samples without this trait. Tissues with fewer than 20 samples within each different group of specific trait/disease were excluded due to insufficient data for model creation.\u003c/p\u003e \u003cp\u003eFor continuous traits such as age and BMI, GBM models were constructed using 5-fold cross-validation and Root Mean Squared Error (RMSE) was used to select the optimal model with the smallest RMSE value. Mean Absolute Error (MAE) and R-squared were calculated using the \"postResample\" function from the Caret package. The mean RMSE, MAE, and R-squared were computed across 100 iterations of each model.\u003c/p\u003e \u003cp\u003eFor categorical traits (smoking, drinking, sex, and ancestry), as well as models for disease history, GBM models were developed similarly to tissue models but with 4-fold cross-validation to mitigate overfitting. Mean AUROC and mean AUPR were computed across 100 iterations of each model. For the ancestry trait, samples with ancestries other than \"White\" or \"Black or African American\" were removed due to limited representation.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec22\" class=\"Section2\"\u003e \u003ch2\u003eValidating significant ML models with living samples\u003c/h2\u003e \u003cp\u003eTo validate the previously created 1-vs-7 tissues models in living tissues, we analyzed data of the 8 aforementioned tissues from project PRJEB4337, comprising RNA-seq data from living samples of the Swedish Biobank. From this project, only the samples from the 8 tissues were kept for testing the models, unfortunately this project did not contain samples from muscle tissue and blood. The RNA-seq data from this project were generated using a poly-A selection protocol and processed similarly to GTEx samples. Firstly, the fastq files were aligned to the same human genome as GTEx using STAR (version 2.7.10b). Then, the unmapped reads were isolated and quality controlled using Atropos with the same configurations as described in the GTEx analysis. Subsequently, the QCed reads were mapped to the microbial database using Puffaligner, and the microbial reads were remapped to the latest human genome using Bowtie2. Lastly, the taxonomical microbial profile of the 38 samples was assessed using Agamemnon. The identical pipeline was employed for both the GTEx and PRJEB4337 samples to ensure consistency.\u003c/p\u003e \u003cp\u003eNext, the previously created 1-vs-7 tissues models (100 iterations per tissue) from the GTEx dataset were tested using the living samples from PRJEB4337. The AUROC and AUPR were calculated using the PRROC library (version 1.3.1). To verify that the observed high performance in the living tissues was attributed to the specific microbial signatures within these tissues, 100 iterations of GBM models were created using the GTEx data as described in the 1-vs-7 tissues models section, but the tissue labels of the samples in the training dataset were randomly shuffled to create random models without any specific tissue information. These random models were then tested with the living dataset to assess their performance. The performances in the living dataset of the 100 iterations of the original 1-vs-7 tissues models were compared with the performances of the 100 iterations of the random 1-vs-7 tissues models using a Wilcoxon test. This comparison was crucial for determining whether the performance observed in the living tissues was attributable to the specific microbial signatures within those tissues or was merely a result of chance.\u003c/p\u003e \u003cp\u003eLastly, samples from the PRJEB4337 project were also analyzed using the Kaiju tool to create their microbial profiles, following the same pipeline employed for GTEx samples. Subsequently, the 100 iterations of GBM models, created using Kaiju taxonomic profiles of GTEx data as described in the 'ML Models with Kaiju' section, were tested on the samples from the PRJEB4337 project.\u003c/p\u003e \u003cdiv id=\"Sec23\" class=\"Section3\"\u003e \u003ch2\u003eSystem configurations\u003c/h2\u003e \u003cp\u003eThe analyses described in this study were conducted on a computational platform equipped with an Intel\u0026reg; Xeon(R) Gold 6226R CPU @ 2.90GHz and 376 GB of RAM. The operating system used was Ubuntu Linux version 20.04. Data processing and modeling were implemented using Python version 3.9.15 and R version 4.2.2.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e"},{"header":"Declarations","content":"\u003ch2\u003eAvailability of data and materials\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eThe taxonomic profiles of GTEx samples, generated using AGAMEMNON and Kaiju tools, along with their functional profiles generated via Humann3, are accessible in Zenodo (link \u003cu\u003ehttps://zenodo.org/uploads/10980664\u003c/u\u003e). Additionally, taxonomic and functional profiles for the distinctive living dataset (NCBI ID PRJEB4337) are also available on the same link, along with its metadata R object used in the analysis. Although the full metadata for GTEx samples is restricted, a subset is available at https://www.gtexportal.org/home/downloads/adult-gtex/metadata; access to the complete metadata requires approval through a dbGaP application at https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000424.v8.p2.\u003c/p\u003e\n\u003cp\u003eAll programming scripts used to download and analyse the data of the GTEx and NCBI project PRJEB4337 as well as the pipelines for \u0026nbsp;normalization, \u003cem\u003ein silico\u0026nbsp;\u003c/em\u003econtamination, ML models can be found at our GitHub repository \u003cu\u003ehttps://github.com/dianalabgr/GTEx_microbiome_analysis\u003c/u\u003e. These scripts are designed to directly process the summarized count data available on Zenodo. However, analysis of microbial signatures concerning phenotypic traits and medical history necessitates access to the complete metadata. For the other analytical pursuits, the open access metadata suffices; nonetheless, it is essential to accommodate variations in column names between the open access and protected access metadata files.\u003c/p\u003e\n\u003ch2\u003eAcknowledgements\u003c/strong\u003e\u003c/h2\u003e\n\u003cp\u003eThe authors acknowledge the members of DIANA-lab for the very useful comments and ideas, and specifically Spyros Tatsoglou, PhD for his help with the procedures of accessing the data from GTEx Consortium.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp;\u003cstrong\u003eAuthors\u0026rsquo; contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAnargyros Skoulakis: Conceptualization, Methodology, Software, Visualization, Writing \u0026ndash; original draft, Writing \u0026ndash; review \u0026amp; editing. Giorgos Skoufos: Conceptualization, Methodology, Writing \u0026ndash; review \u0026amp; editing. Armen Ovsepian: Methodology, Writing \u0026ndash; review \u0026amp; editing. Artemis G. Hatzigeorgiou: Conceptualization, Supervision, Funding acquisition, Writing \u0026ndash; review \u0026amp; editing\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThis research has been co‐financed by the European Regional Development Fund of the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship and Innovation, and Greece 2.0, under the call \u0026ldquo;RESEARCH \u0026ndash; CREATE \u0026ndash; INNOVATE\u0026rdquo; (ID 16971), with project id: TAEDK-06179. Also, specifically for Anargyros Skoulakis, the research work was supported also by the Hellenic Foundation for Research and Innovation (HFRI) under the 5th Call for HFRI PhD Fellowships (Fellowship Number: 20480)\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent for publication\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConflict of Interest\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eMerino, N., Aronson, H. S., Bojanova, D. P., Feyhl-Buska, J., Wong, M. L., Zhang, S., \u0026amp; Giovannelli, D. (2019). Living at the extremes: extremophiles and the limits of life in a planetary context. Frontiers in microbiology, 10, 780.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eUrsell, L. K., Metcalf, J. L., Parfrey, L. W., \u0026amp; Knight, R. (2012). Defining the human microbiome. Nutrition reviews, 70(suppl_1), S38-S44.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eThe Human Microbiome Project Consortium. A framework for human microbiome research. Nature 486, 215\u0026ndash;221 (2012). \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1038/nature11209\u003c/span\u003e\u003cspan address=\"10.1038/nature11209\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eThe Human Microbiome Project Consortium. Structure, function and diversity of the healthy human microbiome. \u003cem\u003eNature\u003c/em\u003e 486, 207\u0026ndash;214 (2012). \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1038/nature11234\u003c/span\u003e\u003cspan address=\"10.1038/nature11234\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eThe Integrative HMP (iHMP) Research Network Consortium. The Integrative Human Microbiome Project. \u003cem\u003eNature\u003c/em\u003e 569, 641\u0026ndash;648 (2019). \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1038/s41586-019-1238-8\u003c/span\u003e\u003cspan address=\"10.1038/s41586-019-1238-8\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNatalini, J.G., Singh, S. \u0026amp; Segal, L.N. The dynamic lung microbiome in health and disease. \u003cem\u003eNat Rev Microbiol\u003c/em\u003e 21, 222\u0026ndash;235 (2023). \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.1038/s41579-022-00821-x\u003c/span\u003e\u003cspan address=\"10.1038/s41579-022-00821-x\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChen, X., Lu, Y., Chen, T., \u0026amp; Li, R. (2021). The female vaginal microbiome in health and bacterial vaginosis. Frontiers in cellular and infection microbiology, 11, 631972.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePoore, G. D., Kopylova, E., Zhu, Q., Carpenter, C., Fraraccio, S., Wandro, S., \u0026hellip; Knight, R. (2020). Microbiome analyses of blood and tissues suggest cancer diagnostic approach. Nature, 579(7800), 567\u0026ndash;574.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDohlman, A. B., Mendoza, D. A., Ding, S., Gao, M., Dressman, H., Iliev, I. D., \u0026hellip; Shen, X. (2021). The cancer microbiome atlas: a pan-cancer comparative analysis to distinguish tissue-resident microbiota from contaminants. Cell host \u0026amp; microbe, 29(2), 281\u0026ndash;298.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNarunsky-Haziza, L., Sepich-Poore, G. D., Livyatan, I., Asraf, O., Martino, C., Nejman, D., \u0026hellip; Straussman, R. (2022). Pan-cancer analyses reveal cancer-type-specific fungal ecologies and bacteriome interactions. Cell, 185(20), 3789\u0026ndash;3806.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNejman, D., Livyatan, I., Fuks, G., Gavert, N., Zwang, Y., Geller, L. T., \u0026hellip; Straussman, R. (2020). The human tumor microbiome is composed of tumor type\u0026ndash;specific intracellular bacteria. Science, 368(6494), 973\u0026ndash;980.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRiquelme, E., Zhang, Y., Zhang, L., Montiel, M., Zoltan, M., Dong, W., \u0026hellip; McAllister, F. (2019). Tumor microbiome diversity and composition influence pancreatic cancer outcomes. Cell, 178(4), 795\u0026ndash;806.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGhaddar, B., Biswas, A., Harris, C., Omary, M. B., Carpizo, D. R., Blaser, M. J., \u0026amp; De, S. (2022). Tumor microbiome links cellular programs and immunity in pancreatic cancer. Cancer Cell, 40(10), 1240\u0026ndash;1253.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAykut, B., Pushalkar, S., Chen, R., Li, Q., Abengozar, R., Kim, J. I., \u0026hellip; Miller, G. (2019). The fungal mycobiome promotes pancreatic oncogenesis via activation of MBL. Nature, 574(7777), 264\u0026ndash;267.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMahmoudabadi, G., Consortium, T. S., \u0026amp; Quake, S. R. (2022). Single cell transcriptomics reveals the hidden microbiomes of human tissues. bioRxiv, 2022\u0026ndash;10.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHieken, T. J., Chen, J., Hoskin, T. L., Walther-Antonio, M., Johnson, S., Ramaker, S., \u0026hellip; Degnim, A. C. (2016). The microbiome of aseptically collected human breast tissue in benign and malignant disease. Scientific reports, 6(1), 30751.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLluch, J., Servant, F., Pa\u0026iuml;ss\u0026eacute;, S., Valle, C., Valiere, S., Kuchly, C., \u0026hellip; Lelouvier, B. (2015). The characterization of novel tissue microbiota using an optimized 16S metagenomic sequencing pipeline. PloS one, 10(11), e0142334.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLonsdale, J., Thomas, J., Salvatore, M., Phillips, R., Lo, E., Shad, S., \u0026hellip; Moore, H. F. (2013). The genotype-tissue expression (GTEx) project. Nature genetics, 45(6), 580\u0026ndash;585.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSkoufos, G., Almodaresi, F., Zakeri, M., Paulson, J. N., Patro, R., Hatzigeorgiou, A. G., \u0026amp; Vlachos, I. S. (2022). AGAMEMNON: an Accurate metaGenomics And MEtatranscriptoMics quaNtificatiON analysis suite. Genome biology, 23(1), 39.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMenzel, P., Ng, K. L., \u0026amp; Krogh, A. (2016). Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nature communications, 7(1), 11257.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePaulson, J. N., Stine, O. C., Bravo, H. C., \u0026amp; Pop, M. (2013). Differential abundance analysis for microbial marker-gene surveys. Nature methods, 10(12), 1200\u0026ndash;1202.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFriedman, J. H. (2001). Greedy function approximation: a gradient boosting machine. Annals of statistics, 1189\u0026ndash;1232.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAlmodaresi, F., Sarkar, H., Srivastava, A., \u0026amp; Patro, R. (2018). A space and time-efficient index for the compacted colored de Bruijn graph. \u003cem\u003eBioinformatics\u003c/em\u003e, \u003cem\u003e34\u003c/em\u003e(13), i169-i177.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCarlier, J. P., Bedora-Faure, M., K'ouas, G., Alauzet, C., \u0026amp; Mory, F. (2010). Proposal to unify Clostridium orbiscindens Winter et al. 1991 and Eubacterium plautii (S\u0026eacute;guin 1928) Hofstad and Aasjord 1982, with description of Flavonifractor plautii gen. nov., comb. nov., and reassignment of Bacteroides capillosus to Pseudoflavonifractor capillosus gen. nov., comb. nov. International journal of systematic and evolutionary microbiology, 60(3), 585\u0026ndash;590.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMorita, T., Koike, H., Hagiwara, H., Ito, E., Machida, M., Sato, S., \u0026hellip; Kitamoto, D.(2014). Genome and transcriptome analysis of the basidiomycetous yeast Pseudozyma antarctica producing extracellular glycolipids, mannosylerythritol lipids. \u003cem\u003ePLoS One\u003c/em\u003e, \u003cem\u003e9\u003c/em\u003e(2), e86490.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMorita, T., Fukuoka, T., Imura, T., \u0026amp; Kitamoto, D. (2015). Mannosylerythritol lipids: production and applications. Journal of Oleo Science, 64(2), 133\u0026ndash;141.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGonelli, A., Boccia, S., Boni, M., Pozzoli, A., Rizzo, C., Querzoli, P., \u0026hellip; Di Luca,D. (2001). Human herpesvirus 7 is latent in gastric mucosa. Journal of medical virology,63(4), 277\u0026ndash;283.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGarcia-Montojo, M., Doucet-O\u0026rsquo;Hare, T., Henderson, L., \u0026amp; Nath, A. (2018). Human endogenous retrovirus-K (HML-2): a comprehensive review. Critical reviews in microbiology, 44(6), 715\u0026ndash;738.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBroderick, N. A., \u0026amp; Nagy, L. (2022). Bacteria may be in the liver, but the jury is still out. The Journal of Clinical Investigation, 132(8).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFagerberg, L., Hallstr\u0026ouml;m, B. M., Oksvold, P., Kampf, C., Djureinovic, D., Odeberg,J., \u0026hellip; Uhl\u0026eacute;n, M. (2014). Analysis of the human tissue-specific expression by genome-wide integration of transcriptomics and antibody-based proteomics. Molecular \u0026amp; cellular proteomics, 13(2), 397\u0026ndash;406.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBeghini, F., McIver, L. J., Blanco-M\u0026iacute;guez, A., Dubois, L., Asnicar, F., Maharjan,S., \u0026hellip; Segata, N. (2021). Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3. elife, 10, e65088.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCiernikova, S., Sevcikova, A., Stevurkova, V., \u0026amp; Mego, M. (2022). Tumor microbiome\u0026ndash;an integral part of the tumor microenvironment. Frontiers in Oncology, 12, 1063100.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRajilić-Stojanović, M., \u0026amp; De Vos, W. M. (2014). The first 1000 cultured species of the human gastrointestinal microbiota. FEMS microbiology reviews, 38(5), 996\u0026ndash;1047.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVanderHeyden, T. R., Yong, S. L., Breitschwerdt, E. B., Maggi, R. G., Mihalik, A. R., Parada, J. P., \u0026amp; Fimmel, C. J. (2012). Granulomatous hepatitis due to Bartonella henselae infection in an immunocompetent patient. BMC Infectious Diseases, 12, 1\u0026ndash;7.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNatalia, S., Luiza, H., Anna, S. C., Katarzyna, T., Schmidt, M., Łukasz, P., \u0026hellip; Philips,A. (2023). Host Factors Associated with Gut Mycobiome Structure. MSystems, 8(2).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., \u0026hellip; Gingeras,T. R. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29(1), 15\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDanecek, P., Bonfield, J. K., Liddle, J., Marshall, J., Ohan, V., Pollard, M. O.,\u0026hellip; Li, H. (2021). Twelve years of SAMtools and BCFtools. Gigascience, 10(2), giab008.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDidion, J. P., Martin, M., \u0026amp; Collins, F. S. (2017). Atropos: specific, sensitive, and speedy trimming of sequencing reads. PeerJ, 5, e3720.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAlmodaresi, F., Zakeri, M., \u0026amp; Patro, R. (2021). PuffAligner: a fast, efficient and accurate aligner based on the Pufferfish index. Bioinformatics, 37(22), 4048\u0026ndash;4055.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGihawi, A., Ge, Y., Lu, J., Puiu, D., Xu, A., Cooper, C. S., \u0026hellip; Salzberg, S. L. (2023).Major data analysis errors invalidate cancer microbiome findings. Mbio, 14(5), e01607-23.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLangmead, B., Wilks, C., Antonescu, V., \u0026amp; Charles, R. (2019). Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics, 35(3), 421\u0026ndash;432.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLangmead, B., \u0026amp; Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nature methods, 9(4), 357\u0026ndash;359.\u003c/span\u003e\u003c/li\u003e \u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"human tissues, human organs, tissues, microbiome, microbial communities, microorganisms, microbial signatures, machine learning","lastPublishedDoi":"10.21203/rs.3.rs-4610181/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4610181/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003ch2\u003eBackground\u003c/h2\u003e \u003cp\u003eThe presence of microbiomes within healthy human internal organs still remains under question. Our study endeavors to discern distinctive microbial signatures within normal human internal tissues from the Genotype-Tissue Expression (GTEx) consortium through the application of machine learning (ML) models. Machine learning models were developed to classify each tissue type against all others using solely the microbial profile of the tissues. The identification of tissue-specific microbial signatures that enable machine learning models to accurately predict different tissue types would suggest the presence of distinct microbial communities within tissues. To identify these microbial signatures, we analyzed 13,871 normal RNA-seq samples from 28 tissues obtained from the GTEx database.\u003c/p\u003e\u003ch2\u003eResults\u003c/h2\u003e \u003cp\u003eTissue-specific microbial signatures were identified in 11 out of the 28 tissues while the signatures for 8 tissues (Muscle, Heart, Stomach, Colon tissue, Testis, Blood, Liver, and Bladder tissue) demonstrated resilience to \u003cem\u003ein silico\u003c/em\u003e contamination. As the GTEx samples were derived from post-mortem biopsies, healthy living tissue samples were also analyzed to further validate the results. The models for Heart, Colon tissue, and Liver displayed high discriminatory performance in the living dataset, suggesting the presence of a tissue-specific microbiome for these tissues even in a living state. Notably, the most crucial features were the fungus \u003cem\u003eSporisorium graminicola\u003c/em\u003e for the heart, the gram-positive bacterium \u003cem\u003eFlavonifractor plautii\u003c/em\u003e for the colon tissue, and the gram-negative bacterium \u003cem\u003eBartonella machadoae\u003c/em\u003e for the liver. The presence of contamination-resilient tissue-specific microbial signatures implies that microbial colonization is not a random occurrence.\u003c/p\u003e\u003ch2\u003eConclusions\u003c/h2\u003e \u003cp\u003eThe presence of tissue-specific microbial signatures in certain tissues suggests that these organs are not devoid of microorganisms even in healthy conditions and probably they harbor low-biomass microbial communities unique to each tissue. The discoveries presented here confront the enduring dogma positing the sterility of internal tissues, yet further validation through controlled laboratory experiments is imperative to substantiate this hypothesis. Exploring the microbiome of internal tissues holds promise for elucidating the pathophysiology underlying both health and a spectrum of diseases, including sepsis, inflammation, and cancer.\u003c/p\u003e","manuscriptTitle":"Machine learning models reveal microbial signatures in healthy human tissues, challenging the sterility of human organs","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-07-16 09:17:07","doi":"10.21203/rs.3.rs-4610181/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"8ea6f185-b92a-4cb1-b83d-21ca7c3b4674","owner":[],"postedDate":"July 16th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2025-02-04T16:11:29+00:00","versionOfRecord":{"articleIdentity":"rs-4610181","link":"https://doi.org/10.3389/fmicb.2024.1512304","journal":{"identity":"frontiers-in-microbiology","isVorOnly":true,"title":"Frontiers in Microbiology"},"publishedOn":"2025-01-26 00:00:00","publishedOnDateReadable":"January 26th, 2025"},"versionCreatedAt":"2024-07-16 09:17:07","video":"","vorDoi":"10.3389/fmicb.2024.1512304","vorDoiUrl":"https://doi.org/10.3389/fmicb.2024.1512304","workflowStages":[]},"version":"v1","identity":"rs-4610181","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4610181","identity":"rs-4610181","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}