Productivity and salinity structuring of the microplankton revealed by comparative freshwater metagenomics.

Little is known about the diversity and structuring of freshwater microbial communities beyond the patterns revealed by tracing their distribution in the landscape with common taxonomic markers such as the ribosomal RNA. To address this gap in knowledge, metagenomes from temperate lakes were compared to selected marine metagenomes. Taxonomic analyses of rRNA genes in these freshwater metagenomes confirm the previously reported dominance of a limited subset of uncultured lineages of freshwater bacteria, whereas Archaea were rare. Diversification into marine and freshwater microbial lineages was also reflected in phylogenies of functional genes, and there were also significant differences in functional beta-diversity. The pathways and functions that accounted for these differences are involved in osmoregulation, active transport, carbohydrate and amino acid metabolism. Moreover, predicted genes orthologous to active transporters and recalcitrant organic matter degradation were more common in microbial genomes from oligotrophic versus eutrophic lakes. This comparative metagenomic analysis allowed us to formulate a general hypothesis that oceanic- compared with freshwater-dwelling microorganisms, invest more in metabolism of amino acids and that strategies of carbohydrate metabolism differ significantly between marine and freshwater microbial communities.

Introduction

Lakes are systems of enhanced biological activity and are central to many biogeochemical processes (Battin ; Tranvik ). Lakes also represent a critical natural resource for human societies (Downing ). Although bacteria are known to perform many critical biogeochemical processes and thus also have the potential to modify and control water quality in these ecosystems, we have limited understanding of their functional potential, genetic variability and community interactions. This is partly because most abundant lake bacteria are notoriously difficult to culture in isolation (Newton ). The first sequenced genomes of abundant freshwater bacteria (Garcia ; Hahn ) and recent metagenomic characterization of microorganisms from Lake Gatun (Rusch ), Lac du Bourget (Debroas ) and Lake Lanier (Oh ) have provided some first snapshots of the functional diversity of freshwater bacterioplankton in single lake ecosystems. These studies have corroborated findings based on 16S rRNA amplicon surveys with regards to the composition of freshwater bacterial communities and the existence of a phylogenetically distinct freshwater microbiota (reviewed in Newton ). Nevertheless, because of the often substantial genomic variation among even closely related strains, it is challenging to predict community metabolism solely from taxonomic markers and the often rather limited metabolic and functional information derived from reference isolates.

In contrast with such marker gene approaches, metagenomic analysis has the potential to summarize the combined genetic blueprint of all organisms in a given community (Riesenfeld ). By sequencing all genetic information in a community, the relative abundance of all represented genes can, at least in theory, be determined and used to provide a synoptic description of the functional potential of communities under scrutiny (i.e. Fierer ; Rusch ; Debroas ; Oh ). By annotating and comparing multiple such data sets, differences in the metabolic profiles across environments can furthermore be identified (Dinsdale ), and it is also possible to identify specific genomic adaptations to life in contrasting habitats. Such metagenomic studies have previously revealed significant relationships between the environmental conditions and the functional composition of microbial communities in a wide range of habitats (Tringe ; DeLong ; Dinsdale ; Kunin ; Gianoulis ; Raes ) including a first comparison between metagenomes from freshwater lake and marine samples (Oh ).

Here, we use metagenomic sequence data from marine and freshwater systems to identify general differences in functional gene profiles and the variability in metabolic profiles among lakes of different trophic status. Comparative analyses of freshwater bacterial communities based on taxonomic markers have previously revealed differences in bacterial community composition across trophic gradients, where specific lineages respond either positively or negatively to high productivity (Kolmonen ). Microbial community structure is not only determined by environmental characteristics (Newton ) and contemporary biotic interactions (Eiler ) but also by a complex combination of historical factors such as dispersal limitation, past environmental conditions and evolution (Martiny ). In comparison with oceans, inland waters are much more directly influenced by the surrounding terrestrial landscape and coupled to inputs of organisms and chemical constituents from the catchment. Such external influences are likely to have a profound influence on the phylogenetic composition of bacterioplankton communities (see for example Lindström, 2000; Lindström ; Yannarell and Triplett, 2005; Eiler and Bertilsson, 2007; Newton ; Peura ).

To better understand factors controlling and shaping the community-level functional traits of freshwater microplankton, nine planktonic DNA samples from seven different lakes were analysed by pyrosequencing-enabled metagenomics. In addition, three available freshwater metagenome data sets from National Center for Biotechnology Information-Short Read Archive (NCBI-SRA) were included in the analysis, resulting in a combined freshwater data set from altogether 12 freshwater metagenomes. As marine references, we used 13 marine metagenomes comprising samples from the open and coastal ocean. One further aim was to corroborate that lake systems are not only different from marine systems in their phylogenetic but also in their functional gene composition. By comparing lakes of contrasting productivity, we further aimed at revealing functional differences related to nutrient and energy acquisition as well as substrate preferences.

Results and discussion

General description of the sampling sites and sequence data

DNA samples were collected from seven lakes, whereof two lakes were sampled twice (in Spring and Summer) (Table 1). These nine samples were subject to whole-community genome shotgun 454 pyrosequencing using Titanium chemistry. An additional three freshwater lake metagenomes and 13 marine metagenomes were obtained from public databases. The latter included samples from open-ocean and coastal habitats (Table S1). We selected these 16 metagenomes available at the time of analysis because they were of sufficient size to be compared with our data and processed in the most similar fashion to the nine new freshwater metagenomes with regards to sample handling, DNA extraction, library preparation and sequencing. Still, we want to make the reader aware that they were not processed in an identical way, which might influence the comparison and our interpretations (Carrigg ). In addition, with the limited number of samples and shallow sequencing at hand, we can never cover the entire functional diversity dwelling in both marine and freshwater biomes, and this adds some uncertainties to generalizations of the major findings from this study.

Table 1

Description of lakes used in this study

