Why we missed it? Computational analysis reveals distribution patterns of Malassezia furfur, the etiological agent of Pityriasis versicolor, in&nbsp;skin metagenomes

doi:10.21203/rs.3.rs-122326/v1

Download PDF

Research Article

Why we missed it? Computational analysis reveals distribution patterns of Malassezia furfur, the etiological agent of Pityriasis versicolor, in skin metagenomes

https://doi.org/10.21203/rs.3.rs-122326/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Malassezia furfur is the main causative species in pityriasis versicolor infections and is now widely believed to be a part of the skin microbiota, yet it was systematically missed in the early skin metagenomic studies. Here, we curated a specific set of M. furfur sequences and used them to reanalyze publicly available skin metagenomes to computationally investigate the distribution of M. furfur and its relative abundance at different skin sites. To this end, we used BLASTN to match and align these marker genes to the selected metagenomic datasets and estimated M. furfur relative abundance as the number of BLASTN hits per million metagenomic reads. We found a relative enrichment of M. furfur in the retro auricular crease, the antecubital fossa, and the forehead. Among skin categories, sebaceous areas were the most significantly enriched in M. furfur, while in terms of exposure/occlusion, exposed areas had the highest abundance. This work will facilitate and allow the estimation and correction of past estimates of this important fungal species in shot gun metagenomic data sets.

General Microbiology

Infectious Diseases

metagenomics

bioinformatics

high-throughput sequencing data analysis

mycobiomics

microbiome

With a surface area of 1.8 m², the skin is the largest organ in the human body and the first barrier against invasion by pathogens [1].

The skin microbiota is the assemblage of commensal bacteria and fungi that reside on the skin. Like other components of the human microbiome, which are in a dynamic state of change—akin to a cloud of uncertainty [2], the skin microbiome composition differs according to the physiological nature of the skin site, which can be oily, moist or dry [3]. Skin diseases are correlated with the alteration of the normal microbiota, which starts with an imbalance in the microbial homeostasis that is termed dysbiosis [4].

Mycobiomics studies the fungal subpopulation of the microbiome. Mycobiome studies are key to determine the exact role of the human-associated fungal communities and their impact on human health state [5]. Fungal pathogenesis is an increasingly significant cause of human morbidity and mortality [6,7]. Mycobiome analysis for precisely reconstructing the fungal community during fungal infection is highly challenging, but with great clinical value for detection and curation of fungal pathogenesis [8].

While bacteriome studies largely rely on 16S rRNA gene profiling to determine the bacterial and archaeal community composition in a given environment or body site [9], mycobiome analysis cannot rely on 18S rRNA amplicon sequencing with the same efficiency. Being eukaryotic microbes, fungi have similar ribosomes to those in human cells, and thus 18S rRNA amplicon data sets may be overwhelmed with human sequences outnumbering fungal sequences. Instead, the internal transcribed spacer areas (ITS1 and ITS2) of the rRNA-encoding genes are used as better markers for fungal profiling [10].

Malassezia is a lipophilic dimorphic fungus known as pityrosporum. It is a commensal skin fungus, but in some conditions, it behaves as a pathogen. Malassezia species reportedly represent 50%–80% of the skin fungal microbiota, and are most common in the sebaceous areas of skin, such as the trunk, face, scalp and back. Malassezia furfur is the main causative species in pityriasis versicolor infections [11,12].

Surprisingly, M. furfur was systematically missed in the early skin metagenomic studies. Here, we curated a specific list of M. furfur ITS and other ribosomal sequences, and used them to reanalyze publicly available skin metagenomic data sets to computationally investigate the distribution of M. furfur and its relative abundance at different skin sites.

Data sets:

We retrieved publicly available skin metagenomic data sets from the National Center for Biotechnology Information (NCBI) Sequence Read Archive. Specifically, we used the human skin metagenome BioProject PRJNA266117 [13], which includes 675 experiments, BioSamples and sequence runs (URL: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA266117/). We selected a subset of 92 biosamples, representing negative controls (clean swabs) in addition to 11 different skin sites (Figure 1) in the human body. The different biosamples can be categorized according to their physiological nature (Table 1) or their occlusion status. The samples in that study were taken from healthy male and female subjects from different age groups (24-53 years old).

