Diversity Manipulation of Psychrophilic Bacterial Consortia for Improved Biological Treatment of Medium-Strength Wastewater at Low Temperature

Psychrophilic bacteria are valuable biocatalysts to develop robust bioaugmentation formulations for enhanced wastewater treatment at low temperatures or fluctuating temperature conditions. Here, using different biodiversity indices [based on species richness (SR), phylogenetic diversity (PD) and functional diversity (FD)], we studied the effects of microbial diversity of artificial bacterial consortia on the biomass gross yields (measured through OD600) and removal efficiency of soluble chemical oxygen demand (mg sCOD removed/mg sCOD introduced) in synthetic, medium-strength wastewater. We built artificial consortia out of one to six bacterial strains isolated at 4°C through combinatorial biodiversity experiments. Increasing species richness resulted in improved sCOD removal efficiency (i.e., 0.266 ± 0.146, 0.542 ± 0.155, 0.742 ± 0.136, 0.822 ± 0.019 for mono-, tri-, penta-and hexacultures, respectively) and higher biomass gross yields (i.e., 0.065 ± 0.052, 0.132 ± 0.046, 0.173 ± 0.049, 0.216 ± 0.019 for mono-, tri-, penta,- and hexacultures, respectively). This positive relationship between biodiversity, sCOD removal and biomass gross yield was also observed when considering metabolic profiling (functional diversity) or evolutionary relationships (phylogenetic diversity). The positive effect of biodiversity on sCOD removal efficiency could be attributed to the selection of a particular, best-performing species (i.e., Pedobacter sp.) as well as complementary use of carbon resources among consortia members (i.e., complementarity effects). Among the biodiversity indices, PD diversity metrics explained higher variation in sCOD removal than SR and FD diversity metrics. For a more effective bioaugmentation, our results stress the importance of using phylogenetically diverse consortia, with an increased degradation ability, instead of single pure cultures. Moreover, PD could be used as an assembly rule to guide the composition of mixed cultures for wastewater bioaugmentation under psychrophilic conditions.

Psychrophilic bacteria are valuable biocatalysts to develop robust bioaugmentation formulations for enhanced wastewater treatment at low temperatures or fluctuating temperature conditions. Here, using different biodiversity indices [based on species richness (SR), phylogenetic diversity (PD) and functional diversity (FD)], we studied the effects of microbial diversity of artificial bacterial consortia on the biomass gross yields (measured through OD 600 ) and removal efficiency of soluble chemical oxygen demand (mg sCOD removed/mg sCOD introduced) in synthetic, medium-strength wastewater. We built artificial consortia out of one to six bacterial strains isolated at 4 • C through combinatorial biodiversity experiments. Increasing species richness resulted in improved sCOD removal efficiency (i.e., 0.266 ± 0.146, 0.542 ± 0.155, 0.742 ± 0.136, 0.822 ± 0.019 for mono-, tri-, penta-and hexacultures, respectively) and higher biomass gross yields (i.e., 0.065 ± 0.052, 0.132 ± 0.046, 0.173 ± 0.049, 0.216 ± 0.019 for mono-, tri-, penta,-and hexacultures, respectively). This positive relationship between biodiversity, sCOD removal and biomass gross yield was also observed when considering metabolic profiling (functional diversity) or evolutionary relationships (phylogenetic diversity). The positive effect of biodiversity on sCOD removal efficiency could be attributed to the selection of a particular, best-performing species (i.e., Pedobacter sp.) as well as complementary use of carbon resources among consortia members (i.e., complementarity effects). Among the biodiversity indices, PD diversity metrics explained higher variation in sCOD removal than SR and FD diversity metrics. For a more effective bioaugmentation, our results stress the importance of using phylogenetically diverse consortia, with an increased degradation ability, instead of single pure cultures. Moreover, PD could be used as an assembly rule to guide the composition of mixed cultures for wastewater bioaugmentation under psychrophilic conditions.