ID	Sample location	Country	Date	Location	Sample depth	T (°C)	Size fraction (μm)	Habitat type	Tot P
DamariscottaSP	Lake Damariscotta	USA	20090528	44°10′n; 69°29′w	0.5–1	12.1	> 0.2	Mesotrophic lake	10
DamariscottaSU	Lake Damariscotta	USA	20090819	44°10′n; 69°29′w	0.5–1	12.1	> 0.2	Mesotrophic lake	10
Ekoln	Lake Ekoln	Sweden	20070731	59°45′n; 17°36′e	0–2	19.0	0.2–100	Eutrophic lake	50
Erken	Lake Erken	Sweden	20070620	59°25′n; 18°15′e	0–2	18.7	0.2–100	Mesotrophic lake	33
Lanier	Lake Lanier	USA	20090827	34°12′n; 83°59′w	0–5	28.5	0.22–1.6	Mesotrophic lake	30
MendotaSP	Lake Mendota	USA	20090512	43° 6′n; 89°24′w	0.5–1	12.68	> 0.2	Eutrophic lake	118
MendotaSU	Lake Mendota	USA	20090823	43° 6′n; 89°24′w	0.5–1	23.07	> 0.2	Eutrophic lake	100
Spark	Sparkling Lake	USA	20090528	46° 0′n; 89°42′w	0.5–1	13.97	> 0.2	Oligotrophic lake	0.3
Trout	Trout Bog Lake	USA	20090528	46° 2′n; 89°41′w	0.5–1	20.71	> 0.2	Dysotrophic lake	7.8
Vattern	Lake Vättern	Sweden	20070717	58°24′n; 14°36′e	0–2	17.0	0.2–100	Oligotrophic lake	3
Yellowstone1	Yellowstone Lake	USA	20080916	44°28′n; 110°22′w	0–2	46	0.1–0.8	Eutrophic lake	80
Yellowstone2	Yellowstone Lake	USA	20080915	44°28′n; 110°22′w	0–2	12.3	0.1–0.8	Eutrophic lake	80

Tot P, total phosphorus concentration (μg l−1); T, temperature.

The nine lakes included in the analysis represent a wide range of trophic states, including oligotrophic, mesotrophic and eutrophic systems (Table 1). They range from 0.3 to 120 μg l−1 in total phosphorus (TP) and are all situated in the temperate climate zone. On average, 325 000 high-quality reads with mean length of 330 bp were obtained for the lake metagenomes and slightly lower numbers of 280 000 sequences with mean read length of 270 bp for the marine data sets. No quality files were available for the marine data sets, but quality filtering (mean read quality > 21) affected the lake metagenomes very little (1–2% for two data sets, 0% for the majority). To match the quality filtering step as best possible, marine metagenomes had an extra upper length filter added because many long sequences were observed to be of poor quality. Lower length limit (> 150 bp) and clustering to remove artefacts were performed in the same way on all data sets, resulting in over 8.2 million reads in total (Table 2; for detail about the removed sequences in the preprocessing steps, see Table S2). Five samples yielded much lower total sequenced nucleotides than average (84%): marine sample from Sargasso Sea (depth 40 m, 67%) and four lake samples from Yellowstone Lake (sample 1, 79%), Lake Mendota (spring sample, 76%), Trout-Bog Lake (75%) and Sparkling Lake (58%). These samples were also among the most extreme outliers in terms of the eukaryotic content. To ensure robustness of the results, the impact of including/excluding those samples from the statistical analyses was investigated.

Table 2

Characteristics of lake and marine metagenomes

ID	Reference/site	Size after QC (Mb)	Reads after QC	GC content (%)	Isoelectric point	SSU rRNA genes	% Eukaryotic SSU rRNA genes	Reads with STRING hit	Reads with COG assignment	% Reads with COG assignment	Average number of single copy COGs	% Bacteria among single copy COGs	Simple EGS Mb/single copy COG
DamariscottaSP	Martinez-Garcia et al. 2012	121	281 625	48.6	9.75	531 (185/0/5)	2.6	149 906	135 640	42	78	94	1.55
DamariscottaSU	Martinez-Garcia et al. 2012	140	323 939	48.2	9.82	666 (200/0/27)	11.9	156 281	140 701	50	93	93	1.51
Ekoln	this study	115	284 609	46.0	9.49	622 (209/0/13)	5.9	107 593	94 783	33	69	95	1.67
Erken	this study	233	554 862	44.9	9.41	1170 (399/0/5)	1.2	273 058	250 931	45	196	93	1.19
Lanier	Oh et al. 2011	449	1 078 031	47.1	9.77	1989 (714/0/20)	2.7	440 459	399 647	37	252	93	1.78
MendotaSP	Martinez-Garcia et al. 2012	133	319 321	45.7	9.52	1118 (242/0/149)	38.1	124 837	111 654	25	77	97	1.73
MendotaSU	Martinez-Garcia et al. 2012	192	447 054	47.7	9.75	795 (247/0/31)	11.2	173 517	146 222	46	76	95	2.53
Spark	Martinez-Garcia et al. 2012	26	66 160	52.5	10.01	108 (28/0/5)	15.2	22 364	19 857	30	8	87	3.25
Trout	Martinez-Garcia et al. 2012	60	150 515	46.5	9.59	335 (63/0/21)	25.0	46 795	41 628	28	26	88	2.31
Vattern	this study	117	285 637	47.4	9.67	540 (177/0/15)	7.8	116 970	103 047	36	66	93	1.77
Yellowstone1	SRR077348	181	416 139	43.7	9.34	541 (212/0/2)	0.9	152 376	136 972	33	83	93	2.18
Yellowstone2	SRR078855	107	346 239	41.4	9.03	754 (256/1/0)	0.0	91 459	86 132	25	75	97	1.43
FRESHWATER (Mean)		156	379 511	46.6	9.60	764 (244/0/24)	10.2	154 635	138 935	37	92	93	1.91
BATS0	Sargasso Sea	118	478 976	48.0	9.74	1137 (431/0/13)	2.9	142 979	131 449	27	104	97	1.14
BATS200	Sargasso Sea	134	525 891	48.3	9.70	1049 (310/38/13)	3.6	133 259	121 763	23	97	85	1.38
BATS250	Sargasso Sea	115	456 677	46.6	9.63	606 (183/20/9)	4.2	95 919	88 658	19	70	89	1.65
BATS40	Sargasso Sea	95	394 461	48.1	9.78	675 (227/0/17)	7.0	86 262	79 155	20	67	96	1.42
EqDP35155	Equatorial Pacific	56	219 390	45.4	9.70	508 (164/10/3)	1.7	62 135	57 103	26	53	91	1.05
NPTG35179	North Pacific Tropical Gyre	45	181 907	44.8	9.53	656 (253/4/4)	1.5	55 589	51 145	28	45	95	1.00
PNEq35163	Pacific North Equatorial	55	221 925	49.8	9.94	790 (300/6/5)	1.6	59 337	53 915	24	52	92	1.06
PNEqCc35171	Pacific North Equatorial	13	50 267	42.5	9.38	101 (31/3/2)	5.6	15 791	14 620	29	13	92	0.97
SPSG35131	South Pacific Subtropical Gyre	36	155 219	47.7	9.77	583 (225/1/4)	1.7	46 502	42 726	28	39	96	0.94
SPSG35139	South Pacific Subtropical Gyre	16	61 766	41.9	9.33	169 (71/1/0)	0.0	23 083	21 352	35	19	97	0.85
SPSG35147	South Pacific Subtropical Gyre	21	80 088	43.2	9.47	259 (97/3/1)	1.0	28 681	26 504	33	25	93	0.83
WChannelApr	Gilbert et al. 2010	102	278 931	39.2	9.01	317 (64/0/6)	8.6	82 819	67 968	24	35	90	2.91
WChannelJan	Gilbert et al. 2010	208	548 680	38.4	8.97	724 (195/17/6)	2.8	180 844	153 475	28	100	74	2.08
MARINE (Mean)		78	281 091	44.9	9.53	583 (196/8/6)	3.2	77 938	69 987	25	55	91	1.33

The isoelectric point represents the average pH at which predicted genes from a specific metagenome carry no electric charge. In column seven, numbers in parentheses represent the number of SSU rRNA genes annotated to Bacteria, Archaea and Eukaryota, respectively.

COG, clusters of orthologous groups; EGS, effective genome size; QC, quality filtering; mb, megabases; SSU rRNA, small subunit of the ribosomal RNA.

In order to investigate the genomic similarity between and within freshwater and marine samples, DNA sequences were first evaluated for features that did not a priori require any taxonomic or functional annotation. Sequences were evaluated for Guanine and Cytosine (GC) content, isoelectric point and amino-acid usage. The GC content of the freshwater metagenome samples was 46.6% on the average (Table 2), ranging from 35% to 60% for the large majority of reads in each of the individual samples (Fig. S1). This was not significantly different to the average GC content of the marine metagenome samples (Wilcoxon test; P = 0.406) where for example the Sargasso Sea samples (46.6–48.6% on the average) had higher GC content than the Western English Channel (below 40%). The isoelectric points were not significantly different (Wilcoxon test; P = 0.624) between freshwater and marine metagenomes using Open Reading Frames (ORFs) of at least 50 aa in length predicted from six frame translation procedures (Table 2). Nor did the inferred amino acid usage differ between marine and freshwater samples [permutational multivariate analysis of variance (PERMANOVA); P = 0.432]. Specifically, we observed no difference in the usage of sulphur-containing amino acids, methionine and cysteine, for which an increased cost could be expected in freshwater environments. Hence, there was no convincing evidence for ‘elemental sparing’, which has been described as adaptive selection pressure on amino acid usage when cellular maintenance costs for protein synthesis are assumed to affect fitness (Bragg and Wagner, 2009).

Taxonomic composition

The microbial diversity captured in the metagenomic sequences from the 25 different metagenomes was analysed using rRNA hidden Markov models (hmm) and tblastx against Search Tool for the Retrieval of Interacting Genes/Proteins (STRING). Identification and analysis of rRNA genes with hmm identified 16 743 small subunit (SSU) rRNA hits, applying an e-value cut-off of 1e-10 for a hit (Table 2). From these, 33% were of bacterial origin, whereas 2.2% and 0.6% were annotated to eukaryotes and archaea, respectively, with the rest being unclassified (64%) using the SILVA database (Quast ) in combination with the naïve Bayesian classifier (Wang ). Two lake metagenomes (Spring sample from Lake Mendota and Trout Bog Lake) had more than 20% of eukaryotic (18S rRNA) reads annotated as mainly algal-derived. Comparing marine and freshwater metagenomes, archaeal 16S rRNA were more common in marine systems (on average, 3.8% of the annotated SSUs in the marine vs. 0.4% in the freshwater metagenomes) when compared with freshwaters where the proportion of eukaryotic 18S rRNA hits was higher (on average, 3.2% of the annotated SSUs in the marine vs. 10.2% in the freshwater metagenomes). Possible explanations are upwelling events at marine sites that may contribute Archaea to surface communities, but also general physicochemical differences between marine and freshwaters could select for the observed patterns. The taxonomic composition of bacteria in each individual sample was also determined by annotating 16S rRNA genes using a custom curated freshwater database (Newton ) (Fig. 1A). Whatever database used, Proteobacteria was the dominant bacterial phylum in all marine metagenomes. Conversely, all but five of the lake metagenomes instead featured Actinobacteria as the most abundant phylum. In marine environments, alpha-Proteobacteria was the dominant class within the Proteobacteria, whereas beta-Proteobacteria were more abundant in freshwaters. Other abundant phyla in the lakes were Verrucomicrobia, Planctomycetes, Cyanobacteria, and Bacteroidetes. Furthermore, we observed significant differences between marine and freshwater metagenomes in community composition analysed at the phylum level (PERMANOVA, R2 = 0.34, P < 0.001). Resolving sequences to a finer taxonomic level (roughly comparable with genus-level) revealed a dominance of previously identified typical freshwater bacteria in the 12 lake samples, including the freshwater SAR11 (LD12), taxa within the Actinobacterial acI lineage (acI-B1, acI-A6, acI-C2) and the beta-Proteobacterial Polynucleobacter (Fig. 1A; Newton ). Moreover, this reflects previously described patterns between systems of different trophic status where dystrophic (humic) systems such as Trout Bog Lake are lacking most typical freshwater taxa (Peura ).

Fig 1

Heatmap of the 20 most abundant typical freshwater taxa (A) in the metagenomics datasets as inferred from their proportion of SSU rRNA gene sequences. Typical freshwater taxa were defined previously using a well-curated freshwater-specific phylogeny (Newton ). (B) Barplot showing taxonomic classification of bacterial reads into phyla based on the best hit to STRING (Franceschini ).

Using the taxonomic annotations of the best tblastx hit to STRING revealed patterns highly similar to that of the SSU rRNA taxonomy where hits to bacteria dominated (on average 92%) over hits matching archaea (2%) and eukarya (7.6%) (Fig. S2). As for most metagenomes, the dominant portion of the reads had no hits (on average 60% for the lake and 71% for the marine data sets) in the STRING database and could thus not be taxonomically assigned. Still, comparing freshwater and marine metagenomes revealed that hits to the bacterial phylum Actinobacteria were more abundant in freshwater metagenomic libraries (on average 31%) compared with marine metagenomes where hits to Proteobacteria, especially alpha-Proteobacteria, were dominant (on average, 38%; see Fig. 1B), thus corroborating observations made at the SSU rRNA level. Other prominent (sub)phyla in the freshwater metagenomes were beta-Proteobacteria (on average 24%), Bacteroidetes (on average 10%), Cyanobacteria (on average 21%), Verrucomicrobia/Chlamydia (on average 3%) and Planctomycetes (on average 2%). Overall, the metagenomic comparison revealed taxonomic distributions as expected from previous studies based on clone libraries (i.e. Zwart ; Eiler and Bertilsson, 2004) and fluorescence in situ hybridization (Glöckner ).

