Occupancy strongly influences faecal microbial composition of wild lemurs

Occupancy strongly influences faecal microbial composition of wild lemurs Abstract The microbiota of the mammalian gut is a complex ecosystem, the composition of which is greatly influenced by host genetics and environmental factors. In this study, we aim to investigate the influence of occupancy (a geographical area of habitation), species, age and sex on intestinal microbiota composition of the three lemur species: Eulemur fulvus, E. rubriventer and E. rufifrons. Faecal samples were collected from a total of 138 wild lemurs across Madagascar, and microbial composition was determined using next-generation sequencing of PCR-amplified 16S rRNA gene fragments. Consistent with reports from other primate species, the predominant phyla were Firmicutes (43 ± 6.4% [s.d.]) and Bacteroidetes (30.3 ± 5.3%). The microbial composition was strongly associated with occupancy in the E. fulvus population, with up to 19.9% of the total variation in microbial composition being explained by this factor. In turn, geographical differences observed in faecal microbiota of sympatric lemur species were less pronounced, as was the impact of the factors sex and age. Our findings showed that among the studied factors occupancy had the strongest influence on intestinal microbiota of congeneric lemur species. This suggests adaptation of microbiota to differences in forest composition, climate variations and correspondingly available diet in different geographical locations of Madagascar. microbiota, Eulemur, Madagascar, gastro-intestinal tract, environment, multivariate statistics INTRODUCTION The intestinal microbiota of mammals is an integral part of an animal's body that contributes significantly to the overall health of the host through modulation of its immune system, facilitation of food digestion, competition with pathogenic microorganisms and production of metabolites beneficial for the host (Kabat, Srinivasan and Maloy 2014; Patterson et al.2014). Expression of these beneficial properties directly correlates with microbial community diversity and composition (Clemente et al.2012). Hence, identifying the factors and underlying processes that shape the intestinal microbiota is important for a better understanding of its contribution to host health. Previous studies in humans have shown that host genetics (Hansen, Gulati and Sartor 2010), lifestyle (David et al.2014) and food preferences (Maukonen and Saarela 2015) contribute to shaping microbiota composition of an individual within a population. Intestinal microbiota composition can be distinguished between different mammalian species, suggesting co-evolution and adaptation of animals and their microbes (Ley et al.2008; Moeller et al.2016). It is not clear, however, to what extent host genotype and environmental factors influence intestinal microbial composition under natural conditions among closely related animal species dwelling in different biogeographical regions. Although wildlife microbiota has received less attention in comparison with that of humans, farm and rodent model animals, data collected from animals in wild conditions can provide complementary information that contributes to our understanding of processes that shape mammalian intestinal microbiota. For instance, studies that highlight similarities and differences in microbiota between humans and other Homininae species (Ellis et al.2013; Moeller et al.2014; Schnorr et al.2014) provided new insight into evolution of microbiota, suggesting adaptation of human microbes to an animal protein-based diet. Studies on microbiota composition of primates that are evolutionarily more distant from humans, such as yellow baboons (Papio cynocephalus) (Ren et al. 2016), black howler monkey (Alouatta pigra) (Amato et al.2015), black and white colobus (Colubus guereza), red colobus (Piliocolobus tephrosceles) and red-tailed guenon (Cercopithecus ascanius) (Yildirim et al.2010) revealed that microbiota composition of these primates is highly variable, also intra-individually, and mostly depends on the available diet. Correspondingly, the diet of a wild animal directly depends on suitable food availability, which depends on climate, flora and fauna of an area. This statement is also true for wild lemurs of Madagascar, for which several studies showed variation in feeding patterns and diet when comparing areas with different forest composition (Balko and Brian Underwood 2005), as well as different seasons (Overdorff 1993; Overdorff, Strait and Telo 1997). Fogel (2015) compared the microbiota composition of sympatric wild Lemur catta and Propithecus verreauxi across dry and wet seasons and showed that microbiota of both lemur species is variable between individuals and dynamic over time. Researchers observed differences in microbial composition between wild vs captive L. catta as well as wild populations of L. catta and P. verreauxi, albeit only with respect to relative abundance of specific microbial groups rather than their presence or absence (McKenney, Rodrigo and Yoder 2015). Wild rufous mouse lemurs (Microcebus rufus) showed an increase in gut microbial diversity with age and differences in microbiota richness and diversity between sampling. Furthermore, microbial composition was affected by site, sex and year, whereas temporal trends within a year were weak (Aivelo, Laakkonen and Jernvall 2016). The above-mentioned studies of lemur microbiota were focused on a single lemur species (Aivelo, Laakkonen and Jernvall 2016), two sympatric lemur species dwelling in the same area (Fogel 2015) or captive lemurs of different species (McKenney, Rodrigo and Yoder 2015). Taken together, these studies showed that lemurs harbour complex intestinal microbiota, the composition of which fluctuates over time among and within individuals, and is affected by season, captivity, age, site of sampling and sex. In our study we focused on microbiota of three closely related Eulemur species (E. fulvus, E. rufifrons and E. rubriventer). To the best of our knowledge, this is the first comparative study of intestinal microbiota composition of multiple wild Eulemur species across Madagascar exposed to large variations in climate conditions and biogeography. In addition to an explorative assessment of the most important features of intestinal microbial composition in the studied species, we addressed to what extent occupancy, host species, sex and age influence lemur intestinal microbial composition, and which of these factors contribute most strongly to intestinal microbiota differentiation in wild lemurs. To this end, we hypothesised that intestinal microbial composition is similar among congeneric lemur species and that of the remaining factors occupancy, with habitat as a determining factor of food items availability and other climate-related effects, has the strongest influence on intestinal microbiota differentiation. MATERIAL AND METHODS Study design Faecal samples (N = 138) were selected from wild lemurs collected across Madagascar from April to July 2014 (Fig. 1). To investigate the effect of the occupancy on intestinal microbiota, E. fulvus samples from three climatic regions and E. rufifrons samples from two climatic regions were compared with each other. To assess the influence of different species, E. rubriventer samples collected in Ranomafana National Park (NP) were compared to E. rufifrons samples from the same area. The effect of age and sex was estimated based on E. rufifrons samples collected in Kirindy NP and Ranomafana NP (Table S1, Supporting Information). Figure 1. View largeDownload slide Lemur faecal sample collection areas across Madagascar. The map showing main types of land cover and vegetation was adapted from www.wildmadagascar.org, and was produced with data taken from the FAO Country Profiles and Mapping Information System (The United Nations Food and Agricultural Organization; © FAO 2004). Faecal samples were collected at five geographical locations across the island from E. rufifrons (two sites), E. fulvus (three sites) and E. rubriventer (one site). n = number of samples. Figure 1. View largeDownload slide Lemur faecal sample collection areas across Madagascar. The map showing main types of land cover and vegetation was adapted from www.wildmadagascar.org, and was produced with data taken from the FAO Country Profiles and Mapping Information System (The United Nations Food and Agricultural Organization; © FAO 2004). Faecal samples were collected at five geographical locations across the island from E. rufifrons (two sites), E. fulvus (three sites) and E. rubriventer (one site). n = number of samples. Study sites Madagascar experiences a strong variation in climate conditions, resulting in different vegetation zones across the island (Irwin et al.2010). The effect of environmental factors on lemur microbiota composition was investigated at five sites across Madagascar (Fig. 1). Kirindy Mitea NP (20°07΄S, 44°67΄E, 722 km2) and Ankarafantsika NP (16°25΄S, 46°80΄E, 1350 km2) consist of dry deciduous forest and are located on the western- and north-western side of Madagascar, respectively (Goodman and Benstead 2005). Kirindy Mitea NP is characterised by pronounced seasonality. The area contains more than 200 species of trees with a mean canopy height of 12–18 m, containing mostly deciduous trees with adaptations to water stress (Lewis and Bannar-Martin 2012). Ankarafantsika NP is a mosaic of floristically heterogeneous dry deciduous forests dissected by small valleys with abundant Raffia palms (Rakotonirina 1996; Ganzhorn and Schmid 1998; García and Goodman 2003; Sorg, Ganzhorn and Kappeler 2003) Ranomafana NP is located in southeastern Madagascar (21°16΄S, 47°20΄E) and encompasses ∼435 km2 of montane moist forest, ranging from altitudes of 500 m up to 1500 m, and receives an average of 3000 mm rainfall per year (Wright et al.2012). The rainfall in Ranomafana NP differs highly between the wet-warm season (December to March, 482–1170 mm per month) and dry-cold season (April to November, 55–513 mm per month) (Atsalis 2000). Andasibe Mantadia NP (155 km2) is located at the eastern side of Madagascar (18°92΄S, 48°42΄E), and is also characterised by relatively wet rain forests. Nosy Tanikely (13°28΄S, 48°14΄E) is an island in the north-east of Madagascar covered with tropical vegetation. This island comprises less than 0.3 km2 and is located between Nosy Be (8 km) and the mainland of Madagascars (13 km). Elevation ranges from 0 to 47 m above sea level (Vences, Köhler and Glaw 1997). The island's vegetation consists of low forest with planted banana and mango trees, surrounded by a sandy shore with large rock formations (de Winter, personal observation) Studied species This study focused on three Eulemur species: the red-fronted lemur (E. rufifrons), the common brown lemur (E. fulvus) and the red-bellied lemur (E. rubriventer). These species are morphologically alike and are frugivorous, although they may include other food sources, such as leaves and invertebrates, in their diet (Sussman 1974; Overdorff 1996; Pyritz, Kappeler and Fichtel 2011; Sato 2013). The main difference in social organisation between the different Eulemur species is their group size. Eulemur rufifrons and E. fulvus live in multimale, multifemale groups from 4 to 18 individuals (Overdorff 1996; Pyritz, Kappeler and Fichtel 2011; Johnson et al.2005), whereas E. rubriventer lives in small monogamous groups from two up to five individuals (Tecot 2008; Tecot 2010). Sampling and data collection Immediately after defecation, fresh faecal samples (3–4 g) from individual lemurs were collected non-invasively. Within 12 h after collection, the samples were stored at ambient temperature in sterile plastic tubes that were prefilled with 5 ml of 70% ethanol until further analyses at the Laboratory of Microbiology, Wageningen University, The Netherlands. Species, age and sex were recorded. All samples included in this paper were taken in compliance with the laws of the Government of Madagascar, and no animal experimentation was involved. DNA extraction Samples collected in the Ranomafana NP were processed using a modified protocol based on method proposed by Yu and Morrison (2004) with modifications described previously (Salonen et al.2010). For this method, faecal material was air-dried for 15–20 min in a fume hood to remove ethanol from the samples. Subsequently, 0.1–0.17 g of dried samples were added into double autoclaved screw-cap tubes containing 0.3 g of 0.1 mm zirconia beads, three pieces of 2.5 mm glass beads and 700 μl of lysis buffer (500 nM NaCl, 50 mM Tris-HCI (pH = 8), 50 mM EDTA, 4% SDS) in each. Samples were treated for 3 × 1 min at 5.5 × 103 movements per minute in a Precellys 24 beadbeater (Bertin technologies, France). After homogenisation, samples were incubated at 95°C for 15 min in a shaking heating block (Vartemp 56, Labnet International, Edison, NJ, USA) at 100 rpm, then centrifuged at 4°C for 5 min at 13 000 rpm. Clean supernatants were transferred into 2 ml tubes. Three hundred microlitres of fresh lysis buffer was added in the same tubes to the pellets, bead beating/incubation steps were repeated and freshly collected supernatant was pooled with that previously collected. Subsequent steps were performed according to the original protocol (Yu and Morrison 2004). Samples collected in Adnasibe NP, Kirindy Mitea NP, Ankarafantsika NP and Nosy Tanikely were extracted using an automatic system, the Maxwell® 16 Research Instrument (Promega, Madison, USA) and the corresponding RNA extraction kit according to the manufacturer's instructions. To improve DNA yield, samples preserved in 70% ethanol were rehydrated through a series of ethanol solutions with decreasing proportions of ethanol in steps of 10%. For rehydration, 1.5 ml of 70% ethanol with faecal particles was transferred into a fresh 2 ml tube and centrifuged at 13 000 rpm for 5 min to separate solid fractions from the liquid. After centrifugation, part of the supernatant was replaced with the same amount of distilled water to decrease ethanol concentration by 10% point, vortexed and incubated for 10 min at RT. These steps were repeated until the ethanol was replaced by distilled water. Cell disruption and lysis was performed as described above, but instead of lysis buffer we used S.T.A.R buffer (Roche Molecular Systems, USA). DNA quality and concentration was determined spectrophotometrically (Nanodrop Technologies, Wilmington, USA). Comparison of the two DNA extraction methods mentioned above, using human faecal samples, indicated that both methods delivered DNA of essentially equal quality, resulting in comparable results with respect to microbial composition based on analyses with the Human Intestinal Tract Chip (HITChip), a DNA oligonucleotide microarray targeting human intestinal microbiota (Heikamp-de Jong and Hartman, personal communication). Amplification of 16S rRNA gene fragments and library preparation After DNA extraction, regions V1-V2 of the 16S rRNA genes were amplified using an in house two-step PCR protocol. In the first step, regions of interest were amplified using the following primers: 27F–DegS: GTTYGATYMTGGCTCAG (van den Bogert et al.2011) and an equimolar mix of 338R–I: GCWGCCTCCCGTAGGAGT (Daims et al.1999) and 338R–II: GCWGCCACCCGTAGGTGT (van den Bogert et al.2013), with attached UniTag I (forward) and II (reverse) linkers (I—GAGCCGTAGCCAGTCTGC; II—GCCGTGACCGTGACATCG) (Tian et al.2016). The PCR mix for one reaction at step 1 contained 10 μl of 5x HF buffer, 1 μl dNTPs (10 μM), 1U of Phusion Hot start II DNA polymerase (2U/μl), 31.5 μl of nuclease-free water, 2.5 μl of forward (10 μM) and 2.5 μl of reverse primers (10 μM), and 40 ng of DNA template. Amplification was performed in a LabCycler Gradient (SensoQuest, Germany) programmed for initial denaturation at 98°C for 30 s and 25 cycles of denaturation at 98°C for 10 s, annealing at 56°C for 20 s and extension at 72°C for 20 s, followed by final extension at 72°C for 10 min. After amplification, the success of the PCR reaction was checked visually by agarose gel electrophoresis, considering amount and size of the amplicon as quality parameters. Amplicons were subsequently used as template for a second PCR for the introduction of sample-specific barcodes, using individual barcode primers for each sample. In total, we used 48 pairs of forward and reverse barcode primers that target UniTag1 and UniTag2 sequences introduced during the first PCR, respectively, and that were appended with sample-specific barcodes. Composition of PCR reagents and cycling conditions were as described for the first PCR, with 10 μl of PCR products from the first step as template. Reactions were performed in a final volume of 100 μl. PCR products were purified and concentrated using magnetic beads (MagBio, Switzerland) according to the HighPrep protocol with adaptation for 2 ml tubes. Purified products were quantified using the Qubit dsDNA BR Assay Kit (Life Technologies, USA) following the manufacturer's protocol. PCR products were pooled in equimolar amounts into libraries of 48 samples each, and sequenced on an Illumina MiSeq platform in 300 bp paired-end mode at GATC Biotech (Constance, Germany). Data processing and statistical analysis Initial analysis of raw 16S rRNA gene sequencing data was performed using the NG-Tax pipeline (Ramiro-Garcia et al.2016). Sequences were separated into sample-specific bins based on the barcodes, after initial filtering of paired-end libraries to contain only read pairs with perfectly matching barcodes. Operational taxonomic units (OTUs) were defined using an open reference approach, and taxonomy was assigned using a SILVA 16S rRNA gene reference database (Quast et al.2013). Microbial composition plots were generated using a workflow based on Quantitative Insights Into Microbial Ecology (QIIME) v1.9.1 (Caporaso et al.2010). Reads assigned to OTUs of plant origin such as chloroplast and plant mitochondrial DNA were removed from the dataset used for downstream analyses. OTU counts were normalised using cumulative sum scaling (Paulson et al.2013). To get an overview of species composition, a normalised OTU matrix was exported to Microsoft Excel, and the relative contribution based on normalised OTU numbers per taxa was calculated. Median values of taxa relative abundance in a group of samples were used to compare groups with each other (e.g. male vs female). The OTU matrix was filtered to exclude OTUs that were present only in a small number of samples. More specifically, for each dataset, OTUs were removed that were present in less than five samples (50% of the smallest group size). Measures of alpha and beta diversity and initial multivariate analysis using principal coordinate analysis (PCoA) were performed on the rarefied matrix (depth–1650 observations) using weighted and unweighted UniFrac as distance measures as implemented in QIIME. Significance of differences in relative abundances of OTUs between individual samples was determined using Kruskal-Wallis tests when comparing more than two groups, and non-parametric t-tests with 500 Monte Carlo permutations in case of comparisons of two groups, using normalised, summarised and filtered OTU tables. False discovery rate (FDR) correction of P-values was used to reduce the chance of type I statistical errors, when multiple statistical hypotheses were tested. To identify strength and statistical significance of sample groupings with weighted and unweighted UniFrac as distance measures, we used the ‘adonis’ test as implemented in the R package ‘vegan’. Canoco 5.0 was used for multivariate statistical analysis and visualisation of correlations between microbial composition of samples and explanatory factors. Redundancy analysis (RDA) was performed as described previously (Šmilauer and Lepš 2014). As input dataset for RDA, we used the taxonomy summary table at genus level after removal of taxa which were found in less than eight samples (applied for each dataset individually) with addition of the corresponding sample metadata. No transformation or normalisation of the data was done. The significance of observed community variations was evaluated using a Monte Carlo permutation test. To identify microbial species most strongly correlated with investigated factors such as area or host species, we used the LefSe (Linear discriminant analysis of Effect Size) algorithm for biomarker identification (Segata et al.2011). Data were processed using tools developed by the Huttenhower laboratory implemented in the Galaxy environment (www.huttenhower.sph.harvard.edu/galaxy/). Preparation of input data and analysis were performed according to the standard workflow, using default settings (0.05—alpha value for the factorial Kruskal-Wallis test among classes; threshold on the logarithmic LDA score for discriminative features was 2.0; the strategy for multiclass analysis was ‘all-against-all’). After raw data processing and initial analysis, samples were organised into five sets (samples are used in several sets), allowing us to perform separate analyses and statistics while focusing on a particular research question: ‘All_samples’—initial set of all 138 samples; ‘E_fulvus’—samples obtained from E. fulvus (n = 45) from three different areas (Andasibe n = 14; Ankarafantsika n = 21; Nosy Tanikely n = 10); ‘Kirindy’—set of samples from E. rufifrons collected in Kirindy NP (n = 44); ‘Ranomafana’—samples from E. rufifrons (n = 27) and E. rubriventer (n = 22) collected in Ranomafana NP (total n = 49); ‘E_rufifrons’—samples from E. rufifrons collected from two different areas (Ranomafana NP n = 27 and Kirindy NP n = 44). Availability of data and materials Datasets generated in this study are available in the public read archive EBI, study name ‘Area of habitation strongly influences faecal microbial composition of wild lemurs’, with accession number PRJEB20007. RESULTS In this study, we analysed the faecal microbiota of a total of 138 individuals belonging to three different Eulemur species, using Illumina MiSeq sequencing of PCR-amplified 16S rRNA gene fragments covering the V1-V2 variable region. In total, we obtained 6220 515 reads, ranging from 1652 to 178 522 reads per sample (r/s) with a median of 36 092 r/s. Obtained reads were assigned to 1053 OTUs using NG-Tax, an in-house developed pipeline (Ramiro-Garcia et al. 2016). Across all samples, OTUs belonged to 12 bacterial phyla, i.e. Acidobacteria, Actinobacteria, Bacteroidetes, ‘Candidatus Saccharibacteria’ (TM7), Cyanobacteria, Firmicutes, Lentisphaerae, Proteobacteria, Spirochaetes, Synergistetes, Tenericutes and Verrucomicrobia. The fraction of non-assigned to any taxonomic level OTUs varied from 2.6% to 16.7% with an average of 9.3 ± 2.5% [s.d.] in all analysed samples. Predominant phyla, regardless of lemur species and sampling location, were Firmicutes 43.3 ± 6.4%, Bacteroidetes 30.3 ± 5.3%, Cyanobacteria 5.2 ± 3.3% and Proteobacteria 7.4 ± 3.1%. At the genus level, a total of 59 taxa were identified, 15 of which had an average relative abundance across all samples of more than 1% and comprised more than 80% of all sequences. Overall, 34% of all OTUs could be assigned at genus level and 63% at family level. Phylogenetic clustering based on relative abundance at the genus level showed that the most abundant genera could be clustered into two groups: one group consisting of the two most abundant genera (unidentified genus (UG_1 (Clostridiales) 24.9% ± 5.4%, UG_2 (Bacteroidales) 14.9 ± 3.6%) and another group consisting of another five genera (UG_3 (Prevotellaceae) 7.7 ± 2.3%, UG_4 (Cyanobacteria) 5.2 ± 3.3%, UG_5 (Bacteroidaceae) 4.6% ± 2.4%, UG_6 (Lachnospiraceae) 4.9% ± 2.8%, UG_7 (Ruminococcaceae) 4.2% ± 2%) (Fig. S1; Table S2, Supporting Information. Taxonomy of UGs). Figure 2. View largeDownload slide Differences in bacterial composition between lemur species (E. fulvus, E. rubriventer, E. rufifrons); (A) relative abundance at the phylum level of faecal microbiota of the different lemur species using dataset ‘All_samples’; (B) taxa identified by LefSe as potential biomarkers for the discrimination of studied Eulemur species; (C) taxa identified by LefSe as potential biomarkers for the discrimination of faecal samples taken in Ranomafana NP and Kirindy Forest. LDA—linear discriminant analysis. Figure 2. View largeDownload slide Differences in bacterial composition between lemur species (E. fulvus, E. rubriventer, E. rufifrons); (A) relative abundance at the phylum level of faecal microbiota of the different lemur species using dataset ‘All_samples’; (B) taxa identified by LefSe as potential biomarkers for the discrimination of studied Eulemur species; (C) taxa identified by LefSe as potential biomarkers for the discrimination of faecal samples taken in Ranomafana NP and Kirindy Forest. LDA—linear discriminant analysis. Effect of lemur species on faecal microbiota composition In order to address the influence of the host species on observed variation in microbial composition, two different datasets (see information regarding the composition of datasets in Experimental Procedures) were analysed:‘All_samples’, i.e. the entire dataset of 138 samples, and ‘Ranomafana’, the latter of which allowed us to minimise the influence of explanatory variables other than host species. Analysis of ‘All_samples’ revealed that samples from E. fulvus (n = 45) had significantly lower alpha diversity (P = 0.003) in comparison with samples taken from E. rufifrons (n = 22) (Fig. S2, Supporting Information). The relative abundance of Proteobacteria, Lentisphaerae, Synergistetes, Bacteroidetes, Actinobacteria, Firmicutes, Tenericutes and Candidate division TM7 differed among studied lemur species (FDR-corrected P < 0.012, Fig. 2A). At the genus level, 26 taxa differed in relative abundance (corrected P < 0.036), with some present only within one lemur species. The genera Anaeroplasma and UG_8 (Desulfovibrionaceae) were found only in samples belonging to E. rufifrons, albeit with relative abundances below 1%. Eight genera were identified by LefSe analysis as potential biomarkers of the different lemur species: Bacteroides and Phascolarctobacterium were identified as microbial biomarkers of E. rufifrons; UG_6 (Lachnospiraceae), Campylobacter, UG_9 (Synergistales), UG_10 (Clostridiales) and UG_11 (Xanthomonadales) were biomarkers for E. rubriventer; and UG_12 (Anaeroplasmatales) was associated with E. fulvus (Fig. 2B). No clear grouping of samples by lemur species was observed based either on weighted or unweighted UniFrac distance. This was confirmed by the ‘adonis’ test that revealed only a weak linear correlation between samples, with an R2 of 0.11 and 0.13 for weighted and unweighted distances, respectively. RDA with lemur species as the only explanatory variable showed that this variable significantly (P = 0.008) contributed to the observed variation in faecal microbiota composition (Fig. 3A). Furthermore, when comparing the faecal microbiota of E. rufifrons and E. rubriventer in Ranomafana NP, RDA showed that although ‘lemur species’ was a significant explanatory variable (P = 0.024), it only explained 6.8% of the observed variation in microbial community composition (Fig. 3B). All phyla observed in the full dataset were also present in Ranomafana NP, and only the phylum Firmicutes showed nearly significant differences in relative abundance between E. rufifrons (44.3 ± 7.1%) and E. rubriventer (39.7 ± 4.7%) (P = 0.008; corrected P = 0.1). At the genus level, UG_13 (Porphyromonadaceae), UG_5 (Bacteroidaceae) and UG_19 (Bacillales) differed in relative abundance when comparing microbial composition in both lemur species (corrected P < 0.04). Similar to dataset ‘All_samples’, no separation or grouping was observed among samples from the ‘Ranomafana’ dataset in weighted or unweighted UniFrac matrix-based PCoA plots (R2 = 0.06 and R2 = 0.05 for unweighted and weighted distance matrices, respectively; data not shown). Figure 3. View largeDownload slide Ordination triplots based on RDA with lemur species as explanatory variables. (A) In dataset ‘All_samples’, 7.2% of the variation is captured by the first two canonical axes; (B) in dataset ‘Ranomafana’, 53.7% of variation is captured by the first two canonical axes, and both host species significantly (P = 0.002) contributed to explaining the observed variation in microbiota composition. Figure 3. View largeDownload slide Ordination triplots based on RDA with lemur species as explanatory variables. (A) In dataset ‘All_samples’, 7.2% of the variation is captured by the first two canonical axes; (B) in dataset ‘Ranomafana’, 53.7% of variation is captured by the first two canonical axes, and both host species significantly (P = 0.002) contributed to explaining the observed variation in microbiota composition. Figure 4. View largeDownload slide PCoA three-dimensional (first tree PCoA axes) plots based on weighted (A, C, E) and unweighted (B, D, F) UniFrac distance matrices. Samples are represented by dots, colour-coded by sampling location. Plots contain all samples (A, B), or are species specific (E. fulvus—C, D; E. rufifrons—E, F). Figure 4. View largeDownload slide PCoA three-dimensional (first tree PCoA axes) plots based on weighted (A, C, E) and unweighted (B, D, F) UniFrac distance matrices. Samples are represented by dots, colour-coded by sampling location. Plots contain all samples (A, B), or are species specific (E. fulvus—C, D; E. rufifrons—E, F). Variation of the lemur microbiota in contrasting regions of Madagascar To determine to what extent occupancy (i.e. area of habitation) can explain the observed variation in faecal microbiota, we analysed three datasets: ‘All_samples’, ‘E_fulvus’ and ‘E_rufifrons’. Dataset ‘All_samples’ allowed us to identify the influence of ‘area’ among all other variables such as host species, sex and age. With dataset ‘E_fulvus’, we focused on a single lemur species, E. fulvus, which was sampled at three different locations, allowing us to more specifically address variation in faecal microbiota composition found within one species exposed to different environmental conditions. We constructed dataset ‘E_rufifrons’ for the same purpose as dataset ‘E_fulvus’, taking into account that individuals of E. rufifrons were sampled in two distinct areas. Dataset’ All_samples’(n = 138) showed that samples from Nosy Tanikely (n = 10) area had lower alpha diversity as compared to all other areas (P = 0.01, Fig. S3, Supporting Information). The relative abundances of 10 out of 12 phyla, except for Bacteroidetes and Acidobacteria, were different (corrected P < 0.026) between sampling sites. It should be noted that members of the phylum Acidobacteria were found in only a few samples (10 out of 138). Furthermore, at the genus level 54 out of 59 taxa that were observed in more than five samples showed significant differences in relative abundance between areas (corrected P < 0.05). Among these genera, some were found only within Kirindy NP, namely Anaeroplasma, Rhizobium and UG_8 (Desulfovibrionaceae). Members of the genus Bacteroides, UG_13 (Anaeroplasmatales), UG_17 (Clostridiales) order and UG_30 (Anaeroplasmataceae) were found exclusively in samples from the relatively dry areas (Kirindy and Ankarafantsika). The most abundant genus (UG_1 (Clostridiales)) did not vary significant between areas. Samples showed a slight visual grouping according to area in PCoA plots based on weighted UniFrac distances (R2 = 0.29), with better group separation being observed in the case of unweighted UniFrac (R2 = 0.34) (Fig. 4A and B). Furthermore, constrained analysis (RDA) showed that all areas included as explanatory variable significantly (P < 0.05) contributed to the observed variation in faecal microbiota composition (Fig. 5A). Figure 5. View largeDownload slide Ordination triplot based on RDA with areas of sampling as explanatory variables. (A) In dataset ‘All_samples’, 14.2% of the variation was captured by the first two canonical axes, and a statistically significant (P < 0.05) effect of areas as explanatory factors was observed; (B) in dataset ‘E_fulvus’, 35.7% of variation was captured by the first two canonical axes, and all three areas significantly contributed to explaining the observed variation in microbiota composition (P = 0.002); (C) in dataset ‘E_rufifrons’, 67.6% of variation was captured by the first two canonical axes. Explanatory variable ‘areas’ was statistically significant (P = 0.004) as a conditional effect. Figure 5. View largeDownload slide Ordination triplot based on RDA with areas of sampling as explanatory variables. (A) In dataset ‘All_samples’, 14.2% of the variation was captured by the first two canonical axes, and a statistically significant (P < 0.05) effect of areas as explanatory factors was observed; (B) in dataset ‘E_fulvus’, 35.7% of variation was captured by the first two canonical axes, and all three areas significantly contributed to explaining the observed variation in microbiota composition (P = 0.002); (C) in dataset ‘E_rufifrons’, 67.6% of variation was captured by the first two canonical axes. Explanatory variable ‘areas’ was statistically significant (P = 0.004) as a conditional effect. Figure 6. View largeDownload slide Heat map of relative microbial abundance at genus level in the dataset ‘E_fulvus’, with samples placed on the Y-axis, and genera with relative abundance of more than 2.5% across the dataset on the X-axis. The yellow colour indicates high relative abundance values; dark blue indicates low relative abundance values. Sample clustering and dendrogram were produced using the ‘Bray’ method as implemented in the ‘vegan’ R package. Side bar (left) indicates area of sample collection: green—Nosy Tanikely, brown—Andasibe NP, blue—Ankarafantsika. Figure 6. View largeDownload slide Heat map of relative microbial abundance at genus level in the dataset ‘E_fulvus’, with samples placed on the Y-axis, and genera with relative abundance of more than 2.5% across the dataset on the X-axis. The yellow colour indicates high relative abundance values; dark blue indicates low relative abundance values. Sample clustering and dendrogram were produced using the ‘Bray’ method as implemented in the ‘vegan’ R package. Side bar (left) indicates area of sample collection: green—Nosy Tanikely, brown—Andasibe NP, blue—Ankarafantsika. In dataset ‘E_fulvus’, 8 out of 12 phyla (Tenericutes, Cyanobacteria, Spirochaetes, Lentisphaerae, Firmicutes, Candidate division TM7, Proteobacteria, Actinobacteria) differed in relative abundance between areas (corrected P < 0.01). In line with the extensive differences observed at the phylum level, 34 out of 55 detected genera showed significant differences in relative abundance between areas. UG_17 (Clostridiales), UG_35 (Enterobacteriales), UG_29 (Spirochaetaceae), UG_12 (Anaeroplasmatales) and the genus Bacteroides were found only in the samples from Ankarafantsika NP. UG_31 (Veillonellaceae) was found exclusively in the samples from Andasibe NP, and Mesorhizobium only in Nosy Tanikely. Several genera were absent in one out of three areas: UG_19 (Bacillales), Helicobacter and Thalassospira were not detected in samples taken in Nosy Tanikely, whereas Pseudobutyrivibrio was not found in Andasibe. Samples in the ‘E_fulvus’ dataset clustered into three groups, correlating with the three different sampling sites based on the relative abundance of bacterial genera (Fig. 6). Furthermore, samples formed separated groups in PCoA plots based on weighted and unweighted UniFrac distances (R2 = 0.35 and R2 = 0.41, respectively; Fig. 4C and D). RDA showed that among all factors only area significantly (P = 0.002) contributed to explaining observed differences in faecal microbial composition of E. fulvus, with Ankarafantsika having the highest explanatory value (19.9%) (Fig. 5B). In dataset ‘E. rufifrons’, the relative abundance of the phyla Actinobacteria, Candidate division TM7 and Proteobacteria was different between the two locations where this lemur species was found, i.e. Kirindy NP and Ranomafana NP (corrected P < 0.009). In total, 26 genera were different in relative abundance when comparing both locations (corrected P < 0.047). Members of the genus Bacteroides were completely absent from the samples collected from Ranomafana NP, while their mean relative abundance was 1.4 ± 1% among samples collected from Kirindy NP. All remaining genera that were found exclusively in samples taken in Kirindy NP (UG_8 (Desulfovibrionaceae), UG_12 (Anaeroplasmatales), Anaeroplasma) had relative abundances below 0.3%. The genus Phascolarctobacterium was found in all samples from Kirindy NP (mean abundance 3 ± 1.3%), but only in 12 out of 27 samples (average abundance 0.2 ± 0.2%) from Ranomafana NP. UG_21 (Bacteroidales, 0.2 ± 0.4%), Mesorhizobium (0.1 ± 0.2%), UG_33 (Xanthobacteraceae, 0.4 ± 0.2%), Stenotrophomonas (0.1 ± 0.2%), UG_32 (Rhizobiales, 0.3 ± 0.9%) and UG_10 (Clostridiales, 0.1 ± 0.1%) were found exclusively in Ramonofana NP, albeit not in all samples collected in that area, and at low relative abundance. Additional differences included UG_5 (Bacteroidaceae) (6.2 ± 2.2% in Kirindy NP vs 3.3 ± 2% in Ranomafana NP) and UG_6 (Lachnospiraceae) (6.2 ± 2.3% in Ranomafana NP vs 2.9 ± 1.9% in Kirindy NP). Twenty one genera for Ranomofana NP and 12 for Kirindy NP were identified by LefSe as microbial ‘biomarkers’ (Fig. 2C). Multivariate analyses supported the separation of samples according to sampling location, with a clear grouping being observed in PCoA plots based on both weighted (R2 = 0.19) and unweighted (R2 = 0.23) UniFrac distance matrices (Fig. 4E and F). Furthermore, RDA showed that from all explanatory variables (area, age, sex) only variables in the group ‘area’ significantly contributed to explaining the observed variation in microbial community composition (P < 0.004), with both areas (Kirindy and Ranomafana) explaining 8.2% of variation (Fig. 5C). Variation in microbiota composition associated with sex and age of animals Influence of sex and age on faecal microbiota composition was investigated using datasets ‘All_samples’ and ‘Kirindy’. No significant differences were observed in microbiota between samples collected from males and females in any of the datasets at phylum or genus level, and no grouping was observed with multivariate analysis using PCoA and RDA (data not shown). Similarly, no significant variation in relative abundance of detected phyla was observed between age groups, and only the relative abundance of the genus Phascolarctobacterium (corrected P = 0.01) differed in age groups in dataset ‘All_samples’. Multivariate analysis (RDA and PCoA) confirmed that age did not significantly contribute to explaining the variation in faecal microbiota composition. DISCUSSION In the current study, we characterised the faecal microbiota of three frugivorous Eulemur species, and assessed to what extent the naturally occurring variation in intestinal microbiota composition is associated with occupancy, species, age and sex of individuals. Findings presented here showed that the gut microbial community of these animals is dominated by members of the phylum Firmicutes and to a lesser extent Bacteroidetes. It has been previously reported that predominance of Firmicutes or Bacteroidetes is different among animal species and mostly correlated with dietary mode and taxonomic lineage of a given species (Ley et al.2008). Our results confirmed the high proportion of Firmicutes that was previously observed in other species of frugivorous and omnivorous primates (Ochman et al.2010; Gomez et al.2015), including humans, despite the fact that phylogenetically lemurs are one of the most distinct and ancient groups within the primates (Ni et al.2013). Notably, human studies showed that the Firmicutes to Bacteroidetes ratio is not static and can be largely influenced by the presence of carbohydrates in the diet, although it is not clear which of the two phyla has a leading role as key degrader of complex carbohydrates in the human intestine. For example, on one hand an increase in relative abundance of Firmicutes was correlated with consumption of whole grains and total carbohydrate intake (Martínez et al.2013), and several species belonging to this phylum are viewed as key degraders of resistant starch (Ze et al.2012). On the other hand, it has been shown that the depletion of Firmicutes and increase in Bacteroidetes in African children from a rural area in comparison with European children were related to consumption of a traditional African diet rich in fibres and polysaccharides (De Filippo et al.2010). Such seemingly conflicting evidence might be related to the high phylogenetic and functional diversity within both phyla, including a large number of fibre- and carbohydrate-degrading species. Consequently, one can speculate that specific aspects of the diet of lemurs will result in a shift of the Firmicutes/Bacteroidetes ratio, the direction of which might not be predictable based on general characteristics of the diet. Our study showed that a large fraction of Firmicutes associated sequences was assigned to a single genus-level taxon, UG_1 (Clostridiales), accounting for 24.9 ± 5.7% of the total bacterial community. It is tempting to speculate that members of this genus have an important role in intestinal function of the three Eulemur species covered in our study; however, due to lack of physiological and ecological data, conclusions regarding their function and place in gut ecology remain speculative, awaiting isolation and/or (meta)omics analyses (Gutleben et al.2018). Notably, members of the Proteobacteria had relatively high abundance (7.4 ± 3.1%) in all three studied Eulemur species. In humans, high relative abundance of this phylum (9.7%–14.9%) has been associated with gastric bypass, metabolic disorders, inflammation and cancer, whereas its relative abundance in healthy individuals amounts to only about 4.5% (Shin, Whon and Bae 2015). On the other hand, previous research showed host species related differences in abundance of Proteobacteria among primates. For instance, Bello González et al. (2015) observed that faecal samples of humans and chimpanzees had similar relative abundances of Proteobacteria (1% and 1.2%, respectively), whereas in gorilla samples this phylum reached a relative abundance of 7% (Bello González et al.2015). Furthermore, a similar relative abundance of Proteobacteria (9.1%) was reported in faecal microbiota of L. catta (McKenney, Rodrigo and Yoder 2015). Hence, we suggest that the high relative abundance observed in this study is not necessarily a sign of a health problem of the investigated population of lemurs, but rather a feature of the normal microbial composition of frugivorous lemurs. We found that on average 5.2 ± 3.3% of all reads were assigned to OTUs belonging to the Cyanobacteria phylum. Latest research shows that members of this phylum are indeed a genuine part of the human intestinal microbiota (Di Rienzi et al.2013). Furthermore, the presence of this phylum was observed in previous studies which characterised the intestinal microbial composition of other primates, including L. catta (McKenney, Rodrigo and Yoder 2015). Furthermore, we found that 66% of genus-level taxa could not be confidently classified to a particular genus in the Silva v111 database, including several of the most predominant taxa. This is in line with the limited attention that the intestinal microbiota of lemurs has received to date, and hence there is a lack of knowledge regarding specific taxa present in the intestine of these animals. Research on intestinal microbiota of other poorly studied animals showed similar findings. Roggenbuck et al. (2014) found that only 28% of the observed genera in the giraffe rumen could be assigned to known taxa (Roggenbuck et al.2014). Similar observations have even been made for less well-characterised human populations. In a recent study, Schnorr et al. (2014) found that 22% of the total microbial community of central Tanzanian Hadza individuals could not be assigned at family and genus level, whereas this was not the case for the Italian control population (Schnorr et al.2014). To assess the role of different natural environmental factors in shaping the intestinal microbiota and how these factors relate to the influence of the different host species, we divided samples into subsets. This approach allowed us to gain a better insight into the effect of specific factors on microbiota composition, in addition to a more generic analysis of all factors at the same time in a relatively heterogeneous dataset. The value of this tiered approach was confirmed by the fact that for all factors of interest, a first insight into potential effects that could be obtained with the whole dataset received additional, more robust, support from the analysis of specific subsets of samples. It should also be noted that, due to the nature of wildlife sampling under natural conditions, it remains challenging to obtain balanced sample sets with equal numbers of samples in each group. We discovered that in samples analysed in our study, the most influential factor contributing to shaping microbiota composition was the area of sampling. When we applied PCoA based on either weighted or unweighted UniFrac distances, separation into areas was obvious in all datasets when included. Remarkably, clustering of samples was tighter with better separation of samples when unweighted UniFrac was used as a distance measure. This observation suggests, taking into account the nature of the UniFrac distance calculation, that the faecal microbiota of lemurs from different areas is more distinct with respect to microbial species composition than in relative abundance of prevalent taxa. Constrained multivariate analysis (RDA) confirmed that occupancy is the most influential explanatory variable with respect to the observed variation in lemur intestinal microbiota composition. Madagascar is known to have different environmental conditions and biodiversity within relatively small areas (Wilmé, Goodman and Ganzhorn 2006). Furthermore, sampling locations were positioned with considerable distance from each other, and were characterised by major climatic differences such as amount of precipitation and forest composition. Hence, it is likely that the availability of food items during the year, in particular fruits and flowers, is the driving force that leads to differences in the intestinal microbiota. These food items are scarce during the dry season in the dry deciduous forest areas such as Kirindy NP (Sorg and Rohner 1996; Ganzhorn 2002), whereas in the areas characterised by wet rainforests such as Ranomafana NP and Andasibe NP, as well as in Nosy Tanikely, these food items are abundantly available almost all year around. Surprisingly, relatively low microbial alpha diversity was observed in the tropical rainforest. It should be noted, however, that samples from this area only included E. fulvus. It is tempting to speculate that the lower alpha diversity in E. fulvus might be explained by adaptation of the microbiota to sugar-rich food that is prevalent in Nosy Tanikely due to the abundant presence of mango trees. We also observed differences in microbiota composition related to the species of lemurs; however, these differences were secondary to those observed between different areas of habitation. Host genetic differences are among the major driving forces shaping intestinal microbiota composition (Khachatryan et al.2008), including different animal species with similar dietary habits (Moeller et al.2014). In our study, this was not the case. One explanation for this finding could be that this study has been conducted on congeneric species which by definition are genetically close (Markolf and Kappeler 2013), and have almost identical digestive systems (Tattersall 1993), leading to the moderate effect of genetics on gut microbiota composition observed here. We did not find any evidence for an influence of sex or age on microbiota composition. It should be noted, however, that in this study different age classes were not equally represented in the datasets, as we mostly sampled adult individuals (74.6% of all samples). Furthermore, all of the non-adults were at or beyond juvenile stage. Based on knowledge about microbiota development of human infants, the transformation of microbiota to an adult-like mature, composition occurs before reaching the juvenile stage (Koenig et al.2011; Wopereis et al.2014). Surprisingly, we did not find any differences in beta- and alpha diversities of microbiota between males and females. Many studies showed an influence of this factor in different species (Bolnick et al.2014), including lemurs (Aivelo, Laakkonen and Jernvall 2016). However, it was pointed out before that other factors, such as host genetics, can outweigh the influence of this factor (Kovacs et al.2011). Hence, it is tempting to speculate that in the present study the effect of sex on the faecal microbiota of the different Eulemur species might have been obscured by more influential factors as well as relatively large variation in microbial composition between individual animals. In conclusion, we showed that intestinal microbiota in three genetically close species of lemurs was most strongly influenced by their occupancy, whereas the influence of genetic differences was minor, and influence of sex and age was not detectable. All three lemur species had similar bacterial composition in terms of predominant and prevalent bacterial taxa. The findings reported here contribute to our knowledge about the intestinal microbiota in non-human primates, and factors that shape the bacterial composition in wild lemur populations, which can be extrapolated into general rules of intestinal microbiota assembly. Furthermore, the high fraction of poorly assigned taxa reinforces the notion that microbiota of non-humanoid primates has so far received little attention, harbouring a broad range of potentially novel bacterial species and genera that deserve attention in future studies. SUPPLEMENTARY DATA Supplementary data are available at FEMSEC online. Acknowledgements For their advice in the field, the authors want to thank Dr P. Wright, Dr E. Larney and John E. Cadle, and for their logistical support we want to thank the Centre Valbio staff (CVB), the Institute for the Conservation of Tropical Environments (ICTE) and the Madagascar Institute pour la Conservation des Ecosystèmes Tropicaux (MICET), especially B. Andriamihaja, who facilitated arrangements on necessary permits. Thanks go to the Malagasy government and Madagascar National Parks (MNP), Centre de Formation Professionelle Forestière (CFPF), and Mitsinjo Association for their formal permission. For their fieldwork assistance and indispensable knowledge, we want to thank all field technicians and guides that helped as at the various national parks. FUNDING This work was supported by an ‘Nederlandse Organisatie voor Wetenschappelijk Onderzoek’ and ‘The Graduate School for Production Ecology and Resource Conservation’ grant [Grant nr. 1208; project nr. 5120873001], and the European Union Erasmus-Mundus PhD fellowship [EM Action 2- Partnership IAMONET- 20132520]. Other necessary funds to conduct the fieldwork of this study were provided by the Koninklijke Nederlandse Akademie van Wetenschappen (KNAW) Academy Ecology Fund, the Dr J. L. Dobberke Foundation, the Treub Foundation, the Dutch Fund for Research on Nature Conservation, the Dutch Royal Botanical Society and the Dutch Royal Zoological Society. Conflict of interest. None declared. REFERENCES Aivelo T, Laakkonen J, Jernvall J. Population- and individual-level dynamics of the intestinal microbiota of a small primate. Appl Environ Microb  2016; 82: 3537– 45. Google Scholar CrossRef Search ADS   Amato KR, Leigh SR, Kent A et al.   The gut microbiota appears to compensate for seasonal diet variation in the wild black howler monkey (Alouatta pigra). Microb Ecol  2015; 69: 434– 43. Google Scholar CrossRef Search ADS PubMed  Atsalis S. Spatial distribution and population composition of the brown mouse lemur (Microcebus rufus) in Ranomafana National Park, Madagascar, and its implications for social organization. Am J Primatol  2000; 51: 61– 78. Google Scholar CrossRef Search ADS PubMed  Balko EA, Brian Underwood H. Effects of forest structure and composition on food availability forVarecia variegata at Ranomafana National Park, Madagascar. Am J Primatol  2005; 66: 45– 70. Google Scholar CrossRef Search ADS PubMed  Bello González T, van Passel M, Tims S et al.   Application of the Human Intestinal Tract Chip to the non-human primate gut microbiota. Benef Microbes  2015; 6: 271– 6. Google Scholar CrossRef Search ADS PubMed  Bolnick DI, Snowberg LK, Hirsch PE et al.   Individual diet has sex-dependent effects on vertebrate gut microbiota. Nat Commun  2014; 5: 4500. Google Scholar CrossRef Search ADS PubMed  Caporaso JG, Kuczynski J, Stombaugh J et al.   QIIME allows analysis of high-throughput community sequencing data. Nat Methods  2010; 7: 335– 6. Google Scholar CrossRef Search ADS PubMed  Clemente JC, Ursell LK, Parfrey LW et al.   The impact of the gut microbiota on human health: an integrative view. Cell  2012; 148: 1258– 70. Google Scholar CrossRef Search ADS PubMed  Daims H, Brühl A, Amann R et al.   The domain-specific probe EUB338 is insufficient for the detection of all Bacteria: development and evaluation of a more comprehensive probe set. Syst Appl Microbiol  1999; 22: 434– 44. Google Scholar CrossRef Search ADS PubMed  David LA, Materna AC, Friedman J et al.   Host lifestyle affects human microbiota on daily timescales. Genome Biol  2014; 15: R89. Google Scholar CrossRef Search ADS PubMed  De Filippo C, Cavalieri D, Di Paola M et al.   Impact of diet in shaping gut microbiota revealed by a comparative study in children from Europe and rural Africa. P Natl Acad Sci USA  2010; 107: 14691– 6. Google Scholar CrossRef Search ADS   Di Rienzi SC, Sharon I, Wrighton KC et al.   The human gut and groundwater harbor non-photosynthetic bacteria belonging to a new candidate phylum sibling to Cyanobacteria. Elife  2013; 2: e01102. Google Scholar CrossRef Search ADS PubMed  Ellis RJ, Bruce KD, Jenkins C et al.   Comparison of the distal gut microbiota from people and animals in Africa. PLoS One  2013; 8: e54783. Google Scholar CrossRef Search ADS PubMed  Fogel AT. The gut microbiome of wild lemurs: a comparison of sympatric Lemur catta and Propithecus verreauxi. Folia Primatol  2015; 86: 85– 95. Google Scholar CrossRef Search ADS PubMed  Ganzhorn JU. Distribution of a folivorous lemur in relation to seasonally varying food resources: integrating quantitative and qualitative aspects of food characteristics. Oecologia  2002; 131: 427– 35. Google Scholar CrossRef Search ADS PubMed  Ganzhorn JU, Schmid J. Different population dynamics of Microcebus murinus in primary and secondary deciduous dry forests of Madagascar. Int J Primatol  1998; 19: 785– 96. Google Scholar CrossRef Search ADS   García G, Goodman SM. Hunting of protected animals in the Parc National d’Ankarafantsika, north-western Madagascar. Oryx  2003; 37: 115– 8. Google Scholar CrossRef Search ADS   Gomez A, Petrzelkova K, Yeoman CJ et al.   Gut microbiome composition and metabolomic profiles of wild western lowland gorillas (Gorilla gorilla gorilla) reflect host ecology. Mol Ecol  2015; 24: 2551– 65. Google Scholar CrossRef Search ADS PubMed  Goodman SM, Benstead JP. Updated estimates of biotic diversity and endemism for Madagascar. Oryx  2005; 39: 73– 77. Google Scholar CrossRef Search ADS   Gutleben J, Chaib De Mares M, van Elsas JD et al.   The multi-omics promise in context: from sequence to microbial isolate. Crit Rev Microbiol  2018; 44: 212– 29. Google Scholar CrossRef Search ADS PubMed  Hansen J, Gulati A, Sartor RB. The role of mucosal immunity and host genetics in defining intestinal commensal bacteria. Curr Opin Gastroenterol  2010; 26: 564– 71. Google Scholar CrossRef Search ADS PubMed  Irwin MT, Wright PC, Birkinshaw C et al.   Patterns of species change in anthropogenically disturbed forests of Madagascar. Biol Conserv  2010; 143: 2351– 62. Google Scholar CrossRef Search ADS   Johnson SE, Gordon AD, Stumpf RM et al.   Morphological variation in populations of Eulemur albocollaris and E. fulvus rufus. Int J Primatol  2005; 26: 1399– 416. Google Scholar CrossRef Search ADS   Kabat AM, Srinivasan N, Maloy KJ. Modulation of immune development and function by intestinal microbiota. Trends Immunol  2014; 35: 507– 17. Google Scholar CrossRef Search ADS PubMed  Khachatryan ZA, Ktsoyan ZA, Manukyan GP et al.   Predominant role of host genetics in controlling the composition of gut microbiota. PLoS One  2008; 3: e3064. Google Scholar CrossRef Search ADS PubMed  Koenig JE, Spor A, Scalfone N et al.   Succession of microbial consortia in the developing infant gut microbiome. P Natl Acad Sci USA  2011; 108: 4578– 85. Google Scholar CrossRef Search ADS   Kovacs A, Ben-Jacob N, Tayem H et al.   Genotype is a stronger determinant than sex of the mouse gut microbiota. Microb Ecol  2011; 61: 423– 8. Google Scholar CrossRef Search ADS PubMed  Lewis RJ, Bannar-Martin KH. The impact of cyclone Fanele on a tropical dry forest in Madagascar. Biotropica  2012; 44: 135– 40. Google Scholar CrossRef Search ADS   Ley RE, Hamady M, Lozupone C et al.   Evolution of mammals and their gut microbes. Science  2008; 320: 1647– 51. Google Scholar CrossRef Search ADS PubMed  Markolf M, Kappeler PM. Phylogeographic analysis of the true lemurs (genus Eulemur) underlines the role of river catchments for the evolution of micro-endemism in Madagascar. Front Zool  2013; 10: 70. Google Scholar CrossRef Search ADS PubMed  Martínez I, Lattimer JM, Hubach KL et al.   Gut microbiome composition is linked to whole grain-induced immunological improvements. ISME J  2013; 7: 269– 80. Google Scholar CrossRef Search ADS PubMed  Maukonen J, Saarela M. Human gut microbiota: does diet matter? Proc Nutr Soc  2015; 74: 23– 36. Google Scholar CrossRef Search ADS PubMed  McKenney EA, Rodrigo A, Yoder AD. Patterns of gut bacterial colonization in three primate species. PLoS One  2015; 10: e0124618. Moeller AH, Li Y, Ngole EM et al.   Rapid changes in the gut microbiome during human evolution. P Natl Acad Sci USA  2014; 111: 16431– 5. Google Scholar CrossRef Search ADS   Moeller AH, Caro-Quintero A, Mjungu D et al.   Cospeciation of gut microbiota with hominids. Science  2016; 353: 380– 2. Google Scholar CrossRef Search ADS PubMed  Ni X, Gebo DL, Dagosto M et al.   The oldest known primate skeleton and early haplorhine evolution. Nature  2013; 498: 60– 4. Google Scholar CrossRef Search ADS PubMed  Ochman H, Worobey M, Kuo C-H et al.   Evolutionary relationships of wild hominids recapitulated by gut microbial communities. PLoS Biol  2010; 8: 2638. Google Scholar CrossRef Search ADS   Overdorff DJ. Similarities, differences, and seasonal patterns in the diets of Eulemur rubriventer and Eulemur fulvus rufus in the Ranomafana National Park, Madagascar. Int J Primatol  1993; 14: 721– 53. Google Scholar CrossRef Search ADS   Overdorff DJ. Ecological correlates to activity and habitat use of two prosimian primates: Eulemur rubriventer and Eulemur fulvus rufus in Madagascar. Am J Primatol  1996; 40: 327– 42. Google Scholar CrossRef Search ADS   Overdorff DJ, Strait SG, Telo A. Seasonal variation in activity and diet in a small-bodied folivorous primate, Hapalemur griseus, in southeastern Madagascar. Am J Primatol  1997; 43: 211– 23. Google Scholar CrossRef Search ADS PubMed  Patterson E, Cryan J, Fitzgerald G et al.   Gut microbiota, the pharmabiotics they produce and host health. Proc Nutr Soc  2014; 73: 477– 89. Google Scholar CrossRef Search ADS PubMed  Paulson JN, Stine OC, Bravo HC et al.   Differential abundance analysis for microbial marker-gene surveys. Nat Methods  2013; 10: 1200– 2. Google Scholar CrossRef Search ADS PubMed  Pyritz LW, Kappeler PM, Fichtel C. Coordination of group movements in wild red-fronted lemurs (Eulemur rufifrons): processes and influence of ecological and reproductive seasonality. Int J Primatol  2011; 32: 1325– 47. Google Scholar CrossRef Search ADS PubMed  Quast C, Pruesse E, Yilmaz P et al.   The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res  2013; 41: D590– 6. Google Scholar CrossRef Search ADS PubMed  Rakotonirina G. Composition and structure of a dry forest on sandy soils near Morondava. Primate Rep  1996: 81– 87. Ramiro-Garcia, J, Hermes, GDA, Giatsis C et al.   NG-Tax, a highly accurate and validated pipeline for analysis of 16S rRNA amplicons from complex biomes. F1000Res  2016; 5. Ren T, Grieneisen LE, Alberts SC et al.   Development, diet and dynamism: longitudinal and cross-sectional predictors of gut microbial communities in wild baboons. Environ Microbiol  2016; 18: 1312– 25. Google Scholar CrossRef Search ADS PubMed  Roggenbuck M, Sauer C, Poulsen M et al.   The giraffe (Giraffa camelopardalis) rumen microbiome. FEMS Microbiol Ecol  2014; 90: 237– 46. Google Scholar CrossRef Search ADS PubMed  Salonen A, Nikkilä J, Jalanka-Tuovinen J et al.   Comparative analysis of fecal DNA extraction methods with phylogenetic microarray: effective recovery of bacterial and archaeal DNA using mechanical cell lysis. J Microbiol Methods  2010; 81: 127– 34. Google Scholar CrossRef Search ADS PubMed  Sato H. Seasonal fruiting and seed dispersal by the brown lemur in a tropical dry forest, north-western Madagascar. J Trop Ecol  2013; 29: 61– 69. Google Scholar CrossRef Search ADS   Schnorr SL, Candela M, Rampelli S et al.   Gut microbiome of the Hadza hunter-gatherers. Nat Commun  2014; 5: 3654. Google Scholar CrossRef Search ADS PubMed  Segata N, Izard J, Waldron L et al.   Metagenomic biomarker discovery and explanation. Genome Biol  2011; 12: R60. Google Scholar CrossRef Search ADS PubMed  Shin N-R, Whon TW, Bae J-W. Proteobacteria: microbial signature of dysbiosis in gut microbiota. Trends Biotechnol  2015; 33: 496– 503. Google Scholar CrossRef Search ADS PubMed  Šmilauer P, Lepš J. Multivariate Analysis of Ecological Data Using CANOCO 5 . New York: Cambridge University Press, 2014. Google Scholar CrossRef Search ADS   Sorg J, Ganzhorn J, Kappeler P. Forestry and research in the Kirindy Forest/Centre de Formation Professionnelle Forestière. In: Goodman SM, Benstead JP. (eds). The Natural History of Madagascar . Chicago: The University of Chicago Press, 2003. Sorg J, Rohner U. Climate and tree phenology of the dry deciduous forest of the Kirindy Forest. Primate Rep  1996; 46: 57– 80. Sussman R. Ecological distinctions in sympatric species of Lemur. In: Martin RD. (ed). Prosimian Biology . London: Gerald Duckworth & Co Ltd, 1974. Tattersall I. Speciation and morphological differentiation in the genus Lemur. In: Kimbel WH, Martin LB. (ed). Species, Species Concepts and Primate Evolution . New York: Springer, 1993, 163– 76. Google Scholar CrossRef Search ADS   Tecot SR. Seasonality and predictability: the hormonal and behavioral responses of the red-bellied lemur, Eulemur rubriventer in Southeastern Madagascar, 2008, https://repositories.lib.utexas.edu/bitstream/handle/2152/3947/tecots23568.pdf?sequence=2&isAllowed=y. Tecot SR. It's all in the timing: birth seasonality and infant survival in Eulemur rubriventer. Int J Primatol  2010; 31: 715– 35. Google Scholar CrossRef Search ADS   Tian L, Scholte J, Borewicz K et al.   Effects of pectin supplementation on the fermentation patterns of different structural carbohydrates in rats. Mol Nutr Food Res  2016; 60: 2256– 66. Google Scholar CrossRef Search ADS PubMed  van den Bogert B, de Vos WM, Zoetendal EG et al.   Microarray analysis and barcoded pyrosequencing provide consistent microbial profiles depending on the source of human intestinal samples. Appl Environ Microb  2011; 77: 2071– 80. Google Scholar CrossRef Search ADS   van den Bogert B, Erkus O, Boekhorst J et al.   Diversity of human small intestinal Streptococcus and Veillonella populations. FEMS Microbiol Ecol  2013; 85: 376– 88. Google Scholar CrossRef Search ADS PubMed  Vences M, Köhler J, Glaw F. First record of Mabuya comorensis (Reptilia: Scincidae) for the Madagascan fauna, with notes on the reptile fauna of the offshore island Nosy Tanikely. Museo Regionale di Scienze Naturali . 1997; 15–N1. Wilmé L, Goodman SM, Ganzhorn JU. Biogeographic evolution of Madagascar's microendemic biota. Science  2006; 312: 1063– 5. Google Scholar CrossRef Search ADS PubMed  Wopereis H, Oozeer R, Knipping K et al.   The first thousand days - intestinal microbiology of early life: establishing a symbiosis. Pediatr Allergy Immunol  2014; 25: 428– 38. Google Scholar CrossRef Search ADS PubMed  Wright PC, Erhart EM, Tecot S et al.   Long-term lemur research at Centre ValBio, Ranomafana National Park, Madagascar. In: Kappler PM, Watts DP. (eds). Long-Term Field Studies of Primates . New York: Springer, 2012. Google Scholar CrossRef Search ADS   Yildirim S, Yeoman CJ, Sipos M et al.   Characterization of the fecal microbiome from non-human wild primates reveals species specific microbial communities. PLoS One  2010; 5: e13963. Google Scholar CrossRef Search ADS PubMed  Yu Z, Morrison M. Improved extraction of PCR-quality community DNA from digesta and fecal samples. BioTechniques  2004; 36: 808– 13. Google Scholar PubMed  Ze X, Duncan SH, Louis P et al.   Ruminococcus bromii is a keystone species for the degradation of resistant starch in the human colon. ISME J  2012; 6: 1535– 43. Google Scholar CrossRef Search ADS PubMed  © FEMS 2018. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png FEMS Microbiology Ecology Oxford University Press

Loading next page...
 
/lp/ou_press/occupancy-strongly-influences-faecal-microbial-composition-of-wild-B0SWmmkDsZ
Publisher
Blackwell
Copyright
© FEMS 2018. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com
ISSN
0168-6496
eISSN
1574-6941
D.O.I.
10.1093/femsec/fiy017
Publisher site
See Article on Publisher Site

Abstract

Abstract The microbiota of the mammalian gut is a complex ecosystem, the composition of which is greatly influenced by host genetics and environmental factors. In this study, we aim to investigate the influence of occupancy (a geographical area of habitation), species, age and sex on intestinal microbiota composition of the three lemur species: Eulemur fulvus, E. rubriventer and E. rufifrons. Faecal samples were collected from a total of 138 wild lemurs across Madagascar, and microbial composition was determined using next-generation sequencing of PCR-amplified 16S rRNA gene fragments. Consistent with reports from other primate species, the predominant phyla were Firmicutes (43 ± 6.4% [s.d.]) and Bacteroidetes (30.3 ± 5.3%). The microbial composition was strongly associated with occupancy in the E. fulvus population, with up to 19.9% of the total variation in microbial composition being explained by this factor. In turn, geographical differences observed in faecal microbiota of sympatric lemur species were less pronounced, as was the impact of the factors sex and age. Our findings showed that among the studied factors occupancy had the strongest influence on intestinal microbiota of congeneric lemur species. This suggests adaptation of microbiota to differences in forest composition, climate variations and correspondingly available diet in different geographical locations of Madagascar. microbiota, Eulemur, Madagascar, gastro-intestinal tract, environment, multivariate statistics INTRODUCTION The intestinal microbiota of mammals is an integral part of an animal's body that contributes significantly to the overall health of the host through modulation of its immune system, facilitation of food digestion, competition with pathogenic microorganisms and production of metabolites beneficial for the host (Kabat, Srinivasan and Maloy 2014; Patterson et al.2014). Expression of these beneficial properties directly correlates with microbial community diversity and composition (Clemente et al.2012). Hence, identifying the factors and underlying processes that shape the intestinal microbiota is important for a better understanding of its contribution to host health. Previous studies in humans have shown that host genetics (Hansen, Gulati and Sartor 2010), lifestyle (David et al.2014) and food preferences (Maukonen and Saarela 2015) contribute to shaping microbiota composition of an individual within a population. Intestinal microbiota composition can be distinguished between different mammalian species, suggesting co-evolution and adaptation of animals and their microbes (Ley et al.2008; Moeller et al.2016). It is not clear, however, to what extent host genotype and environmental factors influence intestinal microbial composition under natural conditions among closely related animal species dwelling in different biogeographical regions. Although wildlife microbiota has received less attention in comparison with that of humans, farm and rodent model animals, data collected from animals in wild conditions can provide complementary information that contributes to our understanding of processes that shape mammalian intestinal microbiota. For instance, studies that highlight similarities and differences in microbiota between humans and other Homininae species (Ellis et al.2013; Moeller et al.2014; Schnorr et al.2014) provided new insight into evolution of microbiota, suggesting adaptation of human microbes to an animal protein-based diet. Studies on microbiota composition of primates that are evolutionarily more distant from humans, such as yellow baboons (Papio cynocephalus) (Ren et al. 2016), black howler monkey (Alouatta pigra) (Amato et al.2015), black and white colobus (Colubus guereza), red colobus (Piliocolobus tephrosceles) and red-tailed guenon (Cercopithecus ascanius) (Yildirim et al.2010) revealed that microbiota composition of these primates is highly variable, also intra-individually, and mostly depends on the available diet. Correspondingly, the diet of a wild animal directly depends on suitable food availability, which depends on climate, flora and fauna of an area. This statement is also true for wild lemurs of Madagascar, for which several studies showed variation in feeding patterns and diet when comparing areas with different forest composition (Balko and Brian Underwood 2005), as well as different seasons (Overdorff 1993; Overdorff, Strait and Telo 1997). Fogel (2015) compared the microbiota composition of sympatric wild Lemur catta and Propithecus verreauxi across dry and wet seasons and showed that microbiota of both lemur species is variable between individuals and dynamic over time. Researchers observed differences in microbial composition between wild vs captive L. catta as well as wild populations of L. catta and P. verreauxi, albeit only with respect to relative abundance of specific microbial groups rather than their presence or absence (McKenney, Rodrigo and Yoder 2015). Wild rufous mouse lemurs (Microcebus rufus) showed an increase in gut microbial diversity with age and differences in microbiota richness and diversity between sampling. Furthermore, microbial composition was affected by site, sex and year, whereas temporal trends within a year were weak (Aivelo, Laakkonen and Jernvall 2016). The above-mentioned studies of lemur microbiota were focused on a single lemur species (Aivelo, Laakkonen and Jernvall 2016), two sympatric lemur species dwelling in the same area (Fogel 2015) or captive lemurs of different species (McKenney, Rodrigo and Yoder 2015). Taken together, these studies showed that lemurs harbour complex intestinal microbiota, the composition of which fluctuates over time among and within individuals, and is affected by season, captivity, age, site of sampling and sex. In our study we focused on microbiota of three closely related Eulemur species (E. fulvus, E. rufifrons and E. rubriventer). To the best of our knowledge, this is the first comparative study of intestinal microbiota composition of multiple wild Eulemur species across Madagascar exposed to large variations in climate conditions and biogeography. In addition to an explorative assessment of the most important features of intestinal microbial composition in the studied species, we addressed to what extent occupancy, host species, sex and age influence lemur intestinal microbial composition, and which of these factors contribute most strongly to intestinal microbiota differentiation in wild lemurs. To this end, we hypothesised that intestinal microbial composition is similar among congeneric lemur species and that of the remaining factors occupancy, with habitat as a determining factor of food items availability and other climate-related effects, has the strongest influence on intestinal microbiota differentiation. MATERIAL AND METHODS Study design Faecal samples (N = 138) were selected from wild lemurs collected across Madagascar from April to July 2014 (Fig. 1). To investigate the effect of the occupancy on intestinal microbiota, E. fulvus samples from three climatic regions and E. rufifrons samples from two climatic regions were compared with each other. To assess the influence of different species, E. rubriventer samples collected in Ranomafana National Park (NP) were compared to E. rufifrons samples from the same area. The effect of age and sex was estimated based on E. rufifrons samples collected in Kirindy NP and Ranomafana NP (Table S1, Supporting Information). Figure 1. View largeDownload slide Lemur faecal sample collection areas across Madagascar. The map showing main types of land cover and vegetation was adapted from www.wildmadagascar.org, and was produced with data taken from the FAO Country Profiles and Mapping Information System (The United Nations Food and Agricultural Organization; © FAO 2004). Faecal samples were collected at five geographical locations across the island from E. rufifrons (two sites), E. fulvus (three sites) and E. rubriventer (one site). n = number of samples. Figure 1. View largeDownload slide Lemur faecal sample collection areas across Madagascar. The map showing main types of land cover and vegetation was adapted from www.wildmadagascar.org, and was produced with data taken from the FAO Country Profiles and Mapping Information System (The United Nations Food and Agricultural Organization; © FAO 2004). Faecal samples were collected at five geographical locations across the island from E. rufifrons (two sites), E. fulvus (three sites) and E. rubriventer (one site). n = number of samples. Study sites Madagascar experiences a strong variation in climate conditions, resulting in different vegetation zones across the island (Irwin et al.2010). The effect of environmental factors on lemur microbiota composition was investigated at five sites across Madagascar (Fig. 1). Kirindy Mitea NP (20°07΄S, 44°67΄E, 722 km2) and Ankarafantsika NP (16°25΄S, 46°80΄E, 1350 km2) consist of dry deciduous forest and are located on the western- and north-western side of Madagascar, respectively (Goodman and Benstead 2005). Kirindy Mitea NP is characterised by pronounced seasonality. The area contains more than 200 species of trees with a mean canopy height of 12–18 m, containing mostly deciduous trees with adaptations to water stress (Lewis and Bannar-Martin 2012). Ankarafantsika NP is a mosaic of floristically heterogeneous dry deciduous forests dissected by small valleys with abundant Raffia palms (Rakotonirina 1996; Ganzhorn and Schmid 1998; García and Goodman 2003; Sorg, Ganzhorn and Kappeler 2003) Ranomafana NP is located in southeastern Madagascar (21°16΄S, 47°20΄E) and encompasses ∼435 km2 of montane moist forest, ranging from altitudes of 500 m up to 1500 m, and receives an average of 3000 mm rainfall per year (Wright et al.2012). The rainfall in Ranomafana NP differs highly between the wet-warm season (December to March, 482–1170 mm per month) and dry-cold season (April to November, 55–513 mm per month) (Atsalis 2000). Andasibe Mantadia NP (155 km2) is located at the eastern side of Madagascar (18°92΄S, 48°42΄E), and is also characterised by relatively wet rain forests. Nosy Tanikely (13°28΄S, 48°14΄E) is an island in the north-east of Madagascar covered with tropical vegetation. This island comprises less than 0.3 km2 and is located between Nosy Be (8 km) and the mainland of Madagascars (13 km). Elevation ranges from 0 to 47 m above sea level (Vences, Köhler and Glaw 1997). The island's vegetation consists of low forest with planted banana and mango trees, surrounded by a sandy shore with large rock formations (de Winter, personal observation) Studied species This study focused on three Eulemur species: the red-fronted lemur (E. rufifrons), the common brown lemur (E. fulvus) and the red-bellied lemur (E. rubriventer). These species are morphologically alike and are frugivorous, although they may include other food sources, such as leaves and invertebrates, in their diet (Sussman 1974; Overdorff 1996; Pyritz, Kappeler and Fichtel 2011; Sato 2013). The main difference in social organisation between the different Eulemur species is their group size. Eulemur rufifrons and E. fulvus live in multimale, multifemale groups from 4 to 18 individuals (Overdorff 1996; Pyritz, Kappeler and Fichtel 2011; Johnson et al.2005), whereas E. rubriventer lives in small monogamous groups from two up to five individuals (Tecot 2008; Tecot 2010). Sampling and data collection Immediately after defecation, fresh faecal samples (3–4 g) from individual lemurs were collected non-invasively. Within 12 h after collection, the samples were stored at ambient temperature in sterile plastic tubes that were prefilled with 5 ml of 70% ethanol until further analyses at the Laboratory of Microbiology, Wageningen University, The Netherlands. Species, age and sex were recorded. All samples included in this paper were taken in compliance with the laws of the Government of Madagascar, and no animal experimentation was involved. DNA extraction Samples collected in the Ranomafana NP were processed using a modified protocol based on method proposed by Yu and Morrison (2004) with modifications described previously (Salonen et al.2010). For this method, faecal material was air-dried for 15–20 min in a fume hood to remove ethanol from the samples. Subsequently, 0.1–0.17 g of dried samples were added into double autoclaved screw-cap tubes containing 0.3 g of 0.1 mm zirconia beads, three pieces of 2.5 mm glass beads and 700 μl of lysis buffer (500 nM NaCl, 50 mM Tris-HCI (pH = 8), 50 mM EDTA, 4% SDS) in each. Samples were treated for 3 × 1 min at 5.5 × 103 movements per minute in a Precellys 24 beadbeater (Bertin technologies, France). After homogenisation, samples were incubated at 95°C for 15 min in a shaking heating block (Vartemp 56, Labnet International, Edison, NJ, USA) at 100 rpm, then centrifuged at 4°C for 5 min at 13 000 rpm. Clean supernatants were transferred into 2 ml tubes. Three hundred microlitres of fresh lysis buffer was added in the same tubes to the pellets, bead beating/incubation steps were repeated and freshly collected supernatant was pooled with that previously collected. Subsequent steps were performed according to the original protocol (Yu and Morrison 2004). Samples collected in Adnasibe NP, Kirindy Mitea NP, Ankarafantsika NP and Nosy Tanikely were extracted using an automatic system, the Maxwell® 16 Research Instrument (Promega, Madison, USA) and the corresponding RNA extraction kit according to the manufacturer's instructions. To improve DNA yield, samples preserved in 70% ethanol were rehydrated through a series of ethanol solutions with decreasing proportions of ethanol in steps of 10%. For rehydration, 1.5 ml of 70% ethanol with faecal particles was transferred into a fresh 2 ml tube and centrifuged at 13 000 rpm for 5 min to separate solid fractions from the liquid. After centrifugation, part of the supernatant was replaced with the same amount of distilled water to decrease ethanol concentration by 10% point, vortexed and incubated for 10 min at RT. These steps were repeated until the ethanol was replaced by distilled water. Cell disruption and lysis was performed as described above, but instead of lysis buffer we used S.T.A.R buffer (Roche Molecular Systems, USA). DNA quality and concentration was determined spectrophotometrically (Nanodrop Technologies, Wilmington, USA). Comparison of the two DNA extraction methods mentioned above, using human faecal samples, indicated that both methods delivered DNA of essentially equal quality, resulting in comparable results with respect to microbial composition based on analyses with the Human Intestinal Tract Chip (HITChip), a DNA oligonucleotide microarray targeting human intestinal microbiota (Heikamp-de Jong and Hartman, personal communication). Amplification of 16S rRNA gene fragments and library preparation After DNA extraction, regions V1-V2 of the 16S rRNA genes were amplified using an in house two-step PCR protocol. In the first step, regions of interest were amplified using the following primers: 27F–DegS: GTTYGATYMTGGCTCAG (van den Bogert et al.2011) and an equimolar mix of 338R–I: GCWGCCTCCCGTAGGAGT (Daims et al.1999) and 338R–II: GCWGCCACCCGTAGGTGT (van den Bogert et al.2013), with attached UniTag I (forward) and II (reverse) linkers (I—GAGCCGTAGCCAGTCTGC; II—GCCGTGACCGTGACATCG) (Tian et al.2016). The PCR mix for one reaction at step 1 contained 10 μl of 5x HF buffer, 1 μl dNTPs (10 μM), 1U of Phusion Hot start II DNA polymerase (2U/μl), 31.5 μl of nuclease-free water, 2.5 μl of forward (10 μM) and 2.5 μl of reverse primers (10 μM), and 40 ng of DNA template. Amplification was performed in a LabCycler Gradient (SensoQuest, Germany) programmed for initial denaturation at 98°C for 30 s and 25 cycles of denaturation at 98°C for 10 s, annealing at 56°C for 20 s and extension at 72°C for 20 s, followed by final extension at 72°C for 10 min. After amplification, the success of the PCR reaction was checked visually by agarose gel electrophoresis, considering amount and size of the amplicon as quality parameters. Amplicons were subsequently used as template for a second PCR for the introduction of sample-specific barcodes, using individual barcode primers for each sample. In total, we used 48 pairs of forward and reverse barcode primers that target UniTag1 and UniTag2 sequences introduced during the first PCR, respectively, and that were appended with sample-specific barcodes. Composition of PCR reagents and cycling conditions were as described for the first PCR, with 10 μl of PCR products from the first step as template. Reactions were performed in a final volume of 100 μl. PCR products were purified and concentrated using magnetic beads (MagBio, Switzerland) according to the HighPrep protocol with adaptation for 2 ml tubes. Purified products were quantified using the Qubit dsDNA BR Assay Kit (Life Technologies, USA) following the manufacturer's protocol. PCR products were pooled in equimolar amounts into libraries of 48 samples each, and sequenced on an Illumina MiSeq platform in 300 bp paired-end mode at GATC Biotech (Constance, Germany). Data processing and statistical analysis Initial analysis of raw 16S rRNA gene sequencing data was performed using the NG-Tax pipeline (Ramiro-Garcia et al.2016). Sequences were separated into sample-specific bins based on the barcodes, after initial filtering of paired-end libraries to contain only read pairs with perfectly matching barcodes. Operational taxonomic units (OTUs) were defined using an open reference approach, and taxonomy was assigned using a SILVA 16S rRNA gene reference database (Quast et al.2013). Microbial composition plots were generated using a workflow based on Quantitative Insights Into Microbial Ecology (QIIME) v1.9.1 (Caporaso et al.2010). Reads assigned to OTUs of plant origin such as chloroplast and plant mitochondrial DNA were removed from the dataset used for downstream analyses. OTU counts were normalised using cumulative sum scaling (Paulson et al.2013). To get an overview of species composition, a normalised OTU matrix was exported to Microsoft Excel, and the relative contribution based on normalised OTU numbers per taxa was calculated. Median values of taxa relative abundance in a group of samples were used to compare groups with each other (e.g. male vs female). The OTU matrix was filtered to exclude OTUs that were present only in a small number of samples. More specifically, for each dataset, OTUs were removed that were present in less than five samples (50% of the smallest group size). Measures of alpha and beta diversity and initial multivariate analysis using principal coordinate analysis (PCoA) were performed on the rarefied matrix (depth–1650 observations) using weighted and unweighted UniFrac as distance measures as implemented in QIIME. Significance of differences in relative abundances of OTUs between individual samples was determined using Kruskal-Wallis tests when comparing more than two groups, and non-parametric t-tests with 500 Monte Carlo permutations in case of comparisons of two groups, using normalised, summarised and filtered OTU tables. False discovery rate (FDR) correction of P-values was used to reduce the chance of type I statistical errors, when multiple statistical hypotheses were tested. To identify strength and statistical significance of sample groupings with weighted and unweighted UniFrac as distance measures, we used the ‘adonis’ test as implemented in the R package ‘vegan’. Canoco 5.0 was used for multivariate statistical analysis and visualisation of correlations between microbial composition of samples and explanatory factors. Redundancy analysis (RDA) was performed as described previously (Šmilauer and Lepš 2014). As input dataset for RDA, we used the taxonomy summary table at genus level after removal of taxa which were found in less than eight samples (applied for each dataset individually) with addition of the corresponding sample metadata. No transformation or normalisation of the data was done. The significance of observed community variations was evaluated using a Monte Carlo permutation test. To identify microbial species most strongly correlated with investigated factors such as area or host species, we used the LefSe (Linear discriminant analysis of Effect Size) algorithm for biomarker identification (Segata et al.2011). Data were processed using tools developed by the Huttenhower laboratory implemented in the Galaxy environment (www.huttenhower.sph.harvard.edu/galaxy/). Preparation of input data and analysis were performed according to the standard workflow, using default settings (0.05—alpha value for the factorial Kruskal-Wallis test among classes; threshold on the logarithmic LDA score for discriminative features was 2.0; the strategy for multiclass analysis was ‘all-against-all’). After raw data processing and initial analysis, samples were organised into five sets (samples are used in several sets), allowing us to perform separate analyses and statistics while focusing on a particular research question: ‘All_samples’—initial set of all 138 samples; ‘E_fulvus’—samples obtained from E. fulvus (n = 45) from three different areas (Andasibe n = 14; Ankarafantsika n = 21; Nosy Tanikely n = 10); ‘Kirindy’—set of samples from E. rufifrons collected in Kirindy NP (n = 44); ‘Ranomafana’—samples from E. rufifrons (n = 27) and E. rubriventer (n = 22) collected in Ranomafana NP (total n = 49); ‘E_rufifrons’—samples from E. rufifrons collected from two different areas (Ranomafana NP n = 27 and Kirindy NP n = 44). Availability of data and materials Datasets generated in this study are available in the public read archive EBI, study name ‘Area of habitation strongly influences faecal microbial composition of wild lemurs’, with accession number PRJEB20007. RESULTS In this study, we analysed the faecal microbiota of a total of 138 individuals belonging to three different Eulemur species, using Illumina MiSeq sequencing of PCR-amplified 16S rRNA gene fragments covering the V1-V2 variable region. In total, we obtained 6220 515 reads, ranging from 1652 to 178 522 reads per sample (r/s) with a median of 36 092 r/s. Obtained reads were assigned to 1053 OTUs using NG-Tax, an in-house developed pipeline (Ramiro-Garcia et al. 2016). Across all samples, OTUs belonged to 12 bacterial phyla, i.e. Acidobacteria, Actinobacteria, Bacteroidetes, ‘Candidatus Saccharibacteria’ (TM7), Cyanobacteria, Firmicutes, Lentisphaerae, Proteobacteria, Spirochaetes, Synergistetes, Tenericutes and Verrucomicrobia. The fraction of non-assigned to any taxonomic level OTUs varied from 2.6% to 16.7% with an average of 9.3 ± 2.5% [s.d.] in all analysed samples. Predominant phyla, regardless of lemur species and sampling location, were Firmicutes 43.3 ± 6.4%, Bacteroidetes 30.3 ± 5.3%, Cyanobacteria 5.2 ± 3.3% and Proteobacteria 7.4 ± 3.1%. At the genus level, a total of 59 taxa were identified, 15 of which had an average relative abundance across all samples of more than 1% and comprised more than 80% of all sequences. Overall, 34% of all OTUs could be assigned at genus level and 63% at family level. Phylogenetic clustering based on relative abundance at the genus level showed that the most abundant genera could be clustered into two groups: one group consisting of the two most abundant genera (unidentified genus (UG_1 (Clostridiales) 24.9% ± 5.4%, UG_2 (Bacteroidales) 14.9 ± 3.6%) and another group consisting of another five genera (UG_3 (Prevotellaceae) 7.7 ± 2.3%, UG_4 (Cyanobacteria) 5.2 ± 3.3%, UG_5 (Bacteroidaceae) 4.6% ± 2.4%, UG_6 (Lachnospiraceae) 4.9% ± 2.8%, UG_7 (Ruminococcaceae) 4.2% ± 2%) (Fig. S1; Table S2, Supporting Information. Taxonomy of UGs). Figure 2. View largeDownload slide Differences in bacterial composition between lemur species (E. fulvus, E. rubriventer, E. rufifrons); (A) relative abundance at the phylum level of faecal microbiota of the different lemur species using dataset ‘All_samples’; (B) taxa identified by LefSe as potential biomarkers for the discrimination of studied Eulemur species; (C) taxa identified by LefSe as potential biomarkers for the discrimination of faecal samples taken in Ranomafana NP and Kirindy Forest. LDA—linear discriminant analysis. Figure 2. View largeDownload slide Differences in bacterial composition between lemur species (E. fulvus, E. rubriventer, E. rufifrons); (A) relative abundance at the phylum level of faecal microbiota of the different lemur species using dataset ‘All_samples’; (B) taxa identified by LefSe as potential biomarkers for the discrimination of studied Eulemur species; (C) taxa identified by LefSe as potential biomarkers for the discrimination of faecal samples taken in Ranomafana NP and Kirindy Forest. LDA—linear discriminant analysis. Effect of lemur species on faecal microbiota composition In order to address the influence of the host species on observed variation in microbial composition, two different datasets (see information regarding the composition of datasets in Experimental Procedures) were analysed:‘All_samples’, i.e. the entire dataset of 138 samples, and ‘Ranomafana’, the latter of which allowed us to minimise the influence of explanatory variables other than host species. Analysis of ‘All_samples’ revealed that samples from E. fulvus (n = 45) had significantly lower alpha diversity (P = 0.003) in comparison with samples taken from E. rufifrons (n = 22) (Fig. S2, Supporting Information). The relative abundance of Proteobacteria, Lentisphaerae, Synergistetes, Bacteroidetes, Actinobacteria, Firmicutes, Tenericutes and Candidate division TM7 differed among studied lemur species (FDR-corrected P < 0.012, Fig. 2A). At the genus level, 26 taxa differed in relative abundance (corrected P < 0.036), with some present only within one lemur species. The genera Anaeroplasma and UG_8 (Desulfovibrionaceae) were found only in samples belonging to E. rufifrons, albeit with relative abundances below 1%. Eight genera were identified by LefSe analysis as potential biomarkers of the different lemur species: Bacteroides and Phascolarctobacterium were identified as microbial biomarkers of E. rufifrons; UG_6 (Lachnospiraceae), Campylobacter, UG_9 (Synergistales), UG_10 (Clostridiales) and UG_11 (Xanthomonadales) were biomarkers for E. rubriventer; and UG_12 (Anaeroplasmatales) was associated with E. fulvus (Fig. 2B). No clear grouping of samples by lemur species was observed based either on weighted or unweighted UniFrac distance. This was confirmed by the ‘adonis’ test that revealed only a weak linear correlation between samples, with an R2 of 0.11 and 0.13 for weighted and unweighted distances, respectively. RDA with lemur species as the only explanatory variable showed that this variable significantly (P = 0.008) contributed to the observed variation in faecal microbiota composition (Fig. 3A). Furthermore, when comparing the faecal microbiota of E. rufifrons and E. rubriventer in Ranomafana NP, RDA showed that although ‘lemur species’ was a significant explanatory variable (P = 0.024), it only explained 6.8% of the observed variation in microbial community composition (Fig. 3B). All phyla observed in the full dataset were also present in Ranomafana NP, and only the phylum Firmicutes showed nearly significant differences in relative abundance between E. rufifrons (44.3 ± 7.1%) and E. rubriventer (39.7 ± 4.7%) (P = 0.008; corrected P = 0.1). At the genus level, UG_13 (Porphyromonadaceae), UG_5 (Bacteroidaceae) and UG_19 (Bacillales) differed in relative abundance when comparing microbial composition in both lemur species (corrected P < 0.04). Similar to dataset ‘All_samples’, no separation or grouping was observed among samples from the ‘Ranomafana’ dataset in weighted or unweighted UniFrac matrix-based PCoA plots (R2 = 0.06 and R2 = 0.05 for unweighted and weighted distance matrices, respectively; data not shown). Figure 3. View largeDownload slide Ordination triplots based on RDA with lemur species as explanatory variables. (A) In dataset ‘All_samples’, 7.2% of the variation is captured by the first two canonical axes; (B) in dataset ‘Ranomafana’, 53.7% of variation is captured by the first two canonical axes, and both host species significantly (P = 0.002) contributed to explaining the observed variation in microbiota composition. Figure 3. View largeDownload slide Ordination triplots based on RDA with lemur species as explanatory variables. (A) In dataset ‘All_samples’, 7.2% of the variation is captured by the first two canonical axes; (B) in dataset ‘Ranomafana’, 53.7% of variation is captured by the first two canonical axes, and both host species significantly (P = 0.002) contributed to explaining the observed variation in microbiota composition. Figure 4. View largeDownload slide PCoA three-dimensional (first tree PCoA axes) plots based on weighted (A, C, E) and unweighted (B, D, F) UniFrac distance matrices. Samples are represented by dots, colour-coded by sampling location. Plots contain all samples (A, B), or are species specific (E. fulvus—C, D; E. rufifrons—E, F). Figure 4. View largeDownload slide PCoA three-dimensional (first tree PCoA axes) plots based on weighted (A, C, E) and unweighted (B, D, F) UniFrac distance matrices. Samples are represented by dots, colour-coded by sampling location. Plots contain all samples (A, B), or are species specific (E. fulvus—C, D; E. rufifrons—E, F). Variation of the lemur microbiota in contrasting regions of Madagascar To determine to what extent occupancy (i.e. area of habitation) can explain the observed variation in faecal microbiota, we analysed three datasets: ‘All_samples’, ‘E_fulvus’ and ‘E_rufifrons’. Dataset ‘All_samples’ allowed us to identify the influence of ‘area’ among all other variables such as host species, sex and age. With dataset ‘E_fulvus’, we focused on a single lemur species, E. fulvus, which was sampled at three different locations, allowing us to more specifically address variation in faecal microbiota composition found within one species exposed to different environmental conditions. We constructed dataset ‘E_rufifrons’ for the same purpose as dataset ‘E_fulvus’, taking into account that individuals of E. rufifrons were sampled in two distinct areas. Dataset’ All_samples’(n = 138) showed that samples from Nosy Tanikely (n = 10) area had lower alpha diversity as compared to all other areas (P = 0.01, Fig. S3, Supporting Information). The relative abundances of 10 out of 12 phyla, except for Bacteroidetes and Acidobacteria, were different (corrected P < 0.026) between sampling sites. It should be noted that members of the phylum Acidobacteria were found in only a few samples (10 out of 138). Furthermore, at the genus level 54 out of 59 taxa that were observed in more than five samples showed significant differences in relative abundance between areas (corrected P < 0.05). Among these genera, some were found only within Kirindy NP, namely Anaeroplasma, Rhizobium and UG_8 (Desulfovibrionaceae). Members of the genus Bacteroides, UG_13 (Anaeroplasmatales), UG_17 (Clostridiales) order and UG_30 (Anaeroplasmataceae) were found exclusively in samples from the relatively dry areas (Kirindy and Ankarafantsika). The most abundant genus (UG_1 (Clostridiales)) did not vary significant between areas. Samples showed a slight visual grouping according to area in PCoA plots based on weighted UniFrac distances (R2 = 0.29), with better group separation being observed in the case of unweighted UniFrac (R2 = 0.34) (Fig. 4A and B). Furthermore, constrained analysis (RDA) showed that all areas included as explanatory variable significantly (P < 0.05) contributed to the observed variation in faecal microbiota composition (Fig. 5A). Figure 5. View largeDownload slide Ordination triplot based on RDA with areas of sampling as explanatory variables. (A) In dataset ‘All_samples’, 14.2% of the variation was captured by the first two canonical axes, and a statistically significant (P < 0.05) effect of areas as explanatory factors was observed; (B) in dataset ‘E_fulvus’, 35.7% of variation was captured by the first two canonical axes, and all three areas significantly contributed to explaining the observed variation in microbiota composition (P = 0.002); (C) in dataset ‘E_rufifrons’, 67.6% of variation was captured by the first two canonical axes. Explanatory variable ‘areas’ was statistically significant (P = 0.004) as a conditional effect. Figure 5. View largeDownload slide Ordination triplot based on RDA with areas of sampling as explanatory variables. (A) In dataset ‘All_samples’, 14.2% of the variation was captured by the first two canonical axes, and a statistically significant (P < 0.05) effect of areas as explanatory factors was observed; (B) in dataset ‘E_fulvus’, 35.7% of variation was captured by the first two canonical axes, and all three areas significantly contributed to explaining the observed variation in microbiota composition (P = 0.002); (C) in dataset ‘E_rufifrons’, 67.6% of variation was captured by the first two canonical axes. Explanatory variable ‘areas’ was statistically significant (P = 0.004) as a conditional effect. Figure 6. View largeDownload slide Heat map of relative microbial abundance at genus level in the dataset ‘E_fulvus’, with samples placed on the Y-axis, and genera with relative abundance of more than 2.5% across the dataset on the X-axis. The yellow colour indicates high relative abundance values; dark blue indicates low relative abundance values. Sample clustering and dendrogram were produced using the ‘Bray’ method as implemented in the ‘vegan’ R package. Side bar (left) indicates area of sample collection: green—Nosy Tanikely, brown—Andasibe NP, blue—Ankarafantsika. Figure 6. View largeDownload slide Heat map of relative microbial abundance at genus level in the dataset ‘E_fulvus’, with samples placed on the Y-axis, and genera with relative abundance of more than 2.5% across the dataset on the X-axis. The yellow colour indicates high relative abundance values; dark blue indicates low relative abundance values. Sample clustering and dendrogram were produced using the ‘Bray’ method as implemented in the ‘vegan’ R package. Side bar (left) indicates area of sample collection: green—Nosy Tanikely, brown—Andasibe NP, blue—Ankarafantsika. In dataset ‘E_fulvus’, 8 out of 12 phyla (Tenericutes, Cyanobacteria, Spirochaetes, Lentisphaerae, Firmicutes, Candidate division TM7, Proteobacteria, Actinobacteria) differed in relative abundance between areas (corrected P < 0.01). In line with the extensive differences observed at the phylum level, 34 out of 55 detected genera showed significant differences in relative abundance between areas. UG_17 (Clostridiales), UG_35 (Enterobacteriales), UG_29 (Spirochaetaceae), UG_12 (Anaeroplasmatales) and the genus Bacteroides were found only in the samples from Ankarafantsika NP. UG_31 (Veillonellaceae) was found exclusively in the samples from Andasibe NP, and Mesorhizobium only in Nosy Tanikely. Several genera were absent in one out of three areas: UG_19 (Bacillales), Helicobacter and Thalassospira were not detected in samples taken in Nosy Tanikely, whereas Pseudobutyrivibrio was not found in Andasibe. Samples in the ‘E_fulvus’ dataset clustered into three groups, correlating with the three different sampling sites based on the relative abundance of bacterial genera (Fig. 6). Furthermore, samples formed separated groups in PCoA plots based on weighted and unweighted UniFrac distances (R2 = 0.35 and R2 = 0.41, respectively; Fig. 4C and D). RDA showed that among all factors only area significantly (P = 0.002) contributed to explaining observed differences in faecal microbial composition of E. fulvus, with Ankarafantsika having the highest explanatory value (19.9%) (Fig. 5B). In dataset ‘E. rufifrons’, the relative abundance of the phyla Actinobacteria, Candidate division TM7 and Proteobacteria was different between the two locations where this lemur species was found, i.e. Kirindy NP and Ranomafana NP (corrected P < 0.009). In total, 26 genera were different in relative abundance when comparing both locations (corrected P < 0.047). Members of the genus Bacteroides were completely absent from the samples collected from Ranomafana NP, while their mean relative abundance was 1.4 ± 1% among samples collected from Kirindy NP. All remaining genera that were found exclusively in samples taken in Kirindy NP (UG_8 (Desulfovibrionaceae), UG_12 (Anaeroplasmatales), Anaeroplasma) had relative abundances below 0.3%. The genus Phascolarctobacterium was found in all samples from Kirindy NP (mean abundance 3 ± 1.3%), but only in 12 out of 27 samples (average abundance 0.2 ± 0.2%) from Ranomafana NP. UG_21 (Bacteroidales, 0.2 ± 0.4%), Mesorhizobium (0.1 ± 0.2%), UG_33 (Xanthobacteraceae, 0.4 ± 0.2%), Stenotrophomonas (0.1 ± 0.2%), UG_32 (Rhizobiales, 0.3 ± 0.9%) and UG_10 (Clostridiales, 0.1 ± 0.1%) were found exclusively in Ramonofana NP, albeit not in all samples collected in that area, and at low relative abundance. Additional differences included UG_5 (Bacteroidaceae) (6.2 ± 2.2% in Kirindy NP vs 3.3 ± 2% in Ranomafana NP) and UG_6 (Lachnospiraceae) (6.2 ± 2.3% in Ranomafana NP vs 2.9 ± 1.9% in Kirindy NP). Twenty one genera for Ranomofana NP and 12 for Kirindy NP were identified by LefSe as microbial ‘biomarkers’ (Fig. 2C). Multivariate analyses supported the separation of samples according to sampling location, with a clear grouping being observed in PCoA plots based on both weighted (R2 = 0.19) and unweighted (R2 = 0.23) UniFrac distance matrices (Fig. 4E and F). Furthermore, RDA showed that from all explanatory variables (area, age, sex) only variables in the group ‘area’ significantly contributed to explaining the observed variation in microbial community composition (P < 0.004), with both areas (Kirindy and Ranomafana) explaining 8.2% of variation (Fig. 5C). Variation in microbiota composition associated with sex and age of animals Influence of sex and age on faecal microbiota composition was investigated using datasets ‘All_samples’ and ‘Kirindy’. No significant differences were observed in microbiota between samples collected from males and females in any of the datasets at phylum or genus level, and no grouping was observed with multivariate analysis using PCoA and RDA (data not shown). Similarly, no significant variation in relative abundance of detected phyla was observed between age groups, and only the relative abundance of the genus Phascolarctobacterium (corrected P = 0.01) differed in age groups in dataset ‘All_samples’. Multivariate analysis (RDA and PCoA) confirmed that age did not significantly contribute to explaining the variation in faecal microbiota composition. DISCUSSION In the current study, we characterised the faecal microbiota of three frugivorous Eulemur species, and assessed to what extent the naturally occurring variation in intestinal microbiota composition is associated with occupancy, species, age and sex of individuals. Findings presented here showed that the gut microbial community of these animals is dominated by members of the phylum Firmicutes and to a lesser extent Bacteroidetes. It has been previously reported that predominance of Firmicutes or Bacteroidetes is different among animal species and mostly correlated with dietary mode and taxonomic lineage of a given species (Ley et al.2008). Our results confirmed the high proportion of Firmicutes that was previously observed in other species of frugivorous and omnivorous primates (Ochman et al.2010; Gomez et al.2015), including humans, despite the fact that phylogenetically lemurs are one of the most distinct and ancient groups within the primates (Ni et al.2013). Notably, human studies showed that the Firmicutes to Bacteroidetes ratio is not static and can be largely influenced by the presence of carbohydrates in the diet, although it is not clear which of the two phyla has a leading role as key degrader of complex carbohydrates in the human intestine. For example, on one hand an increase in relative abundance of Firmicutes was correlated with consumption of whole grains and total carbohydrate intake (Martínez et al.2013), and several species belonging to this phylum are viewed as key degraders of resistant starch (Ze et al.2012). On the other hand, it has been shown that the depletion of Firmicutes and increase in Bacteroidetes in African children from a rural area in comparison with European children were related to consumption of a traditional African diet rich in fibres and polysaccharides (De Filippo et al.2010). Such seemingly conflicting evidence might be related to the high phylogenetic and functional diversity within both phyla, including a large number of fibre- and carbohydrate-degrading species. Consequently, one can speculate that specific aspects of the diet of lemurs will result in a shift of the Firmicutes/Bacteroidetes ratio, the direction of which might not be predictable based on general characteristics of the diet. Our study showed that a large fraction of Firmicutes associated sequences was assigned to a single genus-level taxon, UG_1 (Clostridiales), accounting for 24.9 ± 5.7% of the total bacterial community. It is tempting to speculate that members of this genus have an important role in intestinal function of the three Eulemur species covered in our study; however, due to lack of physiological and ecological data, conclusions regarding their function and place in gut ecology remain speculative, awaiting isolation and/or (meta)omics analyses (Gutleben et al.2018). Notably, members of the Proteobacteria had relatively high abundance (7.4 ± 3.1%) in all three studied Eulemur species. In humans, high relative abundance of this phylum (9.7%–14.9%) has been associated with gastric bypass, metabolic disorders, inflammation and cancer, whereas its relative abundance in healthy individuals amounts to only about 4.5% (Shin, Whon and Bae 2015). On the other hand, previous research showed host species related differences in abundance of Proteobacteria among primates. For instance, Bello González et al. (2015) observed that faecal samples of humans and chimpanzees had similar relative abundances of Proteobacteria (1% and 1.2%, respectively), whereas in gorilla samples this phylum reached a relative abundance of 7% (Bello González et al.2015). Furthermore, a similar relative abundance of Proteobacteria (9.1%) was reported in faecal microbiota of L. catta (McKenney, Rodrigo and Yoder 2015). Hence, we suggest that the high relative abundance observed in this study is not necessarily a sign of a health problem of the investigated population of lemurs, but rather a feature of the normal microbial composition of frugivorous lemurs. We found that on average 5.2 ± 3.3% of all reads were assigned to OTUs belonging to the Cyanobacteria phylum. Latest research shows that members of this phylum are indeed a genuine part of the human intestinal microbiota (Di Rienzi et al.2013). Furthermore, the presence of this phylum was observed in previous studies which characterised the intestinal microbial composition of other primates, including L. catta (McKenney, Rodrigo and Yoder 2015). Furthermore, we found that 66% of genus-level taxa could not be confidently classified to a particular genus in the Silva v111 database, including several of the most predominant taxa. This is in line with the limited attention that the intestinal microbiota of lemurs has received to date, and hence there is a lack of knowledge regarding specific taxa present in the intestine of these animals. Research on intestinal microbiota of other poorly studied animals showed similar findings. Roggenbuck et al. (2014) found that only 28% of the observed genera in the giraffe rumen could be assigned to known taxa (Roggenbuck et al.2014). Similar observations have even been made for less well-characterised human populations. In a recent study, Schnorr et al. (2014) found that 22% of the total microbial community of central Tanzanian Hadza individuals could not be assigned at family and genus level, whereas this was not the case for the Italian control population (Schnorr et al.2014). To assess the role of different natural environmental factors in shaping the intestinal microbiota and how these factors relate to the influence of the different host species, we divided samples into subsets. This approach allowed us to gain a better insight into the effect of specific factors on microbiota composition, in addition to a more generic analysis of all factors at the same time in a relatively heterogeneous dataset. The value of this tiered approach was confirmed by the fact that for all factors of interest, a first insight into potential effects that could be obtained with the whole dataset received additional, more robust, support from the analysis of specific subsets of samples. It should also be noted that, due to the nature of wildlife sampling under natural conditions, it remains challenging to obtain balanced sample sets with equal numbers of samples in each group. We discovered that in samples analysed in our study, the most influential factor contributing to shaping microbiota composition was the area of sampling. When we applied PCoA based on either weighted or unweighted UniFrac distances, separation into areas was obvious in all datasets when included. Remarkably, clustering of samples was tighter with better separation of samples when unweighted UniFrac was used as a distance measure. This observation suggests, taking into account the nature of the UniFrac distance calculation, that the faecal microbiota of lemurs from different areas is more distinct with respect to microbial species composition than in relative abundance of prevalent taxa. Constrained multivariate analysis (RDA) confirmed that occupancy is the most influential explanatory variable with respect to the observed variation in lemur intestinal microbiota composition. Madagascar is known to have different environmental conditions and biodiversity within relatively small areas (Wilmé, Goodman and Ganzhorn 2006). Furthermore, sampling locations were positioned with considerable distance from each other, and were characterised by major climatic differences such as amount of precipitation and forest composition. Hence, it is likely that the availability of food items during the year, in particular fruits and flowers, is the driving force that leads to differences in the intestinal microbiota. These food items are scarce during the dry season in the dry deciduous forest areas such as Kirindy NP (Sorg and Rohner 1996; Ganzhorn 2002), whereas in the areas characterised by wet rainforests such as Ranomafana NP and Andasibe NP, as well as in Nosy Tanikely, these food items are abundantly available almost all year around. Surprisingly, relatively low microbial alpha diversity was observed in the tropical rainforest. It should be noted, however, that samples from this area only included E. fulvus. It is tempting to speculate that the lower alpha diversity in E. fulvus might be explained by adaptation of the microbiota to sugar-rich food that is prevalent in Nosy Tanikely due to the abundant presence of mango trees. We also observed differences in microbiota composition related to the species of lemurs; however, these differences were secondary to those observed between different areas of habitation. Host genetic differences are among the major driving forces shaping intestinal microbiota composition (Khachatryan et al.2008), including different animal species with similar dietary habits (Moeller et al.2014). In our study, this was not the case. One explanation for this finding could be that this study has been conducted on congeneric species which by definition are genetically close (Markolf and Kappeler 2013), and have almost identical digestive systems (Tattersall 1993), leading to the moderate effect of genetics on gut microbiota composition observed here. We did not find any evidence for an influence of sex or age on microbiota composition. It should be noted, however, that in this study different age classes were not equally represented in the datasets, as we mostly sampled adult individuals (74.6% of all samples). Furthermore, all of the non-adults were at or beyond juvenile stage. Based on knowledge about microbiota development of human infants, the transformation of microbiota to an adult-like mature, composition occurs before reaching the juvenile stage (Koenig et al.2011; Wopereis et al.2014). Surprisingly, we did not find any differences in beta- and alpha diversities of microbiota between males and females. Many studies showed an influence of this factor in different species (Bolnick et al.2014), including lemurs (Aivelo, Laakkonen and Jernvall 2016). However, it was pointed out before that other factors, such as host genetics, can outweigh the influence of this factor (Kovacs et al.2011). Hence, it is tempting to speculate that in the present study the effect of sex on the faecal microbiota of the different Eulemur species might have been obscured by more influential factors as well as relatively large variation in microbial composition between individual animals. In conclusion, we showed that intestinal microbiota in three genetically close species of lemurs was most strongly influenced by their occupancy, whereas the influence of genetic differences was minor, and influence of sex and age was not detectable. All three lemur species had similar bacterial composition in terms of predominant and prevalent bacterial taxa. The findings reported here contribute to our knowledge about the intestinal microbiota in non-human primates, and factors that shape the bacterial composition in wild lemur populations, which can be extrapolated into general rules of intestinal microbiota assembly. Furthermore, the high fraction of poorly assigned taxa reinforces the notion that microbiota of non-humanoid primates has so far received little attention, harbouring a broad range of potentially novel bacterial species and genera that deserve attention in future studies. SUPPLEMENTARY DATA Supplementary data are available at FEMSEC online. Acknowledgements For their advice in the field, the authors want to thank Dr P. Wright, Dr E. Larney and John E. Cadle, and for their logistical support we want to thank the Centre Valbio staff (CVB), the Institute for the Conservation of Tropical Environments (ICTE) and the Madagascar Institute pour la Conservation des Ecosystèmes Tropicaux (MICET), especially B. Andriamihaja, who facilitated arrangements on necessary permits. Thanks go to the Malagasy government and Madagascar National Parks (MNP), Centre de Formation Professionelle Forestière (CFPF), and Mitsinjo Association for their formal permission. For their fieldwork assistance and indispensable knowledge, we want to thank all field technicians and guides that helped as at the various national parks. FUNDING This work was supported by an ‘Nederlandse Organisatie voor Wetenschappelijk Onderzoek’ and ‘The Graduate School for Production Ecology and Resource Conservation’ grant [Grant nr. 1208; project nr. 5120873001], and the European Union Erasmus-Mundus PhD fellowship [EM Action 2- Partnership IAMONET- 20132520]. Other necessary funds to conduct the fieldwork of this study were provided by the Koninklijke Nederlandse Akademie van Wetenschappen (KNAW) Academy Ecology Fund, the Dr J. L. Dobberke Foundation, the Treub Foundation, the Dutch Fund for Research on Nature Conservation, the Dutch Royal Botanical Society and the Dutch Royal Zoological Society. Conflict of interest. None declared. REFERENCES Aivelo T, Laakkonen J, Jernvall J. Population- and individual-level dynamics of the intestinal microbiota of a small primate. Appl Environ Microb  2016; 82: 3537– 45. Google Scholar CrossRef Search ADS   Amato KR, Leigh SR, Kent A et al.   The gut microbiota appears to compensate for seasonal diet variation in the wild black howler monkey (Alouatta pigra). Microb Ecol  2015; 69: 434– 43. Google Scholar CrossRef Search ADS PubMed  Atsalis S. Spatial distribution and population composition of the brown mouse lemur (Microcebus rufus) in Ranomafana National Park, Madagascar, and its implications for social organization. Am J Primatol  2000; 51: 61– 78. Google Scholar CrossRef Search ADS PubMed  Balko EA, Brian Underwood H. Effects of forest structure and composition on food availability forVarecia variegata at Ranomafana National Park, Madagascar. Am J Primatol  2005; 66: 45– 70. Google Scholar CrossRef Search ADS PubMed  Bello González T, van Passel M, Tims S et al.   Application of the Human Intestinal Tract Chip to the non-human primate gut microbiota. Benef Microbes  2015; 6: 271– 6. Google Scholar CrossRef Search ADS PubMed  Bolnick DI, Snowberg LK, Hirsch PE et al.   Individual diet has sex-dependent effects on vertebrate gut microbiota. Nat Commun  2014; 5: 4500. Google Scholar CrossRef Search ADS PubMed  Caporaso JG, Kuczynski J, Stombaugh J et al.   QIIME allows analysis of high-throughput community sequencing data. Nat Methods  2010; 7: 335– 6. Google Scholar CrossRef Search ADS PubMed  Clemente JC, Ursell LK, Parfrey LW et al.   The impact of the gut microbiota on human health: an integrative view. Cell  2012; 148: 1258– 70. Google Scholar CrossRef Search ADS PubMed  Daims H, Brühl A, Amann R et al.   The domain-specific probe EUB338 is insufficient for the detection of all Bacteria: development and evaluation of a more comprehensive probe set. Syst Appl Microbiol  1999; 22: 434– 44. Google Scholar CrossRef Search ADS PubMed  David LA, Materna AC, Friedman J et al.   Host lifestyle affects human microbiota on daily timescales. Genome Biol  2014; 15: R89. Google Scholar CrossRef Search ADS PubMed  De Filippo C, Cavalieri D, Di Paola M et al.   Impact of diet in shaping gut microbiota revealed by a comparative study in children from Europe and rural Africa. P Natl Acad Sci USA  2010; 107: 14691– 6. Google Scholar CrossRef Search ADS   Di Rienzi SC, Sharon I, Wrighton KC et al.   The human gut and groundwater harbor non-photosynthetic bacteria belonging to a new candidate phylum sibling to Cyanobacteria. Elife  2013; 2: e01102. Google Scholar CrossRef Search ADS PubMed  Ellis RJ, Bruce KD, Jenkins C et al.   Comparison of the distal gut microbiota from people and animals in Africa. PLoS One  2013; 8: e54783. Google Scholar CrossRef Search ADS PubMed  Fogel AT. The gut microbiome of wild lemurs: a comparison of sympatric Lemur catta and Propithecus verreauxi. Folia Primatol  2015; 86: 85– 95. Google Scholar CrossRef Search ADS PubMed  Ganzhorn JU. Distribution of a folivorous lemur in relation to seasonally varying food resources: integrating quantitative and qualitative aspects of food characteristics. Oecologia  2002; 131: 427– 35. Google Scholar CrossRef Search ADS PubMed  Ganzhorn JU, Schmid J. Different population dynamics of Microcebus murinus in primary and secondary deciduous dry forests of Madagascar. Int J Primatol  1998; 19: 785– 96. Google Scholar CrossRef Search ADS   García G, Goodman SM. Hunting of protected animals in the Parc National d’Ankarafantsika, north-western Madagascar. Oryx  2003; 37: 115– 8. Google Scholar CrossRef Search ADS   Gomez A, Petrzelkova K, Yeoman CJ et al.   Gut microbiome composition and metabolomic profiles of wild western lowland gorillas (Gorilla gorilla gorilla) reflect host ecology. Mol Ecol  2015; 24: 2551– 65. Google Scholar CrossRef Search ADS PubMed  Goodman SM, Benstead JP. Updated estimates of biotic diversity and endemism for Madagascar. Oryx  2005; 39: 73– 77. Google Scholar CrossRef Search ADS   Gutleben J, Chaib De Mares M, van Elsas JD et al.   The multi-omics promise in context: from sequence to microbial isolate. Crit Rev Microbiol  2018; 44: 212– 29. Google Scholar CrossRef Search ADS PubMed  Hansen J, Gulati A, Sartor RB. The role of mucosal immunity and host genetics in defining intestinal commensal bacteria. Curr Opin Gastroenterol  2010; 26: 564– 71. Google Scholar CrossRef Search ADS PubMed  Irwin MT, Wright PC, Birkinshaw C et al.   Patterns of species change in anthropogenically disturbed forests of Madagascar. Biol Conserv  2010; 143: 2351– 62. Google Scholar CrossRef Search ADS   Johnson SE, Gordon AD, Stumpf RM et al.   Morphological variation in populations of Eulemur albocollaris and E. fulvus rufus. Int J Primatol  2005; 26: 1399– 416. Google Scholar CrossRef Search ADS   Kabat AM, Srinivasan N, Maloy KJ. Modulation of immune development and function by intestinal microbiota. Trends Immunol  2014; 35: 507– 17. Google Scholar CrossRef Search ADS PubMed  Khachatryan ZA, Ktsoyan ZA, Manukyan GP et al.   Predominant role of host genetics in controlling the composition of gut microbiota. PLoS One  2008; 3: e3064. Google Scholar CrossRef Search ADS PubMed  Koenig JE, Spor A, Scalfone N et al.   Succession of microbial consortia in the developing infant gut microbiome. P Natl Acad Sci USA  2011; 108: 4578– 85. Google Scholar CrossRef Search ADS   Kovacs A, Ben-Jacob N, Tayem H et al.   Genotype is a stronger determinant than sex of the mouse gut microbiota. Microb Ecol  2011; 61: 423– 8. Google Scholar CrossRef Search ADS PubMed  Lewis RJ, Bannar-Martin KH. The impact of cyclone Fanele on a tropical dry forest in Madagascar. Biotropica  2012; 44: 135– 40. Google Scholar CrossRef Search ADS   Ley RE, Hamady M, Lozupone C et al.   Evolution of mammals and their gut microbes. Science  2008; 320: 1647– 51. Google Scholar CrossRef Search ADS PubMed  Markolf M, Kappeler PM. Phylogeographic analysis of the true lemurs (genus Eulemur) underlines the role of river catchments for the evolution of micro-endemism in Madagascar. Front Zool  2013; 10: 70. Google Scholar CrossRef Search ADS PubMed  Martínez I, Lattimer JM, Hubach KL et al.   Gut microbiome composition is linked to whole grain-induced immunological improvements. ISME J  2013; 7: 269– 80. Google Scholar CrossRef Search ADS PubMed  Maukonen J, Saarela M. Human gut microbiota: does diet matter? Proc Nutr Soc  2015; 74: 23– 36. Google Scholar CrossRef Search ADS PubMed  McKenney EA, Rodrigo A, Yoder AD. Patterns of gut bacterial colonization in three primate species. PLoS One  2015; 10: e0124618. Moeller AH, Li Y, Ngole EM et al.   Rapid changes in the gut microbiome during human evolution. P Natl Acad Sci USA  2014; 111: 16431– 5. Google Scholar CrossRef Search ADS   Moeller AH, Caro-Quintero A, Mjungu D et al.   Cospeciation of gut microbiota with hominids. Science  2016; 353: 380– 2. Google Scholar CrossRef Search ADS PubMed  Ni X, Gebo DL, Dagosto M et al.   The oldest known primate skeleton and early haplorhine evolution. Nature  2013; 498: 60– 4. Google Scholar CrossRef Search ADS PubMed  Ochman H, Worobey M, Kuo C-H et al.   Evolutionary relationships of wild hominids recapitulated by gut microbial communities. PLoS Biol  2010; 8: 2638. Google Scholar CrossRef Search ADS   Overdorff DJ. Similarities, differences, and seasonal patterns in the diets of Eulemur rubriventer and Eulemur fulvus rufus in the Ranomafana National Park, Madagascar. Int J Primatol  1993; 14: 721– 53. Google Scholar CrossRef Search ADS   Overdorff DJ. Ecological correlates to activity and habitat use of two prosimian primates: Eulemur rubriventer and Eulemur fulvus rufus in Madagascar. Am J Primatol  1996; 40: 327– 42. Google Scholar CrossRef Search ADS   Overdorff DJ, Strait SG, Telo A. Seasonal variation in activity and diet in a small-bodied folivorous primate, Hapalemur griseus, in southeastern Madagascar. Am J Primatol  1997; 43: 211– 23. Google Scholar CrossRef Search ADS PubMed  Patterson E, Cryan J, Fitzgerald G et al.   Gut microbiota, the pharmabiotics they produce and host health. Proc Nutr Soc  2014; 73: 477– 89. Google Scholar CrossRef Search ADS PubMed  Paulson JN, Stine OC, Bravo HC et al.   Differential abundance analysis for microbial marker-gene surveys. Nat Methods  2013; 10: 1200– 2. Google Scholar CrossRef Search ADS PubMed  Pyritz LW, Kappeler PM, Fichtel C. Coordination of group movements in wild red-fronted lemurs (Eulemur rufifrons): processes and influence of ecological and reproductive seasonality. Int J Primatol  2011; 32: 1325– 47. Google Scholar CrossRef Search ADS PubMed  Quast C, Pruesse E, Yilmaz P et al.   The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res  2013; 41: D590– 6. Google Scholar CrossRef Search ADS PubMed  Rakotonirina G. Composition and structure of a dry forest on sandy soils near Morondava. Primate Rep  1996: 81– 87. Ramiro-Garcia, J, Hermes, GDA, Giatsis C et al.   NG-Tax, a highly accurate and validated pipeline for analysis of 16S rRNA amplicons from complex biomes. F1000Res  2016; 5. Ren T, Grieneisen LE, Alberts SC et al.   Development, diet and dynamism: longitudinal and cross-sectional predictors of gut microbial communities in wild baboons. Environ Microbiol  2016; 18: 1312– 25. Google Scholar CrossRef Search ADS PubMed  Roggenbuck M, Sauer C, Poulsen M et al.   The giraffe (Giraffa camelopardalis) rumen microbiome. FEMS Microbiol Ecol  2014; 90: 237– 46. Google Scholar CrossRef Search ADS PubMed  Salonen A, Nikkilä J, Jalanka-Tuovinen J et al.   Comparative analysis of fecal DNA extraction methods with phylogenetic microarray: effective recovery of bacterial and archaeal DNA using mechanical cell lysis. J Microbiol Methods  2010; 81: 127– 34. Google Scholar CrossRef Search ADS PubMed  Sato H. Seasonal fruiting and seed dispersal by the brown lemur in a tropical dry forest, north-western Madagascar. J Trop Ecol  2013; 29: 61– 69. Google Scholar CrossRef Search ADS   Schnorr SL, Candela M, Rampelli S et al.   Gut microbiome of the Hadza hunter-gatherers. Nat Commun  2014; 5: 3654. Google Scholar CrossRef Search ADS PubMed  Segata N, Izard J, Waldron L et al.   Metagenomic biomarker discovery and explanation. Genome Biol  2011; 12: R60. Google Scholar CrossRef Search ADS PubMed  Shin N-R, Whon TW, Bae J-W. Proteobacteria: microbial signature of dysbiosis in gut microbiota. Trends Biotechnol  2015; 33: 496– 503. Google Scholar CrossRef Search ADS PubMed  Šmilauer P, Lepš J. Multivariate Analysis of Ecological Data Using CANOCO 5 . New York: Cambridge University Press, 2014. Google Scholar CrossRef Search ADS   Sorg J, Ganzhorn J, Kappeler P. Forestry and research in the Kirindy Forest/Centre de Formation Professionnelle Forestière. In: Goodman SM, Benstead JP. (eds). The Natural History of Madagascar . Chicago: The University of Chicago Press, 2003. Sorg J, Rohner U. Climate and tree phenology of the dry deciduous forest of the Kirindy Forest. Primate Rep  1996; 46: 57– 80. Sussman R. Ecological distinctions in sympatric species of Lemur. In: Martin RD. (ed). Prosimian Biology . London: Gerald Duckworth & Co Ltd, 1974. Tattersall I. Speciation and morphological differentiation in the genus Lemur. In: Kimbel WH, Martin LB. (ed). Species, Species Concepts and Primate Evolution . New York: Springer, 1993, 163– 76. Google Scholar CrossRef Search ADS   Tecot SR. Seasonality and predictability: the hormonal and behavioral responses of the red-bellied lemur, Eulemur rubriventer in Southeastern Madagascar, 2008, https://repositories.lib.utexas.edu/bitstream/handle/2152/3947/tecots23568.pdf?sequence=2&isAllowed=y. Tecot SR. It's all in the timing: birth seasonality and infant survival in Eulemur rubriventer. Int J Primatol  2010; 31: 715– 35. Google Scholar CrossRef Search ADS   Tian L, Scholte J, Borewicz K et al.   Effects of pectin supplementation on the fermentation patterns of different structural carbohydrates in rats. Mol Nutr Food Res  2016; 60: 2256– 66. Google Scholar CrossRef Search ADS PubMed  van den Bogert B, de Vos WM, Zoetendal EG et al.   Microarray analysis and barcoded pyrosequencing provide consistent microbial profiles depending on the source of human intestinal samples. Appl Environ Microb  2011; 77: 2071– 80. Google Scholar CrossRef Search ADS   van den Bogert B, Erkus O, Boekhorst J et al.   Diversity of human small intestinal Streptococcus and Veillonella populations. FEMS Microbiol Ecol  2013; 85: 376– 88. Google Scholar CrossRef Search ADS PubMed  Vences M, Köhler J, Glaw F. First record of Mabuya comorensis (Reptilia: Scincidae) for the Madagascan fauna, with notes on the reptile fauna of the offshore island Nosy Tanikely. Museo Regionale di Scienze Naturali . 1997; 15–N1. Wilmé L, Goodman SM, Ganzhorn JU. Biogeographic evolution of Madagascar's microendemic biota. Science  2006; 312: 1063– 5. Google Scholar CrossRef Search ADS PubMed  Wopereis H, Oozeer R, Knipping K et al.   The first thousand days - intestinal microbiology of early life: establishing a symbiosis. Pediatr Allergy Immunol  2014; 25: 428– 38. Google Scholar CrossRef Search ADS PubMed  Wright PC, Erhart EM, Tecot S et al.   Long-term lemur research at Centre ValBio, Ranomafana National Park, Madagascar. In: Kappler PM, Watts DP. (eds). Long-Term Field Studies of Primates . New York: Springer, 2012. Google Scholar CrossRef Search ADS   Yildirim S, Yeoman CJ, Sipos M et al.   Characterization of the fecal microbiome from non-human wild primates reveals species specific microbial communities. PLoS One  2010; 5: e13963. Google Scholar CrossRef Search ADS PubMed  Yu Z, Morrison M. Improved extraction of PCR-quality community DNA from digesta and fecal samples. BioTechniques  2004; 36: 808– 13. Google Scholar PubMed  Ze X, Duncan SH, Louis P et al.   Ruminococcus bromii is a keystone species for the degradation of resistant starch in the human colon. ISME J  2012; 6: 1535– 43. Google Scholar CrossRef Search ADS PubMed  © FEMS 2018. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com

Journal

FEMS Microbiology EcologyOxford University Press

Published: Mar 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial