Decreased microbial abundance in MGD meibum and eyelid skin
In all, we enrolled 61 MGD patients and 15 volunteers as healthy controls (HC) (Fig. 1, Table S1). They provided DNA samples from swabs of meibum from the eyelid skin and conjunctiva following our previous procedure (Fig. 1, Method). After performing quality control and removing host sequences, 117 metagenome datasets were obtained. Each dataset contained over 100 K sequencing reads for metagenome assembly and subsequent annotations (Table S2). In taxonomic profiling, taxa were determined based on nucleic acid annotation and marker gene presence, respectively (Figure S1). In de novo assembly, each dataset had 32,692 contigs on average, with their mean length ranging from 225 bp to 528 bp (Table S3).
In MGD meibum, microbial populations having the same phylum were preponderant at the three ocular sites: Proteobacteria, Actinobacteria, Firmicutes, Bacteroidetes, and Negarnaviricota (Fig. 2A). At the genus level, the predominant genera at the three ocular sites are Pseudomonas, Cutibacterium, Campylobacter, Corynebacterium, Rubrobacter. At the species level, Cutibacterium acnes, Pseudomonas azotoformans, Rubrobacter xylanophilus, Campylobacter coli, Pseudomonas fluorescens were preponderant (Fig. 2A). Samples were clustered into distinct groups according to different levels of taxonomic classification and their disease status (Fig. 2B, Figure S1). As expected, the eyelid skin microbiome had the highest community richness, whereas conjunctival flora was sparse. In MGD patients, their meibum microbiome community population was much smaller than in the HC group (chao1, Fig. 2C), while their community diversities were similar (Shannon, Fig. 2D; Simpson, Fig. 2E). This population decline was attributable to decreases in the abundance of the more preponderant microbes in these patients. Interestingly, alpha diversity was similar in the meibum and eyelid skin microbiomes (Fig. 2C-E).
Changes in the microbial community were observed in the microbiomes of MGD meibum and eyelid skin at different taxonomic classification levels (Fig. 3A-C). Comparing to HC the most significant microbial change at the phylum level in MGD meibum was the significantly decreasing abundance of Proteobacteria (log2FoldChange = -4, adjusted p-value = 2.8 × 10− 15, Table S4), which was also seen in the MGD eyelid skin microbiome (Fig. 3D, Table S5). A comparison of the MGD and HC groups shows that the top ten genera whose abundance was not the same, and they were present at higher levels in MGD meibum. These genera included Rubrobacter, Novibacillus, Campylobacter, Geobacillus, Sphingomonas, Corynebacterium, Sphingobium, Pedobacter, Fictibacillus, and Enterococcus (log2FoldChange = 3 ~ 14, adjusted p-value < 1 × 10− 6, Table S4). Among these genera, the same changes were observed in the top three and Corynebacterium in the eyelid skin microbiome, while greater abundance of Geobacillus, Sphingomonas, Sphingobium, Pedobacter, Fictibacillus, and Enterococcus was specific to the microbiome of the MGD meibum. At the species level, the top ten were Rubrobacter xylanophilus, Novibacillus thermophilus, Sphingomonas sp. SL9, Campylobacter coli, Campylobacter jejuni, Sphingomonas panacis, Sphingomonas hengshuiensis, Rubrobacter radiotolerans, Sphingobium sp. SYK-6, Sphingomonas wittichii, all of them were more abundant in MGD samples (log2FoldChange = 5 ~ 14, adjusted p-value < 1 × 10− 9, Table S3). Among them, only Sphingobium sp. SYK-6 was only more abundant in the MGD meibum whereas the other 9 species were also more preponderant in the microbiomes of the MGD eyelid skin (Fig. 3F).
At the genus and species levels, there were other changes that were highly significant (adjusted p-value < 0.005, Benjamini-Hochberg, Fig. 3E-F) in the MGD meibum (genus,47; species,185) than in the eyelid skin (genus, 31; species, 99). Specifically, all of the taxa that were common to the two aforementioned sites were more abundant in MGD patients than HC volunteers. However, in some strains of the Pseudomonas genus there were species with altered abundance that were specific to MGD meibum. These strains account for 88% (44/50, green, Fig. 3F) of the species with decreased abundance in the MGD meibum at this significance level. Nevertheless, two other species of Pseudomonas, P. yamanorum and P. virus SM1 (bold italic, Fig. 3F), are exceptions, whose abundance in meibum was greater in MGD than in HC.
Pathogen Preponderance In Meibum And Disease Status
In HC samples, 968 different strains of pathogens were detected (Method), whereas MGD samples contained 2400 pathogens (pathogen, Table S1). However, on the average, the population of each individual pathogen in the MGD group at the two sites (conjunctiva and eyelid skin) was smaller than in the HC group (Mean, MGD 26 vs HC 37, p = 0.02, Welch t-test). Interestingly, the MGD meibum also had far fewer pathogens than that in the HC group (Mean, MGD 13 vs HC 36, p = 0.0014, Welch t-test, Figure S2A). In addition, in the MGD patients, the magnitude of decreases in the number of different types of pathogens were slightly correlated with increases in disease severity (p = 0.065, Welch t-test, Figure S2B). Additionally, there was either no gender or age difference between the numbers of pathogens in the meibum and those at two ocular sites (conjunctiva and eyelid skin) between the MGD and HC groups (Figure S2C and Figure S3).
Twenty pathogens were identified whose positive rate was different in MGD patients from that in the HC (p < 0.05, Welch t-test, Table 1, Figure S2D). Furthermore, 28 pathogens were detected whose positive rate was greater than 10% in the meibum samples (Table 2, Figure S2E). The most preponderant pathogens in all meibum samples were Pseudomonas fluorescens, in over 90% of the samples (Table 2, Figure S2E), Cupriavidus metallidurans (76%) and Pseudomonas putida (52%). Nevertheless, only Pseudomonas fluorescens was present in 90% of the MGD and 100% of the HC samples (p = 0.044, Table 2). Pathogens whose positive rate was much higher in MGD than in HC meibum were Campylobacter coli (p = 4.4 × 10− 10, Welch t-test, Table 1), Campylobacter jejuni (p = 1.2 × 10− 8, Welch t-test), and Enterococcus faecium (p = 8.3 × 10− 8, Welch t-test). Furthermore, their abundance was more than 16-fold higher in the MGD than in the HC meibum (adjusted p-value < 2 × 10− 8, Table S4). Additionally, Pseudomonas aeruginosa (p = 4.8 × 10− 6), Pseudomonas mosselii (p = 1.4 × 10− 4, Welch t-test), Escherichia coli (p = 0.001, Welch t-test), Stenotrophomonas maltophilia (p = 0.008, Welch t-test), and Neisseria sicca (p = 0.016, Welch t-test) had a lower positive rate in the MGD meibum and samples from the eyelid and conjunctiva.
Table 1
Pathogens with differential positive rate in meibum between MGD and HC
Pathogen | MGD (n = 47*) | HC (n = 11*) | P-value |
Pseudomonas fluorescens (G−) | 42 (89%) | 11 (100%) | 0.024 |
Campylobacter coli (G−) | 27 (57%) | 0 (0%) | 4.4E-10 |
Campylobacter jejuni (G−) | 24 (51%) | 0 (0%) | 1.2E-08 |
Enterococcus faecium (G+) | 22 (47%) | 0 (0%) | 8.3E-08 |
Malassezia globosa | 22 (47%) | 10 (91%) | 0.001 |
Pseudomonas protegens (G−) | 18 (38%) | 8 (73%) | 0.045 |
Escherichia coli (G−) | 12 (26%) | 9 (82%) | 0.001 |
Ralstonia pickettii (G−) | 12 (26%) | 8 (73%) | 0.008 |
Pseudomonas aeruginosa (G−) | 11 (23%) | 10 (91%) | 4.8E-06 |
Pseudomonas mosselii (G−) | 6 (13%) | 9 (82%) | 1.4E-04 |
Bacillus licheniformis (G+) | 6 (13%) | 0 (0%) | 0.013 |
Yarrowia lipolytica | 5 (11%) | 0 (0%) | 0.024 |
Enterobacter cloacae complex (G−) | 3 (6%) | 5 (45%) | 0.034 |
Comamonas testosteroni (G−) | 2 (4%) | 5 (45%) | 0.026 |
Serratia marcescens (G-) | 2 (4%) | 5 (45%) | 0.026 |
Stenotrophomonas maltophilia (G−) | 1 (2%) | 6 (55%) | 0.008 |
Neisseria sicca (G−) | 0 (0%) | 5 (45%) | 0.016 |
Neisseria meningitidis (G−) | 0 (0%) | 4 (36%) | 0.038 |
Roseomonas gilardii (G−) | 0 (0%) | 4 (36%) | 0.038 |
Staphylococcus capitis (G+) | 0 (0%) | 4 (36%) | 0.038 |
* with meibum samples provided. MGD, HC, healthy control; G+,Gram-positive;G−༌Gram-negative. |
Table 2
Pathogens with positive rate over 10% across all meibum samples
Pathogen | #Total meibum samples (n = 58) | MGD (n = 47) | HC (n = 11) | P-value |
Pseudomonas fluorescens (G−) | 53 (91%) | 42 (89%) | 11 (100%) | 0.02 |
Cupriavidus metallidurans (G−) | 44 (76%) | 36 (77%) | 8 (73%) | 0.81 |
Malassezia globosa | 32 (55%) | 22 (47%) | 10 (91%) | 8.9E-04 |
Pseudomonas putida (G−) | 30 (52%) | 22 (47%) | 8 (73%) | 0.12 |
Campylobacter coli (G−) | 27 (47%) | 27 (57%) | 0 (0%) | 4.4E-10 |
Pseudomonas protegens (G−) | 26 (45%) | 18 (38%) | 8 (73%) | 0.05 |
Campylobacter jejuni (G−) | 24 (41%) | 24 (51%) | 0 (0%) | 1.2E-08 |
Enterococcus faecium (G+) | 22 (38%) | 22 (47%) | 0 (0%) | 8.3E-08 |
Corynebacterium jeikeium (G+) | 22 (38%) | 20 (43%) | 2 (18%) | 0.10 |
Alternaria alternata | 22 (38%) | 16 (34%) | 6 (55%) | 0.25 |
Escherichia coli (G−) | 21 (36%) | 12 (26%) | 9 (82%) | 8.6E-04 |
Pseudomonas aeruginosa (G−) | 21 (36%) | 11 (23%) | 10 (91%) | 4.8E-06 |
Citrobacter freundii (G−) | 20 (34%) | 14 (30%) | 6 (55%) | 0.17 |
Ralstonia pickettii (G−) | 20 (34%) | 12 (26%) | 8 (73%) | 0.01 |
Staphylococcus epidermidis(G+) | 19 (33%) | 13 (28%) | 6 (55%) | 0.14 |
Corynebacterium simulans (G+) | 18 (31%) | 15 (32%) | 3 (27%) | 0.77 |
Corynebacterium striatum (G+) | 18 (31%) | 15 (32%) | 3 (27%) | 0.77 |
Corynebacterium aurimucosum (G+) | 16 (28%) | 13 (28%) | 3 (27%) | 0.98 |
Klebsiella pneumoniae (G−) | 15 (26%) | 10 (21%) | 5 (45%) | 0.18 |
Pseudomonas mosselii (G−) | 15 (26%) | 6 (13%) | 9 (82%) | 1.4E-04 |
Corynebacterium diphtheriae (G+) | 14 (24%) | 11 (23%) | 3 (27%) | 0.81 |
Corynebacterium ureicelerivorans (G+) | 11 (19%) | 7 (15%) | 4 (36%) | 0.21 |
Toxoplasma gondii | 10 (17%) | 7 (15%) | 3 (27%) | 0.43 |
Bacteroides vulgatus (G−) | 10 (17%) | 5 (11%) | 5 (45%) | 0.06 |
Aeromonas hydrophila (G−) | 9 (16%) | 6 (13%) | 3 (27%) | 0.35 |
Corynebacterium riegelii (G+) | 7 (12%) | 5 (11%) | 2 (18%) | 0.57 |
Corynebacterium urealyticum (G+) | 7 (12%) | 5 (11%) | 2 (18%) | 0.57 |
Bacillus licheniformis (G+) | 6 (10%) | 6 (13%) | 0 (0%) | 0.01 |
G+,Gram-positive;G−༌Gram-negative. | | | | |
Metabolic features of microbiome of MGD meibum and eyelid skin |
The assignation results for function annotation in the KEGG list includes the over-represented gene sets at different classification levels in the MGD meibum and eyelid skin. On the other hand, no differential entry was found for conjunctiva (adjusted p-value > 0.05, Table S6-S7). Similar to taxonomy profiling, the top 10 gene sets with the largest fold abundance change in the MGD meibum microbiome in each annotation category (Fig. 4A-E) also had significant alterations in the MGD eyelid samples, but not in the MGD conjunctival microbiome. The most notable microbiome changes in both the MGD meibum and eyelid skin were those pathways in mediating microbial metabolism. Among the top 10 KEGG pathways, four were involved in the degradation process including xylene (ko00622), dioxin (ko00621), ethylbenzene (ko00642), and bisphenol (ko00363). |
We further explored the molecular processes underlying 509 enzymes, which were mapped into KEGG pathways. Table S7 lists the enzymes in the MGD meibum microbiome whose abundance increased more than two-fold from their levels in the HC meibum (adjusted p-value < 0.05). Over half of them are involved in KEGG Metabolic pathways (63%, 319/509), one hundred of them in biosynthesis of secondary metabolites, 95 others in KEGG pathways of microbial metabolism in diverse environments, and 49 in biosynthesis of antibiotic reference pathways (Table S8). Specifically, all of the 15 pathways involved in carbohydrate metabolism had at least 4 enzymes with increased abundance. The second significant KEGG pathway category is involved with lipid metabolism. Ten out of the 17 pathways contained at least 2 enzymes that were mapped into pathways mediating fatty acid elongation, biosynthesis and degradation, glycerolipid metabolism, glycerophospholipid metabolism, ether lipid metabolism, sphingolipid metabolism, biosynthesis of unsaturated fatty acids, carbon fixation pathways in prokaryotes, and nitrogen metabolism (Table S8). |
For xenobiotic biodegradation and metabolism, 30 enzymes with increased sequencing read abundance mapped into this pathway category. Interestingly, the crucial one was likely mediating benzoate degradation. First, thirteen of these enzymes (43%, 13/30) directly mapped into this pathway (Figure S4A). Second, the subsequent 8 enzymes with increased abundance catalyze the reactions that precede by one or two steps entrance into the benzoate degradation pathway, including E1.14.13.7 (phenol 2-monooxygenase, NADPH) and E1.14.13.82 (vanillate monooxygenase) in aminobenzoate degradation, E3.7.1.8 (2,6-dioxo-6-phenylhexa-3-enoate hydrolase) in dioxin degradation (Figure S4B), E2.3.1.16 (acetyl-CoA C-acyltransferase) in ethylbenzene degradation(Figure S4C), E1.13.11.2 (catechol 2,3-dioxygenase) in styrene degradation (Figure S4D), and E1.14.13.1 (salicylate 1-monooxygenase) in dioxin degradation (Figure S4B), naphthalene degradation (Figure S4E) and polycyclic aromatic hydrocarbon degradation (Figure S4F). In addition, four enzymes involved in the initial step in this process are in a branch leading to benzoate degradation, including E1.1.1.2 (alcohol dehydrogenase, NADP+) and E3.1.1.17 (gluconolactonase) in caprolactam degradation (Figure S4G), E1.14.12.10 (benzoate 1,2-dioxygenase) in fluorobenzoate degradation pathway (Figure S4H) and E4.1.99.11 (benzylsuccinate synthase) in the toluene degradation pathway (Figure S4I). |
Functional features identified in EggNOG database, virulence Factors Database (VFDB) and Antibiotic Resistance Genes Database (ARDB)
EggNOG annotation for contigs obtained in de novo assembly demonstrated significantly more functional entries with decreased abundance in the microbiomes of MGD meibum than those in the MGD eyelid skins (p < 0.0001, chi-square test, Fig. 5A; details in Table S10-11), which is similar with the observation in taxonomy profiling. The featured virulence factor for MGD meibum microbiome was Type IV secretion system (T4SS), achieving about a 5-fold change in abundance (q = 0.017, FDR, Fig. 5B, Table S9). For antibiotic resistance genes, both sites had significantly reduced representation of the FomA gene (Fig. 5C), involved in Fosfomycin resistance.