Comparative functional metagenomics between marine and freshwater systems

Functional assignment was made on the basis of the best tblastx cluster of orthologous genes (COGs) hit using an E-value threshold of 1e−10. To assure the best available taxonomic representation, the STRING database was used (Franceschini ), as it comprises over 1000 genomes of bacteria, archaea and eukaryota compared with 66 genomes in the original COG database. The average percentage of the reads that could be annotated (had a COG annotation) was 37% for lake and 25% for marine metagenomes (range 19–50% per sample). The total number of annotations (COGs) per sample ranged from about 14 500 to almost 400 000 (Table 2). The relative abundance of best hits assigned to each major subsystem (orthologous gene classes, OGCs) in the marine versus freshwater system is summarized in Table 3, showing that ‘Amino acid transport and metabolism’ was the dominant OGC.

Table 3

Summary statistics of each OGC and their comparison between freshwater and marine metagenomes

				Wilcoxon test				Permanova
	Median besta (allb)			Besta		Allb		Besta			Allb
		Fresh	Ocean	W	P value	W	P value	F	R2	P value	F	R2	P value
Translation, ribosomal structure and biogenesis	J	9.4%	10.4%	70	0.098	118	0.030	4.7	0.21	0.002	6.1	0.21	0.001
RNA processing and modification	A	0.0%	0.0%	38	0.473	57	0.270	8.1	0.31	0.001	10.3	0.31	0.001
Transcription	K	3.3%	3.1%	9	0.002	11	0.000	13.0	0.42	0.001	17.1	0.43	0.001
Replication, recombination and repair	L	8.4% (8.6%)	6.4%	4	0.000	6	0.000	10.2	0.36	0.001	13.8	0.38	0.001
Chromatin structure and dynamics	B	0.0%	0.0%	54	0.678	75	0.894	3.7	0.17	0.003	3.9	0.14	0.003
Cell cycle control, cell division, chromosome partitioning	D	1.4% (1.3%)	1.5%	66	0.181	107	0.123	11.4	0.39	0.002	15.1	0.40	0.001
Nuclear structure	Y	0.0%	0.0%	40	0.447	71	0.655	NA	NA	NA	NA	NA	NA
Defence mechanisms	V	2.0%	1.7%	20	0.031	25	0.003	5.7	0.24	0.005	7.9	0.26	0.001
Signal transduction mechanisms	T	1.9%	1.2%	1	0.000	1	0.000	12.3	0.41	0.001	16.9	0.42	0.001
Cell wall/membrane/envelope biogenesis	M	6.3%	5.3%	12	0.004	17	0.000	10.0	0.36	0.001	14.2	0.38	0.001
Cell motility	N	0.2%	0.2%	29	0.157	57	0.270	2.2	0.11	0.075	2.5	0.10	0.047
Cytoskeleton	Z	0.1%	0.3%	81	0.010	127	0.007	3.6	0.17	0.023	3.2	0.12	0.021
Extracellular structures	W	0.0%	0.0%	48	NA	78	NA	NA	NA	NA	NA	NA	NA
Intracellular trafficking, secretion and vesicular transport	U	1.3%	1.4%	73	0.057	114	0.052	5.3	0.23	0.002	7.6	0.25	0.001
Posttranslational modification, protein turnover, chaperones	O	5.0%	5.6%	82	0.007	137	0.001	10.4	0.37	0.001	13.0	0.36	0.001
Energy production and conversion	C	9.5%	11.2% (11.1%)	89	0.001	148	0.000	6.6	0.27	0.001	9.2	0.29	0.001
Carbohydrate transport and metabolism	G	5.5%	5.3%	22	0.047	48	0.110	7.2	0.28	0.001	8.8	0.28	0.001
Amino acid transport and metabolism	E	11.7%	13.7%	89	0.001	146	0.000	11.9	0.40	0.001	14.6	0.39	0.001
Nucleotide transport and metabolism	F	3.8%	4.2%	65	0.208	101	0.225	6.9	0.28	0.001	8.0	0.26	0.001
Coenzyme transport and metabolism	H	4.1%	4.9%	96	0.000	156	0.000	7.9	0.31	0.002	10.0	0.30	0.001
Lipid transport and metabolism	I	4.2%	4.4% (4.3%)	69	0.115	106	0.137	8.6	0.32	0.001	11.3	0.33	0.001
Inorganic ion transport and metabolism	P	4.2% (4.1%)	4.3%	58	0.473	101	0.225	13.0	0.42	0.001	16.8	0.42	0.001
Secondary metabolites biosynthesis, transport and catabolism	Q	1.8% (1.7%)	1.5%	27	0.115	51	0.152	10.5	0.37	0.001	13.8	0.38	0.001
General function prediction only	R	10.0%	9.2% (9.1%)	26	0.098	38	0.030	8.8	0.33	0.001	11.7	0.34	0.001
Function unknown	S	5.5%	3.9%	6	0.000	7	0.000	12.2	0.40	0.001	16.0	0.41	0.001

Smaller data set of eight lake and 12 ocean samples, excluding worst quality samples (see Table S1).

Full data set of 12 lake and 13 ocean samples (see Table S1).

Average and standard deviation are derived from the relative fraction of OCGs averaged over all marine and freshwater metagenomes, respectively. P-values and W statistics from Wilcoxon test on the contribution of ORFs to each OGC as well as results from PERMANOVA to test for differences in functional composition between marine and freshwaters using normalized COGs from each OGC.

NA, Not Assessed.

Counts for 35 marker COGs were used to approximate the average effective genome size in freshwater and marine microbial communities. The estimated average effective genome size for freshwaters (1.91) was slightly higher than for the selected marine systems (1.33, Table 2) (Wilcoxon test; P < 0.003) where the latter estimates are similar to previous estimates for marine plankton (Raes ; Quaiser ). These findings corroborate the widespread assumption that small and streamlined genomes are a more common feature of bacterioplankton from oligotrophic sites (Giovannoni ; Grote ) compared with those that reside in more productive waters such as eutrophic freshwater lakes (i.e. lakes Ekoln, Erken and Mendota; see also Oh ). Discrepancies in estimated genome sizes to previously published estimates (Lake Lanier, our estimate 1.78 vs. published 2.2) are most likely due to differences in databases and quality filtering used.

COGs were normalized against best hits to 35 likely essential and single copy COGs (Table S3; Ciccarelli ; Raes ) without taking read length into account prior to statistical analyses. Each of these single copy COGs had, on average, 77 hits in the 25 metagenomes (range 11–279, representing averages from single metagenomes). To assess whether or not each biome had a distinct functional profile, an ordination was conducted using an occurrence matrix of COGs in nonmetric multidimensional scaling (metaMDS function in R; Oksanen ). PERMANOVA (Anderson, 2001) corroborated the visual impression (Fig. S3) of a significant difference in functional beta-diversity between marine and freshwater systems (PERMANOVA; P < 0.001, R2 = 0.34). These differences were maintained even if low-quality metagenomes were excluded (PERMANOVA; P < 0.001, R2 = 0.34), or when only bacterial COGs where analysed (PERMANOVA; P < 0.001, R2 = 0.35) and when specific OGCs were analysed separately (Table 3). The most pronounced difference in the composition of OGCs was observed for the OGC ‘ion transport and metabolism’ and ‘transcription’, whereas the composition within OGCs ‘Cytoskeleton’ and ‘Cell motility’ were the least separated. Moreover, we also looked for proportional differences at the level of OGC by using Wilcoxon test. Overall, OGCs ‘energy production and conversion’ and ‘coenzyme transport and metabolism’ were under-represented in freshwater metagenomes, whereas core functions involved in ‘transcription’, and ‘replication, recombination and repair’ were over-represented when compared with marine samples (Table 3). The higher proportion of the OGC ‘signal transduction’ in freshwater than marine metagenomes suggest that freshwater microbial communities feature more complex interactions and cellular controls that may involve cell-to-cell communication.

This was also reflected in a more detailed analysis based on the Wilcoxon test where all COGs differing in resampled and normalized occurrence between marine and freshwater systems were tested. Out of 707 COGs identified as significantly different in their prevalence between the marine and freshwater metagenomes (P < 0.01 and false discovery rate < 0.027), and 560 significantly different (P < 0.01) when excluding low-quality metagenomes, limited the list to COGs significant for both all and best data sets to 102 COGs that were over-represented in the marine and 295 in the freshwater metagenomes (Fig. 2, Table S4). For example, core functions belonging to ‘transcription’ such as transcriptional regulators, for example arginine repressor (bacterial) was significantly over-represented in lakes (P < 0.001). ‘Replication, recombination and repair’ was represented by numerous transposases, several helicases, and the recombination repair proteins RecF and RecB, which were all significantly over-represented in the lake metagenomes (all bacterial, P < 0.01). Other COGs over-represented in freshwaters were related to a phosphorus starvation-inducible protein phoH (cog1875, P < 0.001), a growth inhibitor (cog2337, P < 0.002) and two response regulators (cog3707, P < 0.001; cog4566, P < 0.001). Homologues to subunits of archaeal polymerases such as COG1311 (archaeal DNA polymerase II, SSU/DNA polymerase delta, subunit B) and COG1933 (archaeal DNA polymerase II, large subunit) were over-represented in marine metagenomes (P < 0.005 and P < 0.001 respectively). With regards to metabolism, differences between freshwater and marine metagenomes were limited to few key enzymes (Fig. 2). Examples of this is the significant over-representation of malate synthase homologues (cog2225, P < 0.001) and isocitrate lyase (cog2224, P < 0.002) in the marine biome, both coding for enzymes with a central function in the glyoxylate cycle. The isocitrate lyase catalyses the cleavage of isocitrate to succinate and glyoxylate, and the malate synthase feeds glyoxylate into the tricarboxylic acid cycle (TCA) via oxalacetate (known as the glyoxylate shunt). This allows microorganisms to utilize simple carbon compounds as a carbon source when complex sources such as glucose are not available. In the absence of available carbohydrates, the glyoxylate cycle permits the synthesis of carbohydrates needed for cell-wall assembly from lipids via acetate. In contrast, reads annotated as being involved in carbohydrate metabolism (i.e. ‘phosphoenolpyruvate-protein kinase’ cog1080, P < 0.001; ‘Fructose-1-phosphate kinase and related fructose-6-phosphate kinase’ cog1105, P < 0.002) seem to be more common in freshwater as compared with marine metagenomes, where such genes were never significantly over-represented. This included galactose-1-phosphate uridylyltransferase (cog1085, P < 0.001) a putative enzyme central to the Leloir pathway involved in the catalyses between galactose and glucose. Another interesting finding was that homologues of enzymes that hydrolyse glycolipids, glycoproteins, lactose and galactosides to monosaccharides such as alpha- (cog3345, P < 0.001) and beta-galactosidases (cog3250, P < 0.004) were over-represented in freshwater metagenomes. Also, other homologues to enzymes catalysing the hydrolysis of glycosidic linkages were over-represented in the freshwaters metagenomes, including chitinase (cog3179, P < 0.001), glycotransferase (cog438, P < 0.001) and glycosidase (cog2723, P < 0.001; cog366, P < 0.001), known to mediate the production of oligosaccharide and monosaccharide from chitin, cellulose and hemicelluloses. This is consistent with a recent finding that the genomes of the abundant acI-B1 taxon of freshwater Actinobacteria are enriched with glycosidase homologues when compared with other bacterial genomes (Garcia ).

Fig 2

Heatmap of COGs showing only those that were either significantly over- (A) and under-represented (B) in freshwater metagenomes when compared with marine metagenomes after resampling and normalization against single-copy core COGs. Significantly over- and under-represented COGs were identified by Wilcoxon test (P < 0.01) when testing all data sets, as well as the best data sets only, and the subsequent estimation of false discovery rate (q < 0.027). These lists are not exhaustive and only include well-characterized COGs. COGs mentioned in the text are indicated. Dendograms from hierarchical cluster analysis based on displayed COGs are shown at the top of each graph.

