Studies of eukaryotes in the gut have been limited in part by available methods and databases. Eukaryotes can be identified by the variable regions of their 18S rRNA gene, just as prokaryotes can be identified by 16S sequencing. Based on previous findings of dysbiosis and inflammation in ME/CFS, we hypothesized that the gut eukaryotic community may be altered in ME/CFS and could provide additional clues and potential diagnostic markers for this disease. To characterize gut eukaryote composition and investigate a possible link between gut eukaryotes and ME/CFS, we sought to identify eukaryotes present in the gut of both ME/CFS patients and healthy individuals. To accomplish this, we amplified and sequenced the V9 region of the 18S rRNA gene in both patients and healthy controls and analyzed the results for taxa and diversity differences.

ME/CFS has been previously linked to a possible eukaryotic pathogen. In particular, giardiasis was implicated in ME/CFS outbreaks in New Zealand in 1984 and California in 1985 ( Levine et al., 1992 ; Snow, 2002 ). Following an outbreak of giardiasis in Norway in 2004, at least 5% of infected individuals developed ME/CFS ( Naess et al., 2012 ). A link between ME/CFS and a protozoan infection such as Blastocystis or Dientamoeba has been suggested ( Dunwell, 2013 ). An investigation of Candida albicans in stool samples found increased abundance in ME/CFS patients during the acute phase of the illness versus remission ( Evengard et al., 2007 ). A recent study using culture-based methods investigated total yeast abundance in ME/CFS stool samples and found a nonsignificant decrease compared to healthy individuals ( Armstrong et al., 2017 ).

Myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) is a disease without identified causes or mechanisms, characterized by symptoms of profound fatigue, pain, cognitive impairment, post-exertional malaise and inflammation ( Institute of Medicine, 2015 ). Patients are commonly diagnosed with ME/CFS following a flu-like illness, implicating an inciting pathogen in the disease. In some cases, patients became ill following a diagnosed viral, bacterial or protozoan infection. Patients often complain of gastrointestinal symptoms, such as stomach pain, diarrhea, and nausea. Many patients are diagnosed with gastrointestinal disorders such as irritable bowel syndrome (IBS) ( Institute of Medicine, 2015 ; Maes et al., 2014 ).

Sequences were aligned in QIIME using PyNAST and a minimum identity of 60%. The alignment was filtered with a gap filter threshold of 0.8 and an entropy threshold of 0.1. A phylogenetic tree was then constructed at 60% similarity. Alpha and Beta diversity analyses were conducted at a depth of 100 sequences and all samples with sequence counts below 100 were excluded, resulting in 10 controls and 11 patients. Alpha metrics included phylogenetic diversity, the number of observed species and Shannon and Chao 1 indices. Beta diversity was evaluated using both weighted and unweighted UniFrac distances. Statistical significance of diversity comparisons was determined in QIIME using a Wilcoxon rank-sum test.

To compare specific abundances between patients and controls, taxa were considered present in an individual sample only when represented by at least 5 OTUs. Following removal of OTUs represented by fewer than five sequences in a given sample, data remained from 17 ME/CFS patients and 17 healthy controls. Statistical comparisons of taxa abundance were conducted in both QIIME and in R. Normality of data was checked using a Kolmogorov–Smirnov test. A Wilcoxon rank-sum test was used to test differences in taxa abundance between ME/CFS patients and healthy controls.

Sequencing data was analyzed with Quantitative Insights into Microbial Ecology (QIIME) software (1.9.1). Due to poor quality in the second read, analysis was carried out solely using the data from read one. Sequences were split by barcode, the quality threshold was set at 25 and the reverse primer and downstream sequence were removed. Operational Taxonomic Units (OTUs) were picked using the open-reference method at a similarity of 97%. Taxonomy was assigned using BLAST against the SILVA 119 database with a minimum OTU size of 5 and an e-value of e−50. Taxa corresponding to archaea, bacteria, plant or mammal were filtered from the dataset, as were any OTUs identified in negative controls. 19 healthy controls and 23 ME/CFS patient samples contained OTUs and were retained for further analysis.

Genomic DNA was isolated from stool samples of 39 healthy controls and 49 ME/CFS patients according to a modified protocol from Iliev et al. (2012) . Briefly, stool pellets were resuspended in lyticase and lyticase buffer and incubated at 30 °C for 30 min, then spun at 4,000 rpm for 5 min. After removing the supernatant, pellets were re-suspended in stool DNA stabilizer (B Bridge International), homogenized on a bead beater with 0.1 mm and 0.5 mm beads, heated at 95 °C for 10 min, vortexed for 5 s, cooled on ice for 5 min, and spun at 15,000 rpm for 1 min. The remaining extraction steps were conducted with the QIAmp DNA Mini Kit (Qiagen, Valencia, CA, USA). The V9 region of the 18S rRNA gene was amplified using Illumina_Euk_1391f universal forward primer and individually barcoded Illumina EukBr reverse primers, according to the Earth Microbiome protocol ( Earth Microbiome Project, 2017 ). Products were checked via gel electrophoresis for an approximate size of 200 base pairs and purified using Mag-Bind EZ-pure beads (Omega Bio-Tek, Norcross, GA, USA). DNA concentration was determined using QuantIT PicoGreen dsDNA kit (Invitrogen, Carlsbad, CA, USA) and 100 ng of each product was pooled for sequencing. Subjects with concentration too low to pool 100 ng were eliminated, leaving 29 controls and 35 ME/CFS patients. Five water controls from PCR amplification were included in sequencing to serve as negative controls. Pooled products were sequenced in a 2 × 250 bp paired-end run on the Illumina MiSeq platform at the Cornell Biotechnology Resource Center Genomics Facility. The sequence data will be available in the European Nucleotide Archive of the European Bioinformatics Institute under accession number PRJEB23827 .