M. furfur marker sequences:

We used nine M. furfur marker ribosomal sequences (KC141972.1, KC141971.1, KC141965.1, KC152900.1, KC152898.1, KC415103.1, KC415101.1, KC415099.1, and KC415088.1), previously proposed by Jagielski et al. [14], as a query set (Table 2).

Bioinformatics analysis:

We used the basic local alignment search tool (BLAST) to screen metagenomic sequences for M. furfur marker genes. The stand-alone BLAST+ application [15] was used for nucleotide searches and mapping, with the BLASTN program. The nine marker genes listed above were the query set and the metagenomic data sets were used as “subject sequences” in the BLASTN searches. An expect value (E-value) of 10^-5 was used as a threshold.

Table 1: List of abbreviations of each skin site and its physiological nature:

Site	Symbol	Physiological nature	Occlusion status
Antecubital fossa	Ac	Rarely intermittently moist	Intermittently occluded
Axillary vault	Ax	Moist	Occluded
Back	Ba	Sebaceous	Exposed
Forehead	Fh	Sebaceous	Exposed
Palm	Pa	Rarely intermittently moist	Exposed
Plantar heel	Ph	Dry	Exposed
Retro auricular crease	Ra	Sebaceous	Occluded
Subclavius	Sc	Sebaceous	Exposed
Toe web space	Tw	Moist	Occluded
Umbilicus	Um	Moist	Occluded
Volar forearm	Vf	Dry	Exposed
Negative	Neg	Control (clean swabs)	Not applicable (NA)

Table 2: Malassezia furfur marker sequences [14]:

1.	KC141972.1	18S ribosomal RNA gene, partial sequence; internal transcribed spacer 1, 5.8S ribosomal RNA gene, and internal transcribed spacer 2, complete sequence; and 28S ribosomal RNA gene, partial sequence
2.	KC141971.1	18S ribosomal RNA gene, partial sequence; internal transcribed spacer 1, 5.8S ribosomal RNA gene, and internal transcribed spacer 2, complete sequence; and 28S ribosomal RNA gene, partial sequence
3.	KC141965.1	18S ribosomal RNA gene, partial sequence; internal transcribed spacer 1, 5.8S ribosomal RNA gene, and internal transcribed spacer 2, complete sequence; and 28S ribosomal RNA gene, partial sequence
4.	KC152900.1	18S ribosomal RNA gene, partial sequence; internal transcribed spacer 1, 5.8S ribosomal RNA gene, and internal transcribed spacer 2, complete sequence; and 28S ribosomal RNA gene, partial sequence
5.	KC152898.1	18S ribosomal RNA gene, partial sequence; internal transcribed spacer 1, 5.8S ribosomal RNA gene, and internal transcribed spacer 2, complete sequence; and 28S ribosomal RNA gene, partial sequence
6.	KC415103.1	26S ribosomal RNA gene, partial sequence
7.	KC415101.1	26S ribosomal RNA gene, partial sequence
8.	KC415099.1	26S ribosomal RNA gene, partial sequence
9.	KC415088.1	26S ribosomal RNA gene, partial sequence

Normalization and metrics:

BLASTN results, i.e., sequences matching any or all of the M. furfur markers, were tabulated and counted. To estimate the distribution of M. furfur marker sequences in a given metagenomic data set, we normalized the total number of BLASTN hits to the number of reads in each metagenomic sample. For a better presentation of the data, we multiplied the normalized hits by million (10⁶) to calculate the fraction per million sequences (FPM).

For optimization and comparative studies, we calculated different metrics, which we later aimed at comparing: (i) total number of hits to any marker gene; (ii) total number of hits to each marker gene; (iii) number of best hits to any marker gene; (iv) number of best hits to each marker gene. FPMs were computed for all the above metrics (Figure 2) and optimized as indicated in the Results section.

Statistical analysis:

Statistical analysis was performed in the R environment [16], installed from URL: https://www.R-project.org, with the RStudio [17] interface (URL: http://www.rstudio.com), with specific use of the data analysis packages: ggplot2, corrplot, and beanplot. In some instances, the GraphPad Prism program (GraphPad Software Inc., San Diego, CA) was used for confirmation.

To investigate the distribution of M. furfur in different body sites, we considered the normalized hits FPM as an estimate of M. furfur relative abundance. We statistically analyzed the relative abundance of M. furfur in skin samples with different physiological conditions (Figure 1). We also explored the effect of sample occlusion, subject age and subject gender on the relative abundance of M. furfur in the samples.

Optimization of metrics:

We started by testing each of the nine biomarkers to assess their diagnostic ability, reproducibility, and thus reliability in detecting M. furfur sequences within metagenomes.

The results obtained from aligning each of the nine markers, and their combination, as BLAST query sequences were quite reproducible (Figure 3). Pearson correlation analysis of number of BLAST hits, represented as normalized hits FPM, indicated a high correlation near 1 in most cases (Pearson correlation coefficients > 0.9, reaching 10.0 at instances). Two sequences with accession number KC141972.1, KC141971.1 gave exactly the same results (Figure 3).

Since the results obtained from normalized hits (FPM) for all genes were quite correlated with results from those of single genes (Figure 2 and Figure 3), as well as with best hits for all genes, we estimated the abundance as all gene FPM, as these values were comparable for small or large data sets, since they were normalized for the number of sequence reads.

Effect of skin site on the distribution of M. furfur in representative skin metagenomes (Figure 4A)

The type/location of the skin site significantly affected the relative abundance of M. furfur among the samples (Kruskal-Wallis rank, p-value = 1.261x10^-06).

Among the 11 skin sites (Figure 1 and Table 1) included in our analysis, the retro auricular crease, which is a sebaceous area behind the ear, had the highest relative abundance of M. furfur, followed by the antecubital fossa, a moist joint area, then the forehead, a sebaceous and exposed area. Next came the palm, which is an intermittently moist and exposed area (Figure 4A).

The back and the subclavius area, which are both sebaceous areas that are covered most of the time had intermediate abundance. Finally, the axilla, plantar heel, toe web space, umbilicus and volar forearm areas had the lowest M. furfur abundance. Interestingly, all these sites are neither sebaceous nor in a joint area and they are covered by clothes most of the time.

Expectedly and reassuringly, the control samples had the lowest number of hits, confirming the specificity of the methodology. These were included in the published study to normalize the reading errors in the results.

Effect of skin nature on M. furfur abundance:

The relative abundance of M. furfur was significantly higher in the sebaceous skin sites, followed by the rarely intermittently moist sites, then the moist sites, and finally the dry sites (Kruskal-Wallis test, p-value = 1.523x10^-08, Figure 4B).

From another perspective, whether the skin is occluded or exposed significantly affected the relative abundance of M. furfur (Kruskal-Wallis rank, p-value = 0.0004654). The highest abundance was observed in the exposed sites, followed by occluded then intermittently occluded (Figure 4C).

Sex and age had no significant effect on the relative abundance of M. furfur among different skin sites (Kruskal-Wallis rank, p-value = 0.6458, 0.1498, respectively) (Figure 5).

Malassezia species are the causative agents of many skin diseases, and their dysbiosis has been linked with many skin disorders and complications of other infections, e.g., HIV [18]. Some researchers reported that Malassezia species may have a role in skin cancer [19] or to the progression of internal organs cancers (esophageal, gut, pancreatic and prostate); in particular, M. furfur was related to the prognosis of prostate cancer [20,21].