Moreover, freshwater microbial genomes seem to harbour a higher proportion of certain putative genes involved in transport of sugars such as xylose (cog4213, P < 0.001; cog4214, P < 0.001) and various polysaccharides (cog1134, P < 0.001; cog1682, P < 0.001; cog3833 P < 0.001) (Fig. 2). A similar pattern was also observed for genes involved in transport of peptides (cog410, P < 0.001; cog411, P < 0.001; cog4177, P < 0.002). In contrast, ORFs putatively identified as ATP-dependent amino-acid transporters (cog2113, P < 0.00009; cog4160, P < 0.001; cog4175, P < 0.001; cog4176, P < 0.002; cog4215, P < 0.001; cog4597, P < 0.001) were significantly over-represented in marine metagenomes. We propose that the compositional differences in amino acid and carbohydrate metabolism is a consequence of major differences in the overall composition of organic substrates available for heterotrophs in the respective biomes. Freshwater systems, including the temperate systems of this study, are highly influenced by allochthonous organic matter inputs from the catchment as well as plant-derived polysaccharides (e.g. xylose-containing hemicellulose) inputs from the littoral zone, whereas marine systems are less influenced by organic matter loadings from such terrestrial surroundings and littoral fringe zones and instead rely largely on autochthonous organic matter inputs from plankton rich in proteinaceous materials (Duarte and Cebrián, 1996; Bertilsson and Jones, 2003).

ORFs putatively involved in acquisition of phosphate (cog573/cog581, P < 0.003/0.006; cog1117, P < 0.002) including phosphate uptake regulators (cog704, P < 0.003) and sulphate (cog555, P < 0.007; cog1118, P < 0.001; cog1613, P < 0.002; cog4208, P < 0.001) were mostly over-represented in freshwater genomes. The over-representation of exopolyphospatase (cog248, P < 0.001) and polyphosphate kinase (cog855, P < 0.001) homologues supports the previously recognized role of polyphosphates as a form of phosphorus storage in freshwater environments (Broberg and Persson, 1988; Ilikchyan ). We did not observe any significant differences in nitrogen metabolism and uptake between the marine and often more productive freshwater systems. The previously inferred reliance on potassium instead of sodium for osmoregulation was a typical feature of the freshwater metagenomes as well as a higher representation of reads annotated as cobalt, magnesium and nickel transporter systems (Fig. 2A). In contrast, homologues of zinc and manganese transporters were over-represented in marine metagenomes (Fig. 2B). This confirms previously reported differences in osmoregulatory traits between freshwater and marine microorganisms inferred from comparative metagenomics of microbial communities (Oh ). These findings are also consistent with recent results based on comparisons of 16S rRNA gene libraries (Zwart ; Lozupone and Knight, 2007; Logares ; Newton ) where salinity was suggested to represent a strong environmental barrier for microorganisms. Our results also point to the importance of factors other than salinity, at least when comparing marine and freshwater environments with regards to substrate availability and substrate acquisition. As illustrated earlier, microbial communities in these contrasting biomes seem to have different metabolic capabilities as genes involved in amino acid metabolism were over-represented in marine metagenomes when compared with freshwater metagenomes, and clear differences in the strategies of carbohydrate metabolism were observed.

Comparing functional profiles among freshwater systems

When freshwater functional profiles were analysed by non-metric multidimensional scaling, it was apparent that Sparkling Lake and Trout Bog Lake metagenomes were rather distinct from the others (Fig. 3). This can at least partly be attributed to their high amounts of eukaryotic sequences. An additional non-exclusive explanation may be trophic status: Trout Bog Lake was the only humic (dystrophic) system, whereas Sparkling Lake was the most oligotrophic system in the study. Interestingly, we observed a significant correlation between the overall functional composition and TP, a widely used proxy for ecosystem productivity (Schindler, 1978) (R2 = 0.53, P = 0.029; Fig. 3). The correlation was even more significant if only bacteria were taken into account (R2 = 0.52, P = 0.018). The observation that the functional profile of one metagenome from Yellowstone Lake was very different from the others was probably caused by the proximity of this sample to a thermal vent and the associated higher temperature and different ion composition. For a more detailed analysis, we relied on maximal information-based non-parametric exploration (MINE; Reshef ) statistics for identifying and classifying relationships between the proportion of COGs and TP. We used a maximal information coefficient (MIC) > 0.54 (uncorrected P < 0.05 and false discovery rate < 2.56e-07) to identify COGs that were significantly related to productivity (TP) in the sampled lakes. A total of 183 COGs of 3335 COGs tested were identified using these criteria, whereof 34 COGs were positively related to TP (Table S5). An inverse relationship to TP was observed for certain active transporters of phosphonates (cog3454, cog4107) and organic compounds such as amino acids (cog559, cog1147, cog4177). Homologues to other active transporters such as permeases (cog2998, cog4603, cog5265) that facilitate the transport of for example nitrate and sulphate (cog619, cog659) were negatively related with TP. The number of predicted homologues to phosphoserine phosphatase (cog560) and serine acetyltransferase (cog1045) genes involved in amino acid metabolism was negatively correlated with TP as were genes with a crucial role in carbohydrate degradation (cog153, cog1082, cog3250). Other gene products that could be useful for diagnostics of metabolic processes were carbon-monoxide dehydrogenase CoxLMS subunits (CO oxidation) that were significantly negatively related to TP. These genes are involved in the oxidation of CO to CO2 and represent an alternative or supplementary energy source that is widespread in marine bacteria (King and Weber, 2007; Brinkhoff ). CO-dehydrogenase genes were detected at higher relative abundance in three lakes with low levels of TP: Trout Bog, Damariscotta and Vättern.

Fig 3

Non-metric multidimensional scaling plot of microbial functional diversity along a productivity gradient (stress-value = 0.10). This plot is based on Horn–Morisita distances from COGs lists of 12 freshwater metagenomes. Total phosphors (TP) was mapped as en environmental variable vector onto the ordination using R function (TP) ‘envfit’. NMDS, Non-parametric-Multi-Dimensional-Scaling.

The significant relationships observed between TP and COG patterns inferred by MINE mainly provide new genome-level confirmation of earlier empirical findings of how microbial processes such as sugar, amino acid and phosphate acquisition strategies are structured along productivity gradients but also identify variations in the occurrence of response regulators that allow microbes to sense and to react to environmental stress (i.e. cog589).

Phylogenetic analyses of selected functional genes and the correspondence between functional and taxonomic composition