ME/CFS patients and healthy controls were recruited by Dr. Susan Levine, M.D., in New York City. Patients were diagnosed with ME/CFS under the Fukuda Criteria ( Fukuda et al., 1994 ). Additional healthy controls were recruited at Ithaca College by Dr. Betsy Keller. Individuals with a known acute illness or chronic infectious disease were excluded from the study. Written consent was obtained from all participants. Stool samples were collected by subjects at home and refrigerated in RNAlater. Following overnight shipment to Cornell University, the samples were stored in aliquots at −80 °C. For all subjects, BMI, age, sex, and gastrointestinal symptom information was collected. Patients also completed Bell’s disability scale ( Bell, 1994 ). All protocols involving human subjects were approved by the Cornell University Institutional Review Board (#1012001855).

Results

Sample population DNA was extracted from the stool samples of 49 ME/CFS patients and 39 healthy individuals. After amplification of the 18S rRNA gene, sequencing and processing, we identified eukaryotic sequences in 17 patients with ME/CFS and 17 healthy individuals for further analysis. Within the healthy population, there were 16 females and one male, while within the ME/CFS population there were 13 females and four males (Table 1). The average age was similar between both groups at 44.6 ± 10.9 in controls and 52 ± 11.9 in patients (Table 1). Average Body Mass Index (BMI) was nearly identical at 27.4 ± 4.5 in controls and 26.8 ± 4.7 in patients (Table 1). 11 of 17 ME/CFS patients reported gastrointestinal symptoms, compared to six of 17 healthy individuals (Table 1). We additionally had Bell Scale ratings for 15 patients, with scores ranging from 10 to 50 out of 100 in all individuals (Table 1). Taxa abundance comparisons Diversity subgroup Controls (n = 17) ME/CFS (n = 17) Controls (n = 10) ME/CFS (n = 11) Males 1 4 1 4 Females 16 13 9 7 Age 44.6 ± 10.9 52 ± 11.9 43.6 ±13.6 54 ± 11.6 BMI 27.4 ± 4.5 26.8 ± 4.7 27.7 ± 5.2 26.2 ± 4.9 Gastrointestinal symptoms 6 11 3 8 Bell Scale 10 1 0 20 4 3 30 2 1 40 5 3 50 3 2 N/A 17 2 10 2 DOI: 10.7717/peerj.4282/table-1

Sequencing results and OTU picking Sequencing yielded 12,430,061 raw reads. After quality control and library splitting, a total of 3,401,011 reads remained. Open-reference OTU picking resulted in 2015 unique OTUs with a total of 3,384,926 OTUs. After removing bacteria, archaea, plant and mammal OTUs from further analysis, as well as any OTUs found in negative controls, 484 unique OTUs and 68,214 total OTUs remained. Additional OTUs that were not found in the constructed phylogenetic tree or were determined to have errors in classification were removed, resulting in 64,760 total OTUs and 310 unique OTUs. Finally, any OTUs not represented by a minimum of 5 sequences in an individual were removed for the purpose of comparing relative taxa abundance. This left a total of 34 individuals (17 patients and 17 controls) for taxa analysis, as well as 310 unique OTUs and 63,847 total OTUs. The average OTUs per individual was 1,878 and the standard deviation was 6,435.

Identified taxa On average, fungi were the most abundant taxa identified in our population (32.1%). Within fungi, the phylum Ascomycota was the most abundant (17.3%), followed by Basidiomycota (13.3%) and Zygomycota (0.9%). An average of 0.7% were fungi, but belonged either to another phylum or were unclassified at the phyla level. For simplicity, these were grouped as “Other Fungi” (Figs. 1 and 2). A small proportion of OTUs were classified as Stramenopiles (6.6%) and one ME/CFS patient contained 0.8% Excavata (Figs. 1 and 2). An average of 61.4% of OTUs per individual were unclassified at any level after BLAST against the SILVA 119 database and termed “Unknown” (Figs. 1 and 2). Additional OTUs could not be classified to consensus at lower levels of taxonomy. Thus, in many cases, analysis was restricted to higher levels of phylogeny. Figure 1: Individual relative taxa abundances in healthy controls and ME/CFS patients. (A) Relative abundance of eukaryotes observed at the phyla level for 17 healthy controls and (B) 17 ME/CFS patients. Each bar represents one individual. Abundances are expressed as percent of total observed OTUs in an individual. Taxa are considered present if represented by a minimum of 5 OTUs. Unknown refers to OTUs with no match in SILVA 119 database. (C) Relative abundance of eukaryotes observed at the level of order for 17 healthy controls and (D) 17 ME/CFS patients. In cases where OTUs were not classified at the level of order, names indicate lowest identified taxa level followed by unknown. Figure 2: Average relative taxa abundances in healthy controls and ME/CFS patients. (A) Average relative abundance of taxa at the level of phyla in healthy controls versus ME/CFS patients. (B) Average relative abundance of taxa at the level of order in healthy controls versus ME/CFS patients. Nearly all of the fungal OTUs identified as Ascomycota were classified in the Saccharomycetales order (97%). Additional observed Ascomycota were members of class Euromycetes and the orders Pleosporales, Chaetothyriales, Eurotiales and Helotiales. Within Basidiomycota, the majority of OTUs belonged to the Agaricomycetes class (83.7%), while much of the remainder belonged to Tremellomycetes (9.4%), followed by, Microbotryomycetes, Exobasidiomycetes, and Ustilaginomycetes. Most OTUs within Agaricomycetes were members of the Agaricales order, but others belonged to the orders Boletales, Polyporales, and Trechisporales. Similarly, most OTUs from Tremellomycetes were from order Tremellales, but one individual also had OTUs from Cystofilobasidiales (Fig. 1D). Zygomycota consisted mainly of order Mucorales (93.9%), with a small proportion of Entomophthorales (6.1%). Among Stramenopiles, almost all OTUs were of the Blastocystida order (74.9%). Less abundant Stramenopiles belonged to Peronosporales (14.1%), Pleurosigma (2.1%), and Eustigmatales (9%). Analysis of differences in relative abundance at lower levels of taxonomy was limited due to many unclassified reads. In particular, many of the reads assigned to the Saccharomycetales order were not classified further.