INTRODUCTION
Wastewater temperature is a key parameter influencing microbial functions in aerobic wastewater treatment processes (Zhou et al., 2018). Due to climate seasonal changes or differences in geographic area, wastewater temperature can decrease below 10 • C, resulting in a reduction of biological activity and reaction rates, as described by the van't Hoff-Arrhenius equation, and in a lower chemical oxygen demand (COD) removal efficiency (Tchobanoglous et al., 2014). Bioaugmentation using psychrophiles has been proposed as a strategy to improve the performance of specific clean-up operations of municipal/domestic wastewater at low temperatures (e.g., carbon removal through oxidation and biomass growth) and compensate for impaired activities of mesophilic bacteria (Zhou et al., 2018).
Psychrophiles are cold-adapted microorganisms living at temperatures close to the freezing point of water, with an optimum growth temperature below 15 • C and an upper cardinal temperature of about 20 • C (Morita, 1975;Feller and Gerday, 2003). They have the potential to exhibit high metabolic activities at low and moderate temperatures using cold-adaptive traits that compensate for the adverse impact of low temperatures on biochemical reaction rates (Feller and Gerday, 2003;Margesin et al., 2005). The use of psychrophilic strains in bioaugmentation processes was already investigated under cold conditions. For example, a cold-adapted Arthrobacter psychrolactophilus was able to grow in a synthetic wastewater at 10 • C, inducing a complete clarification of the turbid medium with an efficient hydrolysis of proteins, starch and lipids (Gratia et al., 2009). In another study, psychrophilic bacteria and yeasts fully degraded phenol (a typical aromatic toxic pollutant frequently detected in wastewater) at 10 • C under fed-batch cultivation (Margesin et al., 2005). The selection of robust and functionally-active bioaugmentation cultures is of primary importance to optimize the invasion success by mastering the propagule pressure and the establishment of exogenous microorganisms in the invaded ecosystems (El Fantroussi and Agathos, 2005). From product development to final application, bioaugmentation formulations typically progress through five life stages, including (i) capture, (ii) production, (iii) establishment, (iv) function, and (v) downstream impacts (Kaminsky et al., 2019). Bioaugmentation attempts have frequently failed due to poor survival of the inoculated strains. The antagonistic effects of the invaded ecosystem conditions on the active inoculated microbial invaders (i.e., microbiostasis Ho and Ko, 1982) is attributed to both biotic and abiotic factors (van Veen et al., 1997;El Fantroussi and Agathos, 2005). Biotic factors include strain selection criteria based on unique features that confer ecological advantages in the ecosystem, magnitude of propagule pressure, competition between inoculant and indigenous populations or predation by protozoa and bacteriophages. Abiotic factors comprise temperature, pH, substrate availability, or presence of toxic compounds. Both the traits of the invader(s) that allow successful invasion (i.e., invader-centric research Ma et al., 2015;Kinnunen et al., 2016) and the properties of the resident community that determine its susceptibility to invasion (i.e., resident community-centric approach Kinnunen et al., 2016Kinnunen et al., , 2018Mallon et al., 2018) have recently been investigated.
From the invader-centric perspective, the formulation of microbial inoculants based on eco-physiological attributes and the applications of polycultures instead of single strains in isolation (pure cultures) represent two different strategies to improve bioaugmentation, increasing the survival and functional robustness of the introduced strains.
In particular, the use of consortia can ensure higher removal efficiency of contaminants and more robust processes thanks to their higher functional redundancy, diversity and stability (Stenuit and Agathos, 2015;Giri et al., 2020). Indeed, more diverse communities can gain greater benefits from niche opportunities and assimilate a greater proportion of the available resources (Cardinale, 2011). They also show a greater ability to withstand specific disturbances (resistance) or return to their undisturbed structural and functional baseline after perturbation (engineering resilience; Stenuit and Agathos, 2015).
In this study, the effects of increasing diversity of synthetic, psychrophilic consortia was investigated on the removal of soluble chemical oxygen demand (sCOD) and biomass gross yields in synthetic, medium-strength wastewater at 4 • C. The designed synthetic wastewater contained diverse carbon sources available as growth substrates to increase resource heterogeneity for heterotrophic microorganisms (i.e., high niche dimensionality). After isolating six bacterial strains for their ability to remove sCOD from synthetic wastewater at 4 • C, we assembled synthetic psychrophilic consortia with increasing richness level, using a combinatorial experiment. The removal of sCOD and the biomass gross yields were analyzed over time and used as functional proxies to assess the coupling of biodiversity and microbial ecosystem functioning. We hypothesized that increased bacterial diversity can enhance sCOD removal efficiency and wastewater treatment through complementarity between species as the underlying mechanism for function optimization at high niche dimensionality. In addition, different biodiversity metrics (i.e., species richness (SR), phylogenetic metrics [Faith's phylogenetic diversity (FPD), mean pairwise distance (MPD)] and functional diversity metrics [dendrogrambased functional diversity (dFD), functional dissimilarity (FDis)] were tested for their relevance to the selection and assembly rules of efficient consortia for bioaugmentation.

Bacteria Isolation and Identification
The six psychrophilic strains used in this study were isolated from laboratory cold-room facilities and selected for their ability to grow on SW at 4 • C. The isolates were identified using the sequencing of almost full-length 16S rRNA genes. The primer pair used to amplify the 16S rRNA gene included the universal bacterial primer pair 27f (5 ′ -AGAGTTTGATCMTGGCTCAG-3 ′ ; Lane, 1991) and 1492r (5 ′ -GGTTACCTTGTTACGACTT-3 ′ ; Reysenbach et al., 1992). Herculase II Fusion DNA Polymerase (Agilent Technologies) was used following the manufacturer's instructions. The PCR reactions contained 10 µL 5 × Herculase II reaction buffer, 250 µM each dNTP, 250 nM forward and reverse primer (1.25 µL of 10 µM each primer), 1 µL Herculase II Fusion polymerase, 50-200 ng genomic DNA and water to a final volume of 50 µL.
The PCR thermal cycling scheme was: 95 • C for 2 min, 30 cycles (95 • C for 20 s, 48 • C for 20 s, 72 • C for 45 s) and a final extension at 72 • C for 3 min. The PCR products of the expected size were excised from the agarose gel (around 1,500 bp) and purified using the QIAquick Gel Extraction Kit (Qiagen). All purified PCR products were sequenced at Macrogen Europe (Amsterdam, The Netherlands). Taxonomic identification was carried out using the Basic Local Alignment Search Tool (nucleotide BLAST using the Nucleotide collection (nt) consisting of GenBank, EMBL, DDBJ, PDB and RefSeq sequences; Altschul et al., 1997). The genus name with the highest maximal identity percentage was retained. The bacterial species belonged to the following taxa: Rhodococcus sp., Pedobacter sp., Janthinobacterium sp., Brevundimonas sp., Pseudomonas sp., and Arthrobacter sp. The 16S sequences were deposited in GenBank under the accession numbers MN722456-MN722461 (Table S2).

Combinatorial Biodiversity Experiment
To evaluate the individual and combined effects of the six different strains (number of factors k = 6), bacterial species were assembled with a two-level, ½ fractional factorial design with 32 runs (2 6−1 ) ( Table S3). The ½ parameter states that only a fraction (50%) of the runs defined by the full factorial design is considered in this experiment (the half fraction in the design is equal to 2 k−1 ). The presence (+1) or absence (−1) of the strains were used as the levels of the factors and biodegradation activity (COD removal efficiency) and biomass gross yields as the responses. The chosen fractional factorial design resolves all two-factor interactions, with the selection of the negative sign associated with the generating rule. In addition, abiotic controls and hexacultures were included in the design to gain more information without excessively increasing the number of laboratory incubations. Each community was present in triplicate, resulting in the monitoring of 114 ecosystems, including the abiotic controls. A single colony of each bacterial culture was picked and grown in R2A broth (R2B). Cells were harvested at mid-exponential phase, washed (3×) in sterile phosphate buffered saline (PBS) (pH 7.2) and adjusted to a concentration of 8 × 10 8 cells mL −1 , using a particle counter (Multisizer 3 Coulter Counter, Beckman Coulter, CA, USA). Bacterial cultures were left for maximum 6 h at room temperature before assembling the communities. Mixed synthetic consortia at a concentration of 8 × 10 8 cells mL −1 were assembled in 1.5 mL Eppendorf tubes, using an equal amount of individual cell suspensions. The initial bacterial concentration of each pure and mixed culture was fixed at 2.0 × 10 7 total cells mL −1 . To achieve this, 45 µL from each assemblage were inoculated in 1,755 µL of SW in 2 mL microplates (MASTERBLOCK R 96-Well Deep Well Microplates, Greiner Bio-One). To avoid cross-contamination between the wells during the preparation and the sampling, separated empty wells were used (Figure S1). A gas permeable seal (Breathe-Easy R Gas Permeable Sealing Membrane for Microtiter Plates, Diversified Biotech, Dedham, MA) was then placed upon every microplate to enhance gas transfer, keep sterility and avoid contamination in the plate. This seal was replaced after each sampling. Additionally, three negative controls were present in each plate to verify sterile conditions. The total number of microbial ecosystems was 114. The plates were then incubated at 4 • C at 120 rpm for 8 days (192 h). Samples were taken at regular intervals (i.e., every 2 days) to measure the biomass growth. To do so, 150 µL was transferred to microtiter plates and growth was assessed through optical density measurements at 600 nm [OD 600 , (PowerWave HT Microplate Spectrophotometer, BioTek, Winooski, VT, USA)]. At the end of the experiments, 1 mL of culture was filtered and used to quantify sCOD using Spectroquant R test kits (Merck Millipore, Germany).

Single Bacterial Growth Kinetics
Cell growth kinetics of the six bacterial isolates were examined in a batch reactor configuration using SW as the growth medium.
Cultures were incubated at 4 • C and agitated at 180 rpm on a rotary shaker. Before starting the kinetic studies, pre-cultures of each individual strain were carried out at 4 • C in R2B. Cells were harvested at mid-exponential phase by centrifugation (10,000 rpm for 10 min), washed three times in sterile PBS (pH 7.2) and then re-suspended in PBS to obtain a concentrated cell suspension. The strains were individually inoculated in 50 mL SW to give an initial OD 600 of 0.025. Three biotic replicates were performed for each culture in 250 mL Erlenmeyer flasks. Samples were taken at regular intervals over a period of 8 days and growth was monitored using OD 600 measurements, plate counting and particle counter. For plate counting, samples were serially diluted in sterile PBS buffer, and plated on R2A agar. The plates were incubated at 15 • C for 3-5 days before counting. For the bacteria enumeration through the particle counter, bacterial samples were transferred into the ISOTONE R II diluent (Beckman Coulter, CA, USA) and quantified with a 20 µm aperture tube. Particle counts for the blank sample, consisting only in ISOTONE R II diluent, were subtracted from all bacterial samples.

Biodegradation Activity and sCOD Removal Efficiency
The biodegradation of carbonaceous constituents was quantified through measurement of sCOD.
The samples were filtered using 0.2 µm pore-diameter syringe filters (Minisart R , Sartorius AG, Germany) and analyzed with the Spectroquant R test kits (Merck Millipore, Germany). The ranges of measurement used were 10-150 and 25-1,500 mg O 2 L −1 , according to the manufacturer's instructions.

Computation of Diversity Metrics for Bacterial Assemblages
Faith's phylogenetic diversity (FPD) was computed as the sum of all branch lengths of a phylogenetic tree connecting all species in a local community (Faith, 1992;Vellend et al., 2011). The rooted phylogenetic tree was constructed based on almost fulllength 16S rRNA gene sequences, using the Unweighted Pair Group Method with Arithmetic Mean (UPGMA) in Geneious 11.1.4 (https://www.geneious.com). The mean pairwise distance (MPD) was also calculated as the average phylogenetic distance between each pair of species in a community (Webb et al., 2002;Vellend et al., 2011). The FPD and MPD were computed with the picante package in Rstudio (Version 1.0.136; Kembel et al., 2010). Functional diversity metrics based on metabolic profiles were assessed using the Biolog Ecoplate (BIOLOG, Hayward, CA, USA). The Ecoplate contains 31 carbon sources (plus one blank) belonging to different chemical families, repeated three times. Every well of the plate also contains a fluorogenic tetrazolium dye (5 cyano-2,3 ditolyl tetrazolium chloride), which is reduced to a violet-fluorescent formazan molecule when the bacterial cells oxidize the carbon source (Gravel et al., 2011). Briefly, overnight cultures were washed twice in PBS, adjusted to an OD 600 of 0.05 and added to the plates (150 µL per well). The plates were incubated at 4 • C and the color development on each substrate, indicative of bacterial activity, was measured at 590 nm (Jousset et al., 2011). The average OD 590 values for each substrate were first subtracted from the blank (AU substrate ). Each AU substrate value was first normalized by the number of electron equivalents (eeq) of the corresponding substrate in each well (AU substrate /eeq substrate = AUeeq substrate ) and then divided by the average value of the whole plate (AUeeq substrate / AUeeq plate ). Two methods were used to measure functional diversity (Salles et al., 2012): the first approach was based on distances (Functional Dissimilarity, FDis Heemsbergen et al., 2004;Jousset et al., 2011) and the second based on a dendrogram [dFD (Petchey and Gaston, 2002)]. FDis is the mean functional distance between each pair of species in trait space; dFD corresponds to the sum of branch lengths of a functional dendrogram connecting species of a community together (Petchey and Gaston, 2006). FDis was computed by using a Euclidean distance matrix to sum the distances for all pairs of species in a community. The same distance matrix was used to perform a hierarchical clustering and build the functional dendrogram to compute dFD. The functions "dist()" and "hclust()" in R were used to generate the two metrics. In addition a physiological profile heatmap was drawn with the "pheatmap()" function ( Figure S2).

Data Analysis
Two complementary metrics using the log response ratio were computed to quantify potential effects of biodiversity (i.e., complementarity and selection effects) that governed the diversity-sCOD degradation and diversity-biomass gross yield couplings.
LRm = ln sCOD-RE polyculture sCOD-RE best-performing monoculture (4) LRm = ln OD 600,polyculture OD 600,best-performing monoculture (5) LR m compares the value of the response variable (i.e., sCOD removal efficiency or OD 600 ) in a polyculture with the mean value of the response variables of the constituent species grown as monocultures (Cardinale et al., 2006). Positive values of LR m indicate the presence of non-transgressive overyielding, i.e., the community performs better than would be expected from the average performance of the single member species. LRm compares the value of the response variable in a polyculture with the value of the constituent species showing the highest value as monoculture (Cardinale et al., 2006). Measuring values of LRm > 0 corresponds to testing for the presence of transgressive overyielding, i.e., the community outperforms the best-performing monoculture of the component species. The initial OD 600 in all cultures was adjusted to 2.0 × 10 7 total cells mL −1 and the initial value of sCOD in SW was 540 ± 21 mg sCOD L −1 . Continuous lines denote simple linear regression significant at p < 0.05.
All the statistical analyses were performed using GraphPad Prism software (version 8.3.1, San Diego, CA). Simple linear regressions were used to test for the relationships between community functioning and the different diversity metrics. Data were tested for normality of residuals using the D'Agostino-Pearson test. Differences in response variables, LR m and LRm between species richness levels were assessed by means of ordinary One-Way ANOVA or Welch's ANOVA and Kruskal-Wallis test when the ANOVA criteria were not met. For LR m and LRm, post hoc multiple comparison tests were carried out to analyze the statistically significant differences of log response ratios for specific pairs of richness levels.
A screening for the main effects (i.e., individual effect of each factor on a dependent variable) of each strain on the two response variables was computed in JMP R software (version 12.2.0). After fitting a multiple linear regression model including the six strains as predictors, the t-ratio, which is defined as the main effect estimate divided by its standard error, is analyzed to test the null hypothesis. In particular, the t-ratio indicated whether a factor (strain) had a significant influence on the response variables and to which extent. For all statistical analyses, a p < 0.05 was used to determine statistical significance.
Using simple linear regression, a significant positive relationship between biodiversity dimensions (i.e., SR, FPD, and dFD) and both biomass growth and sCOD removal efficiency was observed [OD 600 : Lower percentages of variance explained by the model were observed for the independent variables SR, FPD and dFD when considering the OD 600 as the response variable of community functioning (38, 35, and 36%, respectively) compared to sCOD conversion efficiency (53, 62, and 51%, respectively; Figures 1A-F). FPD was the best predictor of sCOD removal with a percentage of variance of 62% explained by the model. Because the models using sCOD conversion efficiency as system function explained higher variances, we focused on this response variable to analyze in more detail the different effects of the independent variables (i.e., diversity metrics).
The strong correlation between SR-FPD and SR-dFD ( Figure S4) makes it difficult to properly assess the effect of FPD or dFD per se on the sCOD removal efficiency. To disentangle the single effects of FPD and dFD, the data set was first separated into the different SR levels (i.e., SR = 3 and SR = 5) and the individual effects of FPD and dFD were quantified within the two SR levels. FPD had a significant and positive relationship with sCOD removal for both SR levels [SR = 3, F FPD(1,18) = 14.30, p = 0.0014; SR = 5, F FPD(1,4) = 7.97, p = 0.0477; Figures 2A,B]. Conversely, dFD was negatively related to sCOD removal within both SR levels, although neither of these correlations were statistically significant [SR = 3, F dFD(1,18) = 0.1931, p = 0.6656; SR = 5, F dFD(1,4) = 0.6473, p = 0.4662; Figures 2C,D].
To reduce the collinearity with SR, we also used two other metrics that are less correlated to species richness (Swenson, 2014): mean pairwise distance (MPD, from phylogenetic analysis) and functional dissimilarity (FDis, based on metabolic profile). The lack of correlation between MPD-SR and FDis-SR was also confirmed by directly analyzing the relationships in our experiment through Spearman's correlation (r ∼ 0; Figure S4).

Effect of Species Identities
To analyze the effect of species identity on the coupling of biodiversity and functioning, a screening for the main effects (i.e., individual effects of each factor) was realized in JMP R . Each species can contribute differently to the targeted function and this is closely related to the selection effect (i.e., higher probability of including particularly productive species that disproportionately contribute to the targeted function in more diverse communities; Loreau, 1998). The plot of each t-ratio is shown in Figure 3A. A different contribution of species identities was observed depending on which response variable was considered as a measure of system functioning. Using the values of OD 600 after 144 h of incubation as a target function, the strain Pseudomonas sp. showed the highest effect, followed by Rhodococcus sp. The other strains (i.e., Pedobacter sp., Brevundimonas sp., Arthrobacter sp., and Janthinobacterium sp.) did not significantly contribute to biomass gross yield.
On the contrary, concerning the sCOD removal ability [8 days (192 h) after inoculating the synthetic wastewater], a different effect of species identity was observed: Pedobacter sp. had the highest effect on sCOD removal, followed by Rhodococcus, Pseudomonas, and Brevundimonas sp. (Figure 3A).
The activity of Arthrobacter sp., considering both OD 600 and sCOD removal, did not influence the aggregate community property and the whole behavior of the system.
The importance of Pedobacter sp. for sCOD degradation is also reflected by the effect of species composition on degradation activity ( Figure 3B): the highest degradation activity was observed in communities that contained Pedobacter sp.
To explain the differences in species contribution depending on the chosen response variable, we carried out single bacterial growth kinetics on the synthetic wastewater at 4 • C, measuring the OD 600 and cell density (cells mL −1 ). In most of the cases, time courses of cell growth measured through plate counting did not correspond to those obtained with OD 600 measurements ( Figure S6). In other words, when a stationary phase is detected based on the OD 600 curve, an exponential growth phase is observed at the same incubation time for the growth curve based on bacterial cell counting on agar plates [i.e., Brevundimonas sp. (Figure S6C), Pedobacter sp. (Figure S6E)]. Moreover, to a nearly similar value of OD 600 for two individual strains, different cell concentrations were measured indicating different conversion factors of one OD 600 unit into cell concentrations (cells mL −1 ) for each individual strain due to different cell morphologies and dynamics of the total cell volume in a sample ( Figure S7). This behavior could be explained by taking into account the time profiles of the mean cell volume (Figure S7), computed from the ratio of the volume of cells per mL and the number of cells per mL obtained from the Coulter counter. The mean cell volume of Rhodococcus sp. was almost 10-fold higher than the cell volume of the other strains, whereas the mean cell volume of Brevundimonas sp. was the smallest. In the specific case of Pedobacter sp. and Brevundimonas sp., the quasi constant values of OD 600 (Figures S6E,C) were due to the decrease in volume starting from the mid-exponential phase and a concomitant increase in cell concentrations. The same trend was observed in general for all bacteria.

Analysis of Non-transgressive and Transgressive Overyielding (Complementarity Effect)
We computed both non-transgressive (LR m ) and transgressive overyielding (LRm) to identify potential complementarity effects that could underlie the observed positive diversityfunction couplings (Figures 1A-F). For both biomass growth and sCOD removal activities, we obtained positive values of non-transgressive overyielding for each species richness level (Figures 4A,B). Moreover, there was a significant increase in the mean values between the tri-, penta-, and hexacultures (Tukey's test, Figure 4A; Dunn's test, Figure 4B). However, polycultures can exhibit non-transgressive overyielding through complementarity or selection effect. In order to assess the presence of complementarity effects with more accuracy, we computed the transgressive overyielding. Indeed, polycultures should exhibit transgressive overyielding when both complementarity and selection underpin a positive effect of diversity on ecosystem functioning. Considering the biomass growth, the LRm values did not significantly differ from 0 for tri-and pentaculture and the differences between richness levels were not significant (Dunn's test, Figure 4C). Conversely, LRm values associated with sCOD removal efficiency were significantly positive for all richness levels and differed significantly comparing the tri-, penta-, and hexacultures (Dunnett's T3 test, Figure 4D).

DISCUSSION
To study the role of biodiversity on synthetic medium-strength wastewater (SW) treatment under psychrophilic conditions, we constructed artificial microbial communities containing up to six different bacterial isolates. We carried out combinatorial biodiversity experiments using a 2 6−1 fractional factorial design to correlate microbial ecosystem functioning and biodiversity. The investigated components of biodiversity included species richness and identities, phylogenetic diversity metrics (MPD and FPD indices based on 16S rRNA gene analysis) and functional diversity metrics (FDA and dFD indices obtained from metabolic profiles on Biolog Ecoplate). Artificially constructed ecosystems in sterile laboratory devices were used to study the activities of defined mixed cultures growing on synthetic wastewater under psychrophilic conditions. The response variables comprised sCOD and the optical density at 600 nm (OD 600 ) to assess the effects of biodiversity on sCOD conversion efficiency and biomass gross yields. The statistically-significant positive relationships between biodiversity and both sCOD conversion efficiency and biomass gross yields (Figures 1A-F) supported our hypothesis that higher bacterial diversity would result in an improved wastewater treatment, resource utilization efficiency and biomass productivity (Venail and Vives, 2013).
No saturation in functioning was observed with increasing biodiversity (Figures 1A-F), but this trend was probably due to the low diversity and low functional redundancy in the constructed consortia compared to natural ecosystems characterized by a much higher biodiversity difficult to reproduce in laboratory experiments with tractable model communities (Langenheder et al., 2012).
Another limitation of our study is the difficulty to measure individual population density (and so the evolving diversity during the incubation and the final diversity at the end of the experiment in terms of species evenness). Indeed, the relative abundance of populations in the constructed consortia (species evenness) will likely change in the course of incubation. Therefore, our results refer sensu stricto to the coupling of initial bacterial diversity (initial richness and evenness) with sCOD conversion efficiency and bacterial growth (Langenheder et al., 2012).
The effects of community diversity indices (SR, FPD, dFD, MPD, and FDis) were analyzed on the two targeted functions. SR, FPD, and dFD explained lower percentages of variance when considering biomass gross yield instead of sCOD removal efficiency. Focusing on the latter function, FPD was the best predictor of community functioning (i.e., it explained the highest percentage of variance (62%) in sCOD removal efficiency). Due to the high correlation between SR and FPD or dFD (Figure S4), we separated the data set into the different SR levels (i.e., SR = 3 and 5) to properly investigate the single effect of FPD and dFD on the sCOD removal. Only FPD was positively correlated with degradation ability within individual SR levels (Figure 2). Therefore, this suggests that not only species richness is important in explaining higher system functioning but that, at the same SR level, other diversity metrics, such as FPD, might be more valuable determinants of system functioning. Moreover, a significant variation in sCOD removal efficiency was also explained by phylogenetic diversity measured as MPD, which is less dependent on SR than FPD (Swenson, 2014;Venail et al., 2015) (Figures S4A,C, S5A). Our results supported the rationale that species sharing distant common ancestors are more likely to be more functionally unique, increasing the chance for complementarity and higher overall ecosystem functioning. Previous studies already showed a positive effect of biodiversity driven by phylogenetic diversity and its higher predictive power in explaining ecosystem functioning (Jousset et al., 2011;Venail and Vives, 2013;Galand et al., 2015).
The two indices of functional diversity (dFD and FDis) explained a lower ( Figure 1F, Figure S5B) or no variation in sCOD removal (Figures 2C,D). This is probably because the computation of functional diversity metrics is highly dependent on the methodology used to obtain physiological and phenotypic profiles. In this study, we used the Biolog Ecoplates containing 31 different carbon sources to assess the metabolic diversity of individual bacterial isolates. As suggested previously (Wei et al., 2015;Yang et al., 2017), the design of tailor-made resource use profiling plates can be more representative of the ecological niches of the investigated ecosystem by targeting specific relevant resources and phenotypic traits.
Due to the difficulties of measuring individual population density (due to low volume of the system and low biomass), we could not assess the contributions of individual bacterial species to system functioning through the additive partitioning of biodiversity effects (i.e., complementarity and selection effects; Loreau and Hector, 2001). However, we tried to provide a rough quantitative measurement on the relative importance of selection and complementarity effects in our consortia. The analysis of the effect of species identity on the responses allowed us to determine which bacterial species contributed the most to the targeted function and it represented a rough assessment for the presence of a selection effect. Here, the experimental design of combinatorial biodiversity experiments enabled us to study the main effects of each factor used in the fractional factorial design. Species contributed differently to the functionality of the microbial ecosystem depending on which response was considered as a measure of system functioning (Figure 3). Single bacterial growth kinetics on the synthetic wastewater at 4 • C revealed that the biomass gross yield response was biased by OD 600 measurements due to differences in cell volumes between strains. Moreover, the increase in cell concentrations during the kinetics was characterized by a simultaneous decrease in cell volumes that made the measured OD 600 values less representative of the actual growth. These results demonstrate that OD 600 , often used as a proxy for biomass productivity, might lead to a biased interpretation of data due to morphological differences between bacteria and temporalfluctuations in cell sizes.
Here, we show the importance of carefully choosing the target function to link a specific biodiversity dimension and ecosystem functioning, and the necessity to properly select the variable to measure the community function (e.g., total community biomass). The OD 600 has often been used as a common method for estimating the concentration of bacterial cells in microbial BEF experiments (Eisenhauer et al., 2013;Awasthi et al., 2014), but the results might be biased due to different cell volume among the different strains or fluctuating cell volumes in the course of time. To analyze the presence of potential complementarity effects that could underpin, along with selection effects, the observed positive diversity-function couplings, we measured the transgressive overyielding (LRm) for each polyculture. The transgressive overyielding is a widespread measure used in ecology to quantify the co-contribution of both complementarity and selection effect on the increased system functioning. Polycultures can exhibit non-transgressive overyielding through complementarity or selection effect (Schwinning et al., 2014). When both complementarity and selection underlie a positive effect of diversity on ecosystem functioning, polycultures should exhibit transgressive overyielding (Schwinning et al., 2014;Vasseur and Messinger, 2015). The occurrence of transgressive overyielding in all richness levels was observed only for the removal of sCOD (Figure 4D), suggesting that not only selection effect (due to the presence of Pedobacter sp., the best-performing strain) but also complementarity effects between species were common in our experiment. This is in accordance with previous studies where the positive relationship between biodiversity and community productivity was driven by the complementarity effect of diversity (Cardinale, 2011;Venail and Vives, 2013).

CONCLUSIONS
In conclusion, microbial ecosystem functioning and robustness can be optimized for biotechnological applications using the rational design of more diverse synthetic consortia, instead of pure cultures. Positive influences of microbial biodiversity on psychrophilic wastewater treatment could be attributed to both complementarity and selection effects. Despite the low diversity levels used in our experiment, phylogenetic diversity turned to be the major determinant of community performance and a valuable predictor of complementary resource utilization between species. Therefore, it could be used as an assembly rule to formulate consortia able to efficiently degrade the wastewater organic load. Further studies should now investigate the use of the best-performing psychrophilic consortium as a bioaugmentation formulation to treat real domestic wastewater in a reactor configuration that simulates continuous-flow wastewater treatment systems.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and