Phylogenetic trees were constructed for a selected number of proteins including the mmoA, nirK, pstA/B, RuBisCo and the nifH/bchL/chlL family, including Swiss-Prot references and their homologues in the metagenomes (Fig. S4). The selected genes are involved in key biogeochemical processes including methane oxidation, denitrification, phosphorus uptake, CO2 fixation, nitrogen fixation and the synthesis of photopigments. Obtained phylogenies were analysed to infer the phylogenetic structuring between and within freshwater and marine sequences using PYLOCOM (Webb ). Resulting beta nearest taxon indexes (βNTIs; see Experimental procedures) from the functional genes were compared with βNTI derived for the 16S rRNA (Table 4). These comparisons revealed that proteins, similar to the 16S rRNA, exhibit phylogenetic overdispersion between biomes when compared with random phylogenetic structures. This infers that freshwater and marine protein sequences are more different from each other than expected by chance (Webb ). This suggests that these key functional genes from marine and freshwater biomes are usually not closely related and often group into distinct marine and freshwater phylogenetic clusters, similar to what has been reported before for the 16S rRNA marker gene (Logares ).

Table 4

Results from phylogenetic analyses estimating beta-NTI within and between marine and freshwater sequences

BetaNTI		Freshwater	Marine
Freshwater	16S	−7.775
		bcn	−1.989
		nirK	0.036
		RuBisCo	15.345
		pstA/B	−3.156
		mmoA	12.287
Marine	16S	−27.123	21.883
		bcn	−4.866	1.754
		nirK	−3.736	0.278
		RuBisCo	−232.179	30.088
		pstA/B	−1.654	3.144
		mmoA	−90.028	22.273

Genes annotated as 16S rRNA and related to functional genes such as mmoA, nirK, pstA/B, RuBisCo, and the nifH/bchL/chlL family.

Values above +2 indicate phylogenetic clustering, whereas a NTI below −2 indicates overdispersion.

To determine if 16S rRNA-derived taxonomic and functional profiles among the metagenomes were coherent, a procrustes analysis was performed (Oksanen ). Our results demonstrate that the known rRNA-inferred microbial community shifts across the freshwater to marine gradient are reflected also in cohesive shifts in community-level functions observed in the metagenomes. 16S rRNA taxonomy resolved to either genus/typical freshwater taxa or phylum levels were significantly correlated with the functional data based on COG annotations (R = 0.95 and R = 0.83, respectively, P < 0.001 using procrustes analysis). When lake data were analysed separately, the procrustes analyses between 16S rRNA community composition (both phylum and genus composition) and functional COG annotations revealed similarly high coefficient values (R = 0.74 and R = 0.84, respectively, using procrustes analysis), but because these analyses included fewer samples, P-values increased dramatically (P < 0.033 and P < 0.11 respectively). This suggests that the taxonomic composition as inferred by phylogenetic markers (i.e. 16S rRNA gene) and the functional potential of communities are linked through evolutionary history. Still, it remains to be shown whether this implies that differentiation at the fine-scale population level has only minor effects on the overall gene content and potential subsequent ecosystem function, or instead is mainly determined by distribution patterns of broad taxonomic groups.

Outlook

Our metagenomic analyses of pelagic microbial communities in lakes and oceans suggest that many core functions are shared across these two biomes. Although the functional overlap is substantial, our analyses also point to some profound functional differences. Because of the rather shallow coverage of the underlying genetic diversity in the metagenomes analysed here, many genes or gene categories were not sufficiently abundant in the data set to determine with any certainty, whether or not there were significant changes in their relative abundances across the freshwater marine boundary or across the freshwater productivity gradient. This applies to genes associated with less widespread metabolic processes that may nevertheless be of critical importance to carbon and nitrogen cycling in these aquatic systems (including genes associated with N cycling, chitin degradation and ammonia oxidation). Forthcoming deeper metagenomic sequencing will likely capture trends also in these genes across environmental gradients and will help build a more comprehensive understanding of how the functional capabilities of aquatic microbial communities change along salinity and productivity gradients. Nevertheless, the present comparison of freshwater and marine metagenomes based on whole-genome shotgun sequence data did provide functional, phylogenetic and taxonomic trends across these gradients and will help us design biogeochemical experiments to test metagenome-inferred predictions such as differences in substrate preferences. Examples are the inferred prevalence towards amino acids in marine systems and difference in carbohydrate metabolism between marine and freshwaters, and the over-representation of homologues involved in the oxidation of recalcitrant organic matter in oligotrophic lakes compared with eutrophic lakes.

Experimental procedures

Sample characterization and DNA extraction

For Lakes Vättern, Ekoln and Erken, integrated water samples from the upper 2 m were collected with a rinsed 2 m Polyvinyl chloride (PVC) tube. Samples were sieved through an autoclaved 100 μm nylon mesh prior to further processing. Samples were kept dark at near in situ temperature and upon return to the laboratory, microbial cells from between 0.5 and 1 l of water were collected on replicated 0.2 μm membrane filters (Supor 200, 47 mm diameter; Gelman) by vacuum filtration followed by freezing at −80°C until further analyses. Water temperature profiles measured on site at the time of sampling verified that the sampling was limited to the upper mixed layer (epilimnion). TP and dissolved organic carbon was measured using standard methods as previously described (Eiler ). Community DNA was extracted from individual membrane filters using the FASTDNA spin kit for soil (QBiogene, Carlsbad, CA, USA) as recommended by the manufacturer. At least three membrane filters were extracted to recover sufficient DNA for 454 pyrosequencing. The amount and quality of recovered DNA was quantified by spectrophotometry at 260 and 280 nm, and aqarose gel electrophoresis revealed DNA with an average molecular weight exceeding 20 kb. All three metagenome samples were sequenced with 454 pyrosequencing with Titanium chemistry using half a chip for the Lake Erken metagenome and one quarter of a chip for Ekoln and Vättern (separated by sample specific molecular barcodes). Samples were collected from the epilimnia of Damariscotta Lake, Lake Mendota, Sparkling Lake, and Trout Bog Lake, and sequenced as described elsewhere (Martinez-Garcia ). Sequences are publicly available through the European Nucleotide Archive under project PRJEB4844.

Data mining

Metagenome data from Lake Lanier (Oh ) and two samples from Yellowstone Lake (T. McDermott, unpubl. data) were acquired from SRA (fastq-files) and analysed following the quality control and annotation procedures as described later. Fasta files for all selected marine samples were downloaded from Community cyberinfrastructure for Advanced Microbial Ecology Research and Analysis (CAMERA) (Seshadri ). Annotations were performed as described later.

Sequence annotation and functional assignment

Preprocessing was performed to bring all data sets (fasta and quality files) to the same starting point. This procedure included the following steps: length filter (length > 150) and quality filter (mean quality > 21) for lake metagenomes and just a length filter for marine metagenomes (upper length filter as listed in Table S2), clustering artificial duplicates with cd-hit-454 (Beifang ) using 97% identity threshold and 80% of the sequence in the alignment, and finally creation of consensus sequences from the clusters with cdhit-cluster consensus ignoring terminal gaps (-maxlen = 1). Quality-filtered data sets were used in all analyses. Simple six-frame translation with 50 aa length threshold was used for non-annotation-based analyses (aa usage and isoelectric point). COG annotations of the reads were extracted from the best tblastx hit against STRING (Franceschini ), and rRNAs were identified using hmm rRNA to obtain annotations. An E-value threshold of 1e-10 was applied.

We also performed a second preprocessing and annotation procedure, and subsequent statistical analyses in which results supported the main findings presented earlier. In short, a more stringent quality filtering was performed with cutting reads when quality scores dropped below 21 and using a length cut-off of 150 bp. Clustering artificial duplicates was performed as described earlier. The quality-filtered data sets were then submitted to CAMERA using Rapid Analysis of Multiple Metagenomes with a Clustering and Annotation Pipeline (RAMMCAP) (Seshadri ; Weizhong, 2009) with the following parameters: six-frame translation, hmm rRNA and annotation (no clustering), which masks tRNAs and rRNAs before calling ORFs. Subsequently, Reversed Position Specific-Basic Local Alignment Search Tool (RPS-BLAST) was performed against COG (Tatusov ). An E-value threshold of 1e-10 was applied. Fasta files for the marine data sets downloaded from CAMERA were used without any quality filtering, except cd-hit-454 for artificial duplicate removal.

Statistical analyses

To ensure robustness of the statistical tests to outliers, we have compared the results using three types of data sets: all COGs from 12 lake and 13 marine samples; all COGs from 8 lake and 12 marine samples; and only bacterial COGs from 12 lake and 13 marine samples. The smaller number of samples for the second set resulted from excluding samples with the worst quality processing results. Bacterial COGs are used to address the issues of varied eukaryotic content between the samples. The abundance of individual reads matching a particular COG were normalized against the average abundance of 35 likely essential and single copy COGs (Ciccarelli ; Raes ) and used to generate a metabolic profile of the metagenome. This provides a proxy for the number of genomes harbouring a specific COG in the community. Core-gene normalized profiles were then used in statistical analyses such as metaMDS, PERMANOVA and procrustes test with Horn-Morisita distance measure using the functions in the ‘ecodist’ and ‘vegan’ libraries in the R-package (http://www.r-project.org; Goslee and Urban, 2007; Oksanen ). PERMANOVA (Anderson, 2001) was used to determine significant differences between freshwater and marine functional beta-diversity, and procrustes analysis was used to determine correspondence between taxonomic and functional composition. To fit TP as an environmental vector onto the ordination, we used the function ‘envfit’. The fitted vector is an arrow that points to the direction of its most rapid change in the ordination space (direction of the gradient), and its length is proportional to the correlation between community composition and TP. Prior to applying the Wilcoxon test, COGs were resampled and then normalized against the single-copy core COGs to identify over- and under-represented COGs in freshwater compared with marine metagenomes. False discovery rate (q-value) was estimated after Storey (2002). MINE (Reshef ) was used with default settings for identifying and classifying relationships between the resampled and normalized COG abundances, and TP that was used as a proxy of lake productivity (Schindler, 1978). Relationships were defined as significant when the MIC was > 0.54 with a P-value < 0.01 and a false discovery rate < 2.56e-07.

Taxonomic assignments

Protein-based taxonomic assignments for domain and phylum were extracted from the best hit to the STRING database (E-value threshold 10−10). In addition, SSU rRNAs were extracted by hmm rRNA. The Bayesian classifier (Wang ) (using bootstrap cut-off > 60) was used to annotate 16S rRNA genes against a custom curated freshwater-specific database (Newton ) and the SILVA database using taxonomy after SILVA (Quast ). The number of reads annotated to the different bacterial phyla and bacterial ‘genera’ were extracted and ordinated using R (Oksanen ). A procrustes test was used to compare the functional annotations with the taxonomic annotations at the genus level.

Phylogenetic analyses

Reference (master) sequences for mmoA, nirK, pstA/B, RuBisCo and the nifH/bchL/chlL family were obtained from Swiss-Prot. After six-frame-shift translation of the sequences (using a minimum length of 50 aa), homologous ORFs in the 29 metagenomes were identified based on blastp searches using an E-value threshold of 1e-10 and per cent identity of 40%. Alignments of master sequences were obtained for each of the five genes using Multiple sequence comparison by log-expectation (MUSCLE) (default settings; Edgar, 2004). Preliminary multiple sequence alignment were obtained for the metagenomic ORFs by MUSCLE using settings -maxiters 1 and -diags to increase speed. These ‘slave’ alignments were then aligned against the master alignment with muscle using function ‘-profile’. Bootstrapped Random Axelerated Maximum Likelihood (RAxML) trees (Stamatakis ) were computed based on trimmed master alignments using standard model JTTF and default convergence criteria. Trees and alignments were imported into ARB (Ludwig ), and the quick parsimony option was used to add the aligned metagenomic ORFs to the RAxML master trees. For 16S rRNA genes, the procedure outlined in Peura and colleagues (2012) was used to insert metagenomic 16S rRNA homologues into the SILVA106 reference tree. Phylogenetic trees were visualized using iTOL (Letunic and Bork, 2011) and analysed using PHYLOCOM (Webb ). The phylocom function ‘comdistnt’ was used to infer if freshwater and marine sequences were phylogenetically distinct by estimating the βNTI. Here, we used both the marine and freshwater biomes as separate groups. Mean nearest taxon distance (MNTD) was estimated for within each biome and between biomes. To weigh phylogenetic distances by taxa abundances, the average distance among random individuals drawn from each of the two biomes was calculated. The NTI was quantified by the number of standard deviations that the observed MNTD is from the mean of the null distribution (999 randomizations; MNTDnull). MNTDnull is found by randomizing OTUs across the phylogeny and recalculating MNTD 999 times.

NTI = −1*(MNTDobs-mean MNTDnull/sdMNTDnull) (Webb ).

For a single community, NTI greater than +2 indicates that coexisting taxa are more closely related than expected by chance (phylogenetic clustering). NTI less than −2 indicates coexisting taxa are more distantly related than expected by chance (phylogenetic overdispersion). βNTI is the between-group analogue of NTI (Fine and Kembel, 2011; Webb ).