The metagenomic study we chose to analyze here (BioProject PRJNA266117) was mainly concerned with the double stranded-DNA virome and its relationship with the human skin microbiota, and focused to a lesser extent on the bacteriome and mycobiome [13]. In terms of fungi, that study only reported the presence of M. globosa but overlooked other Malassezia species, especially M. furfur, which has higher clinical value than M. globosa [13]. This lack of detection of M. furfur is not surprising and is most likely attributed to the absence of a well-developed molecular biomarker to differentiate between Malassezia species at the time of the study (before 2015), given that the earliest NCBI record of a reference ITS sequence for M. furfur (NCBI accession: NR_149347.1) was released in 2016 (URL: https://www.ncbi.nlm.nih.gov/nuccore/NR_149347.1) and that the first complete reference genome for this fungus (NCBI accession: CP046234.1) has been submitted in 2019 (URL: https://www.ncbi.nlm.nih.gov/nuccore/CP046234.1) and was only recently published. [22].

Our results showed that M. furfur was present in a reasonable abundance in several samples. We used the whole genome of each M. furfur and M. globosa and aligned them with skin metagenome sample SRR1647025, one that had the highest number of hits in our study. The number of hits to M. furfur and M. globosa was 8,955 and 14,386, respectively, which indicated the presence of both M. furfur and M. globose, and the ability of our analysis to differentiate between them.

Because the marker set is rather small in number (nine sequences, Table 2), which would strongly affect E-values, underestimating true positive results (i.e., the E-values would be exaggeratedly high if the BLAST database is as small as nine short sequences), we used the marker sequences as a BLASTN query to screen metagenomes (which were used as a search space or ‘BLAST subject sequences’). Either way, we examined all alignments visually to make sure the results were valid and reliable.

Based on the above findings, we conclude that the physiological nature and occlusion of skin site (as exposer to air and joint bending location) are significant factors affecting M. furfur distribution and its relative abundance. On the other hand, sex and age have a minor effect. We can also conclude that alteration of the skin condition may affect the transformation of M. furfur from commensal to pathogenic; we need more investigations to see if this could be applied for other Malassezia species and the whole fungal comminutes in the skin as well.

Ethics statement:

This study is merely computational and involves no patient records or experimental work on human subjects or laboratory animals. All raw data used and analyzed in this study are publicly available.

Competing interests:

The authors declare no financial or personal competing interests related to this work.

Author contributions:

RKA conceived the study; MTA, MH, and RKA designed the experiments; MTA and RKA collected data and metadata; MTA, MH, and RKA analyzed the data; MH and DAE supervised MTA; MTA prepared the figures; MTA and MH drafted the manuscript; DAE and RKA revised the manuscript; All authors read and agreed with the final format.

Gawkrodger, D.J.; Ardern-Jones, M.R. Dermatology: an illustrated colour text; Elsevier: 2016.
ElRakaiby, M.; Dutilh, B.E.; Rizkallah, M.R.; Boleij, A.; Cole, J.N.; Aziz, R.K. Pharmacomicrobiomics: the impact of human microbiome variations on systems pharmacology and personalized therapeutics. OMICS 2014, 18, 402-414, doi:10.1089/omi.2014.0018.
Byrd, A.L.; Belkaid, Y.; Segre, J.A. The human skin microbiome. Nat. Rev. Microbiol. 2018, 16, 143-155, doi:10.1038/nrmicro.2017.157.
Sanford, J.A.; Gallo, R.L. Functions of the skin microbiota in health and disease. Seminars in immunology 2013, 25, 370-377, doi:10.1016/j.smim.2013.09.005.
Cui, L.; Morris, A.; Ghedin, E. The human mycobiome in health and disease. Genome Med. 2013, 5, 63, doi:10.1186/gm467.
Severo, C.B.; Severo, L.C.; Fischer, G.B. Editorial. Fungal infections have emerged as major causes of morbidity and mortality in the last decade. Paediatr. Respir. Rev. 2009, 10, 159-160, doi:10.1016/j.prrv.2009.06.011.
Pal, M. Morbidity and mortality due to fungal infections. J. Appl. Microbiol. Biochem. 2017, 1, 2, doi:10.21767/2576-1412.100002.
McTaggart, L.R.; Copeland, J.K.; Surendra, A.; Wang, P.W.; Husain, S.; Coburn, B.; Guttman, D.S.; Kus, J.V. Mycobiome sequencing and analysis applied to fungal community profiling of the lower respiratory tract during fungal pathogenesis. Front. Microbiol. 2019, 10, 512, doi:10.3389/fmicb.2019.00512.
Aagaard, K.; Petrosino, J.; Keitel, W.; Watson, M.; Katancik, J.; Garcia, N.; Patel, S.; Cutting, M.; Madden, T.; Hamilton, H., et al. The Human Microbiome Project strategy for comprehensive sampling of the human microbiome and why it matters. FASEB j. 2013, 27, 1012-1022, doi:10.1096/fj.12-220806.
Schoch, C.L.; Seifert, K.A.; Huhndorf, S.; Robert, V.; Spouge, J.L.; Levesque, C.A.; Chen, W. Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for fungi. Proc. Natl Acad. Sci. USA 2012, 109, 6241-6246, doi:10.1073/pnas.1117018109.
Goldsmith, L.; Katz, S.; Gilchrest, B.; Paller, A.; Leffel, D.; Wolff, K. Fitzpatrick's Dermatology in General Medicine, 8th Edition ed.; McGraw-Hill Education New York, United States, 2012.
Findley, K.; Grice, E.A. The skin microbiome: a focus on pathogens and their association with skin disease. PLoS Pathog. 2014, 10, e1004436, doi:10.1371/journal.ppat.1004436.
Hannigan, G.D.; Meisel, J.S.; Tyldsley, A.S.; Zheng, Q.; Hodkinson, B.P.; SanMiguel, A.J.; Minot, S.; Bushman, F.D.; Grice, E.A. The human skin double-stranded DNA virome: topographical and temporal diversity, genetic enrichment, and dynamic associations with the host microbiome. mBio 2015, 6, e01578-01515, doi:10.1128/mBio.01578-15.
Jagielski, T.; Rup, E.; Ziółkowska, A.; Roeske, K.; Macura, A.B.; Bielecki, J. Distribution of Malassezia species on the skin of patients with atopic dermatitis, psoriasis, and healthy volunteers assessed by conventional and molecular identification methods. BMC Dermatol. 2014, 14, 3, doi:10.1186/1471-5945-14-3.
Camacho, C.; Coulouris, G.; Avagyan, V.; Ma, N.; Papadopoulos, J.; Bealer, K.; Madden, T.L. BLAST+: architecture and applications. BMC Bioinformatics 2009, 10, 421, doi:10.1186/1471-2105-10-421.
R Core Team R: A language and environment for statistical computing., R Foundation for Statistical Computing: Vienna, Austria., 2020.
RStudio Team RStudio: Integrated Development for R., RStudio, Inc.: Boston, MA., 2019.
Ramos, E.S.M.; Lima, C.M.; Schechtman, R.C.; Trope, B.M.; Carneiro, S. Superficial mycoses in immunodepressed patients (AIDS). Clin. Dermatol. 2010, 28, 217-225, doi:10.1016/j.clindermatol.2009.12.008.
Gaitanis, G.; Magiatis, P.; Hantschke, M.; Bassukas, I.D.; Velegraki, A. The Malassezia genus in skin and systemic diseases. Clin. Microbiol. Rev. 2012, 25, 106-141, doi:10.1128/cmr.00021-11.
Filip, R.; Davicino, R.; Anesini, C. Antifungal activity of the aqueous extract of Ilex paraguariensis against Malassezia furfur. Phytother. Res. 2010, 24, 715-719, doi:10.1002/ptr.3004.
Leake, I. PD-1 inhibitors for oesophageal cancer. Nat. Rev. Gastroenterol. Hepatol. 2019, 16, 706, doi:10.1038/s41575-019-0237-4.
Sankaranarayanan, S.R.; Ianiri, G.; Coelho, M.A.; Reza, M.H.; Thimmappa, B.C.; Ganguly, P.; Vadnala, R.N.; Sun, S.; Siddharthan, R.; Tellgren-Roth, C., et al. Loss of centromere function drives karyotype evolution in closely related Malassezia species. Elife 2020, 9, doi:10.7554/eLife.53944.

Download PDF

Version 1

posted

You are reading this latest preprint version

Why we missed it? Computational analysis reveals distribution patterns of Malassezia furfur, the etiological agent of Pityriasis versicolor, in skin metagenomes

Status:

Version 1

Abstract

Figures

Introduction

Materials And Methods

Results

Discussion

Conclusion

Declarations

References

Status:

Version 1