extfieldsofstudy:
Biology
Computer Science

text:
Chromosome-level quality scaffolding of brown algal genomes using InstaGRAAL, a proximity ligation-based scaffolder

Hi-C has become a popular technique in recent genome assembly projects. Hi-C exploits contact frequencies between pairs of loci to bridge and order contigs in draft genomes, resulting in chromosome-level assemblies. However, application of this approach is currently hampered by a lack of robust programs that are capable of effectively treating this type of data, particularly open source programs. We developed instaGRAAL, a complete overhaul of the GRAAL program, which has adapted the latter to allow efficient assembly of large genomes. Both GRAAL, and instaGRAAL use a Markov Chain Monte Carlo algorithm to perform Hi-C scaffolding, but instaGRAAL features a number of improvements including a modular polishing approach that optionally integrates independent data. To validate the program, we used it to generate chromosome-level assemblies for two brown algae, Desmarestia herbacea and the model Ectocarpus sp., and quantified improvements compared to the initial draft for the latter. Overall, instaGRAAL is a program able to generate, using default parameters with minimal human intervention, near-complete assemblies.

scaffolder GRAAL (Genome Re-Assembly Assessing Likelihood from 3D), a probabilistic tool that uses a Markov Chain Monte Carlo (MCMC) approach was able to generate the first chromosome-level assembly of an incomplete eukaryotic genome [12] by permuting DNA segments according to their contact frequencies until the most likely scaffold was reached (see also [13]). Since these proof of concept studies, the assemblies of many genomes of various sizes from eukaryotes [e.g. [14][15][16] and prokaryotes [17] have been significantly improved using scaffolding approaches exploiting Hi-C data.
Although GRAAL was effective on medium-sized or small (<100 Mb) eukaryotic genomes such as that of the fungus Trichoderma reesei [18], scalability limitations were encountered when tackling genomes whose complexity and size required significant computer calculation capacity. Furthermore, as was also observed with other Hi-C-based scaffolders, the raw output of GRAAL includes a number of caveats that need to be corrected manually to obtain a finished genome assembly. To tackle these limitations, we developed instaGRAAL, an enhanced, open-source program optimized to reduce the computational load of chromosome scaffolding and that includes polishing steps to automatically complete the assembly process. The polishing, which aims to minimise assembly errors, can exploit available genetic linkage data.
InstaGRAAL was applied to the 214 Mb and 500 Mb haploid genomes of the common brown algae Ectocarpus sp. and Desmarestia herbacea, respectively. The Ectocarpus genome is currently only published in draft form [19] and the Desmarestia genome was previously unpublished. Brown algae are a group of complex multicellular eukaryotes that have been evolving independently from animal and land plants for more than a billion years.
Ectocarpus sp. was the first species within the brown algal group to be sequenced, as a model organism to investigate multiple aspects of brown algal biology including the acquisition of multicellularity, sex determination, life cycle regulation and adaptation to the intertidal [20][21][22][23]. A range of genetic and genomic resources have been established for Ectocarpus sp. including a dense genetic map generated with 3,588 SNP markers [24], which was used to comprehensively validate the instaGRAAL assembly.

From GRAAL to instaGRAAL
The core principles of GRAAL and instaGRAAL are similar: both exploit a MCMC approach to perform a series of permutations (insertions, deletions, inversions, swapping, etc.) of genome fragments (referred to here as a 'bins', see Methods) based on an expected contact distribution [12]. The parameters (A, α and δ) that describe this contact distribution are first initialized using a model inspired by polymer physics [25]. This model describes the The parameters are then iteratively updated directly from the real scaffolds once their size increases sufficiently [12]. Each bin is tested in several positions relative to putative neighbouring fragments. The likelihood of each arrangement is assessed from the simulated or computed contact distribution, and the arrangement is either accepted or rejected [12].
This analysis is carried out in cycles, with a cycle being completed when all the bins of the genome have been processed in this way. Any number of cycles can be run iteratively and the process is usually continued until the genome structure ceases to evolve, as measured by the evolution of the parameters of the model. The core functions of the program use Python libraries as well as the CUDA programming language, and therefore necessitate a NVIDIA graphics card with at least 1 Gb of memory.
The technical limitations of GRAAL were i) high memory usage when handling Hi-C data for large genomes (i.e. over 100 Mb), 2) difficulties when installing the software, and 3) the need to adjust multiple ad hoc parameters to adapt to differences in genome size, read coverage, program that is lighter in resources, more flexible, and more user-friendly.
Other problems encountered with the original GRAAL program included 1) the presence of potential artefacts introduced by the permutation sampler, such as spurious permutations (e.g. local inversions) or incorrect junctions between bins; 2) difficulties with the correct integration of other types of data such as long reads and 3) difficulties with handling sequences that were either too short, highly repeated or with low coverage. We addressed these points by identifying and putting aside these problematic sequences during a filtering step. These sequences are subsequently reinserted into the final scaffolds, whenever possible (see Methods), with the help of linkage data when available. Overall, when compared to the raw GRAAL output, the resulting "polished" assemblies were significantly more complete and more faithful to the actual chromosome structure.

Scaffolding of the Ectocarpus sp. chromosomes with instaGRAAL
To test and validate instaGRAAL, we generated an improved assembly of the genome of the model brown alga Ectocarpus sp.. A reference genome consisting of 1,561 scaffolds generated from Sanger sequence data is available [20]. A Hi-C library was generated from a clonal culture of a haploid partheno-sporophyte carrying the male sex chromosome using a GC-neutral restriction enzyme (DpnII  (Fig 1a) (see Methods for details on the experimental and computational steps).
The resulting contact map in sparse matrix format was then used to initialize instaGRAAL along with the restriction fragments (RFs) of the reference genome (Fig 1a-b) (see Table S1 for an example of sparse file matrix). (a) (from left to right): i) The input data to be processed, paired-end reads to be mapped onto the Ectocarpus. sp. draft assembly; ii) Raw contact map before binning: each pixel is a contact count between two restriction fragments (RF); iii) Raw contact map after binning: each pixel is a contact between a determined numbers of RFs (see b). (b) Schematic description of one iteration of the binning process over 10 restriction fragments (arrows).
From left to right: i) initial contact map, each pixel is a contact count between two RFs; ii) filtering step: RFs either too short or presenting a read coverage below one standard deviation below the mean are discarded; iii) binning step (1 bin = 3RFs): adjacent RFs are pooled by three, with sum-pooling along all pixels in a 3x3 square.

Convergence of the Ectocarpus sp. assembly towards 27 major scaffolds
Given the probabilistic nature of the algorithm, we evaluated the program's consistency by running it three times with different resolutions. Briefly, we filtered out RFs that were shorter than 50 bp and/or whose coverage was one standard deviation below the mean coverage.
Then, we sum-pooled (or binned) the sparse matrix by groups (or bins) of three RFs five times, recursively (Fig 1a-b). Each recursive instance of the sum-pooling is subsequently referred to as a level of the contact map. A level determines the resolution at which permutations are being tested: the higher the level, the lower the resolution, the longer the sequences being permuted and, consequently, the faster the computation. The binning process is shown in Fig. 1b. Regarding Ectocarpus sp., we found that level four (bins of 81 RFs) was an acceptable balance between high resolution and fast computation on a desktop computer with a GeForce GTX TITAN Z graphics card. Moreover, whether instaGRAAL was run at level four, five or six (equivalent to bins of 81, 243 and 729 RFs respectively), all assemblies quickly (~6hrs) converged towards similar genome structures (Fig. 2a). was between four and ten times longer than the combined length of the remaining sequences ( Fig. 3). This strongly suggests that the 27 scaffolds correspond to chromosomes, a number consistent with karyotype analyses [26]. Taken together, these results indicate that instaGRAAL successfully assembled the Ectocarpus sp. genome into chromosome-level scaffolds. As the supplementary movie suggests, scaffold-level convergence is visible after only a few cycles, indicating that instaGRAAL is able to quickly determine the global genome structure most likely to fit the contact data. The remainder of the cycles are devoted to intra-chromosomal refinement. After filtering, and prior polishing, 27 main scaffolds (red bars) or putative chromosomes were obtained. The dotted green horizontal line represents the proportion of the filtered genome that was not integrated into the main 27 scaffolds and represent less than 0.6% of the initial assembly. Each scaffold presents, after normalization, a high quality Hi-C profile with features that are typical of eukaryotic genomes ( Figure S2).

Polishing the chromosome-level assembly of the Ectocarpus sp. genome
InstaGRAAL also includes a number of procedures that aim to correct modifications of the reference assembly contigs introduced during the Hi-C scaffolding (Fig. 4).
Polishing procedure (top to bottom): i) in silico restriction of the genome and binning, yielding a set of bins; ii) reordering of all bins into scaffolds without reference to their contig of origin; typically, groups of bins from the same contig naturally aggregate, but some bins get scattered to other scaffolds (e.g. bin 13, pink arrow), while others will be "flipped" with respect to the original assembly (e.g. bin 4, red arrows); iii) reconstruction of the original contigs by relocating scattered bins next to the biggest bin group; iv) bins in the original contigs are oriented according to their original consensus orientation.
These modifications principally involve discrete inversions or insertions of DNA segments (typically corresponding to single bins or RFs) (see also [12]). Such alterations are inherent to the statistical nature of instaGRAAL, which will occasionally improperly permute neighbouring bins because of the high density of contacts between them. These are part of a broader set of assembly errors (subsequently referred to as 'misassemblies') that we detected by aligning the reference contigs with the instaGRAAL scaffolds and analysing the mapping results using QUAST. We corrected these errors as follows: first, all bins processed by instaGRAAL that belonged to the same contig were constrained to their original orientation (Fig 4). If a contig was split across multiple scaffolds, the smaller parts of this contig were relocated to the largest one, respecting the original order and orientation of the bins. Then we reinserted whenever possible sequences that had been filtered out prior to instaGRAAL processing (e.g. contig extremities with poor read coverage; see Methods and Marie-Nelly et al., 2014 [12]) into the chromosome level scaffold at their original position in the contig of origin. A total of 3,832,980 bp were reinserted into the assembly this way. These simple steps alleviated artificial contig truncations observed with the original GRAAL program.
Some filtered bins had no reliable region to be associated with post scaffolding, because the entire initial contig they belonged to had been filtered. These sequences, which were left as-is and appended at the end of the genome, consisted of 543 scaffolds spanning 3,141,370 bp, i.e. < 2% of the total DNA. Together, these steps removed all the misassemblies detected by QUAST (Table 1). To further validate the assembly, we exploited genetic linkage data to search for potential translocations between scaffold extremities [24]. This analysis, now implemented as an option in instaGRAAL, detected such events in the unpolished version but none in the polished assembly. The polished instaGRAAL assembly is therefore fully consistent with the genetic recombination data, confirming the efficiency of the approach.

Comparisons with previous Ectocarpus sp. genome assemblies and validation of the instaGRAAL assembly
We compared the polished instaGRAAL assembly (subsequently referred to as the polished assembly) with three earlier assemblies of the Ectocarpus sp. genome (Table1 and Table   S2): 1) the reference assembly generated using Sanger sequencing data [20], mentioned above, which was assumed to be highly-accurate but highly fragmented (1,561 scaffolds); 2) an assembly generated by combining genetic recombination data and the Sanger assembly [19,24] (pseudochromosome assembly) and 3) an assembly generated by running the original GRAAL program on the reference genome (GRAAL assembly).
We aligned the instaGRAAL, pseudochromosomal and GRAAL assemblies to the reference assembly to detect misassemblies and determine whether the genome annotations (362,919 features) were conserved. We then validated each assembly using genetic linkage data (Methods). For each assembly, we computed the following metrics: the number of misassemblies, the fraction of conserved annotations, ortholog completeness and cumulative length/Nx distributions (Table 1). These assessments were carried out using BUSCO [27] for orthologue completeness ( Figure S1) and QUAST-LG's validation pipeline [28] for the other tests. QUAST-LG is an updated version of the traditional QUAST pipeline specifically designed for large genomes. We followed the terminology used by both programs, such as the BUSCO definition of ortholog and completeness, as well as QUAST's classification system of contig and scaffolds misassemblies.
The polished instaGRAAL assembly was of better quality than both the pseudochromosome and GRAAL assemblies (Table 1 and Figure S1). The polished assembly incorporated 795 of the reference genome scaffolds (96.8% of the sequence data) into the 27 chromosomes based on the high density genetic map [19], compared to 531 for the pseudochromosomal assembly (90.5% of the sequence data). Moreover, this assembly contained fewer misassemblies, retained more annotations and was more complete (both in terms of k-mers and BUSCO ortholog content). For some metrics the differences were marginal, but always in favour of the instaGRAAL assembly. BUSCO completeness was similar (76.2%, 76.9%, 77.6% for the GRAAL, pseudochromosomal and instaGRAAL assemblies, respectively, Figure S1), and an improvement over the 75.9% of the reference. These absolute numbers remain quite low, presumably because of the lack of a set of orthologs well adapted to brown algae.
All quantitative metrics, such as N50, L50 and cumulative length distribution, increased dramatically when compared with the reference genome ( contigs was integrated into the 27 largest scaffolds after instaGRAAL processing. Overall, the analysis indicated that many of the rearrangements found in the pseudochromosome assembly were potentially errors, and that both GRAAL and instaGRAAL were efficient at placing large regions where they belong in the genome, albeit less accurately for GRAAL and in the absence of polishing. These statistics underline the importance of the postscaffolding polishing steps, and the usefulness of a program that automates these steps.

Comparison between the Ectocarpus sp. instaGRAAL and pseudochromosomal assemblies
Compared to the pseudochromosomal assembly, the instaGRAAL assembly lost 23 scaffolds but gained 287 that the genetic map was unable to anchor to chromosomes (Table   S2). We observed few conflicts between the two assemblies. One major difference is that instaGRAAL was able to link the 4 th and 28 th pseudochromosomes that were considered to be separate by the genetic map [24] because of the limited number of recombination events observed. The fusion in the instaGRAAL assembly is consistent with the fact that the 28 th pseudochromosome is the smallest of the linkage groups, with only 54 markers over 41.8 cM and covering 3.8 Mb. The 28 th pseudochromosome has a very large gap which might reflect uncertainty in the ordering of the markers. Interestingly, this gap is located at one end of the linkage group, precisely where instaGRAAL now detects a fusion with the 4 th pseudochromosome. In addition, the fact that there is no mix between the 4 th and 28 th pseudochromosomes on the merged instaGRAAL chromosome but rather a simple concatenation, suggests that the genetic map was unsuccessful in joining those two linkage groups, but that instaGRAAL correctly assembled the two pseudochromosomes (see Table   S3 for correspondences between pseudochromosomes and instaGRAAL super scaffolds).
InstaGRAAL was also more efficient than the genetic map in orienting scaffolds (Table S2).
Among the scaffolds that were oriented in the pseudochromosomal assembly, about half of the 'plus' orientated were actually 'minus' and vice versa. The limited number of markers detected in the scaffolds anchored to the genetic map was likely the reason for this high level of incorrect orientations.

Scaffolding of the Desmarestia herbacea genome
To test and validate instaGRAAL on a second, larger genome, we generated an assembly of the haploid genome of D. herbacea, a brown alga that had not been sequenced before.
We set up the assembly pipeline and subsequent scaffolding from raw sequencing reads to assess the robustness of instaGRAAL with de novo, non-curated data. currently being sequenced will further refine this analysis, and improve the assembly.

Comparisons with existing methods
To date, only a limited number of Hi-C based scaffolding programs are publicly available, and, as far as we can tell, no detailed comparison has been performed between the existing programs to assess their respective qualities and drawbacks. In an attempt to benchmark instaGRAAL, we ran SALSA2 [31] and 3D-DNA on the same Ectocarpus sp. and Desmarestia herbacea reference genome and Hi-C reads. 3D-DNA is a scaffolder that was hallmarked with the assembly of Aedes aegypti and SALSA2 is a recent program with a promising approach that directly integrates Hi-C weights into the assembly graph. For Ectocarpus sp., SALSA2 ran for nine iterations and yielded 1,042 scaffolds, with an N50 of 6,552,506 (L50 = 11). Its BUSCO-completeness was 77.6%, a level identical to that obtained with instaGRAAL. Overall the metrics were satisfactory but SALSA2 was outperformed by instaGRAAL post-polishing. The contact map of the resulting SALSA2 assembly, displayed noticeably unfinished scaffolds ( Figure S4 and S5). This, coupled with a lower N50 value, suggests that instaGRAAL is more successful at merging scaffolds when appropriate.
We computed similar size and completeness statistics for the final instaGRAAL D. herbacea assembly and compared these to the values obtained with SALSA2 and 3D-DNA. We also mapped the Hi-C reads onto all three final assemblies in order to qualitatively assess the chromosome structure. The results are summarized in Table S5.  Figure S4). Notably, SALSA2 created a number of poorly-supported junctions to generate chromosomes, whereas 3D-DNA failed to converge towards any kind of structure. In contrast, although the instaGRAAL final assembly still contains contigs that are incorrectly positioned, a coherent structure corresponding to 40 scaffolds (including contaminants) emerged ( Figure S4). One possibility is that the de novo MaSuRCA assembly was lowquality, which could have resulted in alignment errors that disrupted the contact distribution and subsequent Hi-C scaffolding. Another possible explanation for these differences is that it remains difficult to dissect all the options and tuneable parameters of these scaffolders, and therefore that we did not find the optimal combination with respect to the D. herbacea draft assembly. Nevertheless, these results highlight the robustness of instaGRAAL which was able to successfully scaffold the D. herbacea genome using default parameters.

Discussion
InstaGRAAL is a Hi-C scaffolding program that can process large eukaryotic genomes.
Below we discuss the improvements made to the program, its remaining limitations and the steps that will be needed to tackle them.

Reference-based polishing
An important improvement of instaGRAAL compared to GRAAL relates to post-scaffolding polishing. Local misassemblies, e.g. local bin inversions or disruptive insertions of small scaffolds within larger ones, are an inevitable consequence of the algorithm's most erratic random walks. These small misassemblies are retained because flipping a bin doesn't markedly change the relative distance of an RFs relative to its neighbours, and because small scaffolds typically carry less signal and therefore exhibit a greater variance in terms of acceptable positions. Depending on the trust put in the initial set of contigs, one may be unwilling to tolerate these changes as well as "partial translocations", i.e. the splitting of an original contig into two scaffolds. The prevalence of such mistakes can be estimated by comparing the orientation of bins relative to their neighbours in the instaGRAAL assembly vs. the original assembly. Our assumption is that, if a single bin was flipped or split by instaGRAAL, this was likely a mistake that needed to be corrected by the polishing procedures. Consequently, we chose to remain faithful to the reference contigs, given that the initial Ectocarpus sp. reference was based on very accurate Sanger reads. Polishing therefore aims at reinstalling the initial contig structure and orientation while preserving to a maximum extent the overall instaGRAAL scaffold structure.
In addition, polishing reintegrates into the assembly the bins removed during the initial filtering process according to their position along the original assembly contigs. Most filtered bins corresponded to the extremities of the original contigs, because their size depended on the position of the restriction sites within the contig, or because they consisted of repeated sequences with little or no read coverage. The tail filtering polishing step inserts these bins back at the extremities of these contigs in the instaGRAAL assembly.
The combination of a probabilistic algorithm with a deterministic polishing step provides robustness to instaGRAAL. First, the MCMC step identifies, with few prior assumptions, a high-likelihood family of genome structures, almost always very close to the correct global scaffolding. The polishing step combines this result with prior assumptions made about the initial contig structures generated through robust, established assembly programs, refining the genomic structure within each scaffold. To give the user a fine-grained degree of control over polishing, the implementation itself into instaGRAAL is split into independent modules that each make an assumption about the initial contig structure necessary to perform the correction: the 'reorient' module assumes that the initial contigs do not display inversions, the 'rearrange' module assumes that there are no relocations within contigs.
We would like to underline that, despite the improvements brought about by these new procedures, instaGRAAL assemblies remain perfectible, notably because of the reliance on the quality of the original contig assembly. However, chromosome-level assemblies consistent with known genetic maps and other datasets were reached every time the program was applied on a new species.

Sparse data handling
The implementation of a sparse data storage method in instaGRAAL allows much more intense computation than with GRAAL. Because the majority of map regions are devoid of contacts, instaGRAAL essentially halves the order of magnitude of both algorithm complexity and memory load, i.e. they increase linearly with the size of the genome instead of geometrically. This improvement potentially allows the assembly of Gb-sized genomes in five to six days using a desktop computer (and faster with a larger computational resource).

Filtering
Variations in GC% along the genome, and/or other genomic features, can lead to variation in Hi-C read coverage and impair interpretation of the Hi-C data. Correction and attenuation procedures that alleviate these biases are therefore commonly used in Hi-C studies [32][33][34]. However, these procedures are not compatible with instaGRAAL's estimation of the contact distribution (for more see [35]). A subset of bins will therefore diverge strongly from the others, displaying little if no coverage. A filtering step is needed to remove these bins as they would otherwise impact the contact distribution and the model parameter estimation.
These disruptive bins represent a negligible fraction of the total genome (< 3% of the total genome size of Ectocarpus sp., for instance) and are reincorporated into the assembly during polishing. On the other hand a subset of bins representing small, individual scaffolds are not reinserted during polishing, and are added to the final assembly as extra-scaffolds (as in all sequencing projects). Additional analyses and new techniques such as long or linked reads are needed to improve the integration of these scaffolds into the genome.

Resolution
The binning procedure will influence the structure of the final assembly as well as its quality.
For example, low level binning (e.g. one bin = three RFs) will lead to an increased number of bins and a large, sparse contact map with a low signal-to-noise ratio, where many of the bins display poor read coverage as on average they will have fewer contacts with their immediate neighbours. Because of the resulting low signal-to-noise ratio, an invalid prior model will be generated and, when referring to this model, the algorithm will fail to scaffold the bins properly, if at all. Moreover, due to its probabilistic nature, the algorithm will generate a number of false positive structural modifications such as erroneous local inversions or permutations of bins. The numerous bins will create more genome structures to explore to handle all the potential combinations, and exploring this space until convergence will take longer and be computationally demanding.
An optimal resolution is therefore a compromise between the bin size, the coverage, and the quality of the original contig assembly. Although a machine powerful enough operating on an extremely contact-rich matrix would be successful at any level, it is unclear whether such resources are necessary. Our present assemblies (e.g. 1 bin = 81 RFs for both; Material and Methods) had good quality metrics after a day's worth of calculation on a standard desktop computer. Moreover convergence was qualitatively obvious after a few cycles. This suggests that more computational power yields diminishing returns, and therefore that appropriate polishing is a more efficient approach to correct any remaining misassemblies.

Binning
The fragmentation of the original assembly used to generate the initial contact map has obviously a substantial effect on the quality of the final scaffolding. Because binning cannot be performed beyond the resolution of individual contigs, however small they may be, there is a fixed upper limit to the a scale at which a given matrix can be binned. A highly fragmented genome with many small contigs will necessarily generate a high-noise, highresolution matrix. Attempts to reassemble a genome based on such a matrix will run into the problems discussed above (resolution). This limitation can be alleviated, to some extent, by discarding the smallest contigs, with the hope that the remaining contigs will cover enough of the genome. The contigs that are removed can be reintegrated into the assembly during the polishing steps. This ensures an improved Nx metric while retaining genome completeness. It should be noted, however, that the size of the contigs is only important insofar as they need to contain sufficient restriction sites, and each of the restriction fragments must have sufficient coverage. The choice of enzyme and the frequency of its corresponding site is thus crucial. For instance, with an average of one restriction site every 600 to 1,000 bp for DpnII, contigs as short as 10 kb may contain enough information to be correctly reassembled. The restriction map therefore strongly influences both the minimum limit on N50 and genome fragmentation.

Benchmarking
In order to test our tool against existing programs, we ran two scaffolders available online (SALSA2 and 3D-DNA) on our two genomic datasets. In all instances, instaGRAAL proved more successful at scaffolding both genomes. However, we have not extensively tested all the combinations of parameters of both programs, and acknowledge the difficulty in designing and implementing Hi-C scaffolding pipelines with extensive dependencies that compound the initial complexity of the task and add yet more configurable options to know in advance. Finding the correct combination of CUDA and Python dependencies to install instaGRAAL on a given machine can be challenging as well. Therefore, our benchmarking attempt should be rather seen as a way to stress the importance of implementing sensible default parameters that readily cover as many use cases as possible for the end user.
There is almost no doubt that both 3D-DNA and SALSA2, with the appropriate parameters and polishing steps, would produce satisfying scaffolding; on the other hand, knowing which input parameters has to be specified in advance is a nontrivial task, especially given the computational resources needed for a single scaffolding run. With instaGRAAL, we wish to combine the simplicity of a default configuration that works in most instances, with the flexibility offered by the power of MCMC methods.

Integrating information from the Hi-C analysis with other types of data
Aggregating data from multiple sources to construct a high-quality assembly remains a challenging problem with no systematic solution. As long read technologies become more affordable, there is an increasing demand to reconcile the scaffolding capabilities of Hi-C based methods with the ability of long reads to span regions that are difficult to assemble, such as repeated sequences. The most intuitive approach would be to perform Hi-C scaffolding on an assembly derived from high-coverage and corrected long reads, as was done for several previous assembly projects [14,36]. Alternative approaches also exist, such as generating Hi-C and long-read-based assemblies separately and merging them InstaGRAAL shows that high quality metrics can still be attained without the help of long reads, but it can nevertheless integrate them when necessary or available. Long reads are not the only type of data that may be used to improve assemblies. Linkage maps, RNAseq, optical mapping and 10X technology all provide independent data sources that can help improve genome structure and polish specific regions. The success of future assembly projects will hinge on the ability to process these various types of data in a seamless and efficient manner.

Preparation of the Hi-C libraries
The Hi-C library construction protocol was adapted from [7,37].

Contact map generation
Contact maps were generated from reads using the hicstuff pipeline for processing generic 3C data, available at https://github.com/koszullab/hicstuff. The backend uses the bowtie2 (version 2.2.5) aligner run in paired-end mode (with the following options: --maxins 5 -verysensitive-local). Alignments with mapping quality lower than 30 were discarded. The output was in the form of a sparse matrix where each fragment of every chromosome was given an unique identifier and every pair of fragments was given a contact count if it was nonzero.
Fragments were then filtered based on their size and total coverage. First, fragments shorter than fifty base pairs were discarded. Then, fragments whose coverage was less than one standard deviation below the mean of the global coverage distribution were removed from the initial contact map. A total of 6,974,350 bp of sequence was removed this way. An initial contact distribution based on a simplified a polymer model [25] with three parameters was first computed for this matrix. Finally, the instaGRAAL algorithm was run using the resulting matrix and distribution.
For the Ectocarpus sp. genome, instaGRAAL was run at level 4 (n = 81 RFs), 5 (n = 243 RFs) and 6 (n = 729 RFs). Levels 5 and 6 were only used to check for genome stability and consistency in the final chromosome count. Level 4 was used for all subsequent analyses.
All runs were performed for 250 cycles. The starting fragments for the analysis were the reference genome entirely fragmented into restriction fragments. The MCMC was run with 3 burn-in cycles. The same parameters were used for the Desmarestia herbacea genome.

Polishing of genome assemblies
The assembled genome generated by instaGRAAL was polished to remove misassemblies using a number of simple procedures that aimed to reinstate the local structure of the initial contigs where possible. Briefly, bins belonging to the same initial contig were juxtaposed in the same relative positions as in the starting assembly contig. Small groups of bins were preferentially moved to the location of larger groups when several such groups were present in the assembly. The orientations of sets of bins that had been regrouped in this manner were modified so that orientation was consistent and matched that of the majority of the group, re-orientating minority bins when necessary. Both steps are illustrated in Fig. 4.
Finally, fragments that had been removed during the filtering steps were reincorporated if they had been adjacent to an already integrated bin in the initial assembly. The remaining sequences that could not be reintegrated this way were appended as non-integrated scaffolds.

Validation metrics
Initial and final assembly metrics (Nx, GC distribution) were obtained using QUAST-LG [28].
Misassemblies were quantified using QUAST-LG with the minimap2 aligner in the back-end.
Ortholog completeness was computed with BUSCO (v3) [27]. Assembly completeness was also assessed with BUSCO. The evolution of genome metrics between cycles was obtained using instaGRAAL's own implementation.

Validation with the genetic map
The validation procedure with respect to linkage data was implemented as part of instaGRAAL. Briefly, the script considers a set of pseudochromosomes where regions are separated by SNP markers, and a set of Hi-C scaffolds where regions are bins separated by restriction sites. It then finds best-matching pairs of pseudochromosomes/scaffolds by counting how many of these regions overlap from one set to the other. Then, for each pair, the bins in the Hi-C scaffold are rearranged so that their order is consistent with that of the corresponding pseudochromosome. Such rearrangements are parsimonious and try to alter as little as possible. Since there isn't a one-to-one mapping from restriction sites to SNP markers, some regions in the Hi-C scaffolds are not present in the pseudochromosomes, in which case they are left unchanged. When the Hi-C scaffolds are altered this way, as was found in the case of the raw GRAAL assembly, the script acts as a correction. When the scaffolds are unchanged, as was the case with the instaGRAAL assembly, the script acts as a validation.

Benchmarking
For each genome, the 3D-DNA program was run using the run-assembly-pipeline.sh entrypoint script with the following options: -i 1000 --polisher-input-size 10000 --splitterinput-size 10000. The Hi-C data was prepared with the Juicer pipeline as recommended by 3D-DNA's documentation. The SALSA2 program was run with the -cutoff=0 option, and misassembly correction with the -clean=yes option. No expected genome size was provided. The program halted after 9 iterations for Ectocarpus sp., and 18 iterations for Desmarestia herbacea. Hi-C data was prepared with the Arima pipeline as recommended by SALSA2's documentation.

Software tool requirements
The instaGRAAL software is written in

Availability of data and materials
The raw, unfiltered datasets generated and analysed during the current study are available in the SRA repository, SRR8550777.

Competing Interests
InstaGRAAL is owned by the Institut Pasteur. The entire program and its source code are freely available under a free software license.

Supplemental Information
Chromosome-level quality scaffolding of brown algal genomes using InstaGRAAL, a proximity ligation-based scaffolder        Movie S1: the iterative scaffolding process can be visualized on a movie accessible through the following link. https://github.com/koszullab/ectocarpus_scripts/blob/master/images/matrix_evolution.gif Each frame corresponds to a cycle during which each fragment has been processed once.


extfieldsofstudy: Chemistry Medicine

s2fieldsofstudy: Biology Engineering

text:
Long-term exposure of activated sludge in chemostats leads to changes in microbial communities composition and enhanced biodegradation of 4-chloroaniline and N-methylpiperazine

Exposure history and adaptation of the inoculum to chemicals have been shown to influence the outcome of ready biodegradability tests. However, there is a lack of information about the mechanisms involved in microbial adaptation and the implication thereof for the tests. In the present study, we investigated the impact of a long-term exposure to N-methylpiperazine (NMP) and 4-chloroaniline (4CA) of an activated sludge microbial community using chemostat systems. The objective was to characterize the influence of adaptation to the chemicals on an enhanced biodegradation testing, following the OECD 310 guideline. Cultures were used to inoculate the enhanced biodegradability tests, in batch, before and after exposure to each chemical independently in chemostat culture. Composition and diversity of the microbial communities were characterised by 16s rRNA gene amplicon sequencing. Using freshly sampled activated sludge, NMP was not degraded within the 28 d frame of the test while 4CA was completely eliminated. However, after one month of exposure, the community exposed to NMP was adapted and could completely degrade it. This result was in complete contrast with that from the culture exposed for 3 months to 4CA. Long term incubation in the chemostat system led to a progressive loss of the initial biodegradation capacity of the community, as a consequence of the loss of key degrading microorganisms. This study highlights the potential of chemostat systems to induce adaptation to a specific chemical, ultimately resulting in its biodegradation. At the same time, one should be critical of these observations as the dynamics of a microbial community are difficult to maintain in chemostat, as the loss of 4CA biodegradation capacity demonstrates. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). msterdam, the Netherlands.

Introduction
The biodegradability of organic compounds is one of the most important properties that determine the environmental risk assessment of chemicals. Biodegradability can be determined using different tests and guidelines, such as the guidelines from the Organisation for Economic Cooperation and Development (OECD) (OECD, 2014). The results of these tests allow a fast categorization of the tested compounds from readily biodegradable to persistent, or recalcitrant, in the environment (ECHA, 2017). Ready biodegradability is an arbitrary classification of chemicals which have passed a specified type of screening tests for ultimate biodegradability (mineralization). Compounds identified as readily biodegradable are considered to be rapidly and completely biodegraded in aquatic environments under aerobic conditions (Comber and Holt, 2010;ECETOC, 2013;OECD, 2014). However, the commonly used ready biodegradability tests suffer from several well-known problems that can affect their outcome (Dick et al., 2016;Federle et al., 1997). One of these is the variability of the lag phase period, which is the time required for a community to start degrading a compound (Tor€ ang and Nyholm, 2005)). In some cases, specific factors such as the temperature or the occurrence of aerobic or anaerobic conditions could explain these variabilities. In many cases, however, the causes for the differences in biodegradation rate or lag-phase period are unknown (Greskowiak et al., 2017).
Adaptation of a microbial community that is used as inoculum to degrade the tested compound is one of the factors that potentially influences the lag-phase period (Itrich et al., 2015;Kowalczyk et al., 2015). A pre-adaptation step prior to a biodegradation test has been shown to decrease the number of false-negative results of the test (Ingerslev and Nyholm, 2000) and to reduce the lag-phase period (Tor€ ang and Nyholm, 2005). It has been proposed to include adaptation to chemicals in enhanced biodegradability tests (ECETOC, 2007), however, pre-exposure prior any standard RBT is not allowed under the current ECHA guidelines. Enhanced biodegradability tests use the same methodology than the RBT, however, with variation in their duration and inoculum density (ECETOC, 2007;Ott et al., 2019). The objective of enhanced biodegradability tests is to improve to the reliability of tests and to decrease the number of false negative and positive results. In general, however, there is a lack of information about the microbial community dynamics and evolution during adaptation periods, both in laboratory systems or in the environment (Birch et al., 2017). Chemostat systems may be used to mimic the natural adaptation process in the laboratory under defined and controlled conditions. In the present study, we investigated the influence of long-term exposure of bacterial communities to chemicals in a chemostat on the biodegradation testing outcome along with the adaptive responses of the bacterial community composition. For this purpose, we exposed an activated sludge community to two different chemicals that show inconsistent behaviour in biodegradability testing.
The first one of these chemicals, 4-chloroaniline (4CA), is an intermediate product formed during the production of various organic chemicals and polymers, such as pharmaceuticals, pesticides, dyes and others. Chloroanilines are known to be relatively toxic for humans and wildlife (Argese et al., 2001), but also for microorganisms (Gosetti et al., 2010). It has been shown that the toxicity of chloroanilines for activated sludge microbial communities can affect and inhibit the degradation of other organic chemicals (Zhu et al., 2011). 4CA was reported to be present in the low mg/l range in Dutch surface waters (Heugens and Verbruggen, 2009;Wickham, 2016), however, due to technical limitations, the concentration in activated sludge was not measured in this project. 4-chloroaniline shows erratic behaviour in ready biodegradability testing (Comber and Holt, 2010) as it may either pass or completely fail the test (Ingerslev and Nyholm, 2000). Hence, 4CA has been proposed as a reference chemical that would normally pass an enhanced ready biodegradability test but would fail any other standard tests (Comber and Holt, 2010). N-Methylpiperazine (NMP), the second test compound in this study, is member of the family of piperazines, which are widely used industrial building blocks for pharmaceuticals, plastics and other products (Cai et al., 2013). As the extent of use and occurrence of NMP is unknown, many uncertainties exist about its biodegradability, though unpublished results, from the ECHA database, indicate that NMP is recalcitrant to biodegradation (ECHA, 2013). Significant degradation was, however, observed after adaptation of the community to this compound in a semi-continuous activated sludge system (ECHA, 2013).
In the present study, the biodegradation capacities for 4CA and NMP of unexposed activated sludge communities and those exposed to the chemicals in chemostat cultures were compared using the OECD 310 ready biodegradability testing protocol (OECD, 2014) and an enhanced version of the test. Differences between the tests output were evaluated in the light of the changes in communities composition in time followed by 16s rRNA gene Illumina sequencing.

Sampling sites
The activated sludge used in this project was sampled in June 2016 from the aeration tank of the Amsterdam West wastewater treatment plant (WWTP) of the city of Amsterdam (The Netherlands). This facility operates since 2005 and has a hydraulic capacity of 8400 m 3 /hours in dry weather flow and 30,000 m 3 / hours in wet weather flow. Activated sludge was maintained under constant shaking (150 rpm) at 22 C (±1 C) until use and for a maximum of 1 week, in order to allow the inoculum to acclimate to the laboratory conditions, according to the OECD 310 guideline (OECD, 2014). After one week of incubation under constant shaking, the dry matter content of the sludge sample was determined by filtration of a known volume through a 1.2 mm glass microfibre filter GF/C (Whatman), after which the filters were dried overnight in an oven at 70 C and weighed.

Chemicals
N-methylpiperazine (CAS Number 109-01-3) and 4chloroaniline (CAS Number 106-47-8) were purchased from Sigma-Aldrich (Zwijndrecht, the Netherlands) with a >99% purity grade. Working and standard solutions were prepared in ultra-pure water (ELGA MK2-Analytic, Ede, the Netherlands) and stored at 4-8 C. Methanol and acetonitrile used for sample pre-treatment and chromatography had an analytical or LC-MS grade and were purchased either from Biosolve, Valkenswaard, the Netherlands.

Chemostat reactor
To cultivate and expose the activated sludge community, we used three chemostat systems of 0.5 L volume each. Chemostats contained a similar sterile growing medium and one additional test compound. Both systems exposed to either 4CA or NMP were fed with medium containing a nominal concentration of 1.5 mg/l of the test compound, while the third one was used as a control and was only fed with growth medium without added substance. A chemical concentration of 1.5 mg/l was chosen in order to represent a stable food source without competing with the other carbon and energy sources.
Fresh medium (10 L) was prepared every 10 d. Stock solutions of vitamins, 4CA or NMP were prepared using ultra-pure water and added under sterile condition to the autoclaved medium via a 0.2 mm silicone membrane filter (Whatman plc). System vessels for media, cultures and effluents were made of glass, tubes were made of stainless steel, air tubes of silicones, and pump tubes of silicon. The system was autoclaved (121 C, 211 kPa, 21 min) before incubation. During the experiment, the pH (7.4) of the effluent was regularly measured with a pH probe. A constant bubbling system using compressed air (filtered through a 0.2 mm PTFE membrane) provided oxygen to the culture. All three chemostat systems were inoculated at the same time with the same activated sludge volume (5 ml) and were maintained in the dark at 22 C (±1 C) under constant stirring (300 rpm). Pumps were immediately activated after inoculation. Flow of the medium influent (20 ml/h) was controlled by a Gilson Minipuls-3 peristatic pump to maintain a dilution rate of 0.04 h À1 . Samples of the culture were taken from the sampling port under sterile conditions.

Biodegradability testing
The biodegradation tests were performed following the OECD 310 ready biodegradability protocol (OECD, 2014) for the unexposed community and an enhanced version of the test for the exposed communities. During these tests, primary biodegradation was followed by the direct quantification of the compounds by LC-MS/MS analysis. Tested substances were dissolved at a nominal concentration of 20 mg/l in the mineral medium buffered at pH 7.4. The mineral medium was a mix of diverse salts and ammonium chloride as nitrogen source: 85 mg/l KH 2 PO 4 ; 217.5 mg/l K 2 HPO 4 ; 334 mg/l Na 2 HPO 4 $2H 2 O; 5 mg/l NH 4 Cl; 36.4 mg/l CaCl 2 $2H 2 O; 22.5 mg/l MgSO 4 $7H 2 O; 0.25 mg/l FeCl 3 $6H 2 O.
Bottles were directly inoculated by adding 10 4 e10 5 CFU/mL (colony forming units) of activated sludge or chemostat samples. CFU numbers were determined on nutrient agar by serial dilution, after 48 h of incubation at room temperature. In parallel, abiotic degradation was assessed by inoculating a test with an autoclaved inoculum, diluted to the same level as the other inocula. Bottles were sealed with rubber seals and aluminium caps and headspace was replaced by synthetic air, composed of 80% N 2 and 20% O 2 . Incubation was performed in the dark at 22 C (±1 C) and under constant shaking (150 rpm).
Data analyses of the biodegradation data was performed using Excel and R studio. When possible, models were used to determine the time necessary to reach 10%, 50% and 90% of degradation. The DRC R package (Ritz et al., 2015, p.) was used to fit a threeparameters log-logistic models to the data and to determine the different degradation times.

LC-MS/MS
To quantify 4CA and NMP in liquid samples, two LC-MS/MS methods were established using a Shimadzu LC20 instrument (Shimadzu, 's-Hertogenbosch, the Netherlands), equipped with a SIL-20ACXR auto sampler and connected in series with a Sciex 4000 QTRAP MS/MS system (AB SCIEX, MA, USA). Separation was performed on a reversed phase C18 column (1.6 mm, 50 Â 2.0 mm, Shim-Pack, Shimadzu) with a security guard ultra-cartridge as guard column for 4CA and a HILIC column (1.7 mm, 50 Â 2.1 mm, amide BEH, Waters, Etten Leur, the Netherlands) for NMP. For both methods, the column oven was set at 35 C and the injection volume was 20 ml, Two pumps (model LC-20ADXR) were used to create a gradient of two solvents. The mobile phase of the reversed phase method, used to quantify 4CA, consisted of (A) ultra-pure water þ 0.1% acetic acid, and (B) methanol. The flow rate of the mobile phase was of 0.3 ml/min. After injection, the B concentration was increased from 10% to 100% within 2.5 min, then B held at 100% for 2 min before returning to the initial concentration of 10% for 5 min before injection of the next sample. For NMP, the flow rate was 0.4 ml/min. The mobile phase consisted of (A) ultrapure water and (B) acetonitrile, both buffered with 5 mM of ammonium formate and 0.075% (v/v) of formic acid. Samples were ionized using electrospray ionization (ESI, Turbospray) operating in positive mode. Multiple-reaction monitoring mode (MRM) was applied with the first transition (TR1) used for quantification and a second transition (TR2) used for qualitative confirmation. All of the MRM parameters can be found in Table 1. The following parameters were used for the mass spectrometer: Collision gas: 6 l/h, curtain gas: 10 l/h, ion-spray voltage: 4000 C, Temperature: 500 C, Ion source gas 1: 40 l/h and ion source gas 2: 50 l/h. Before injection, the samples were diluted in the relevant solvents (e.g. 4CA: water; NMP: acetonitrile) and filtered using a 0.2 mm polypropylene filter (Filter-BIO, China). Calibration standards of both compounds were prepared and injected in order to quantify the substances in the samples. These standards were prepared under the same conditions as the samples. Data acquisition and analysis were performed with AB SCIEX Analyst software (Ver. 1.5.1).

Physiological profiling of bacterial communities
Physiological profiling of the microbial communities was done with Biolog EcoPlates™, which are 96 wells plates that contain 31 different carbon and free energy sources in triplicate. Briefly, microbial communities derived from chemostat or activated sludge were diluted to 10 À3 CFU/ml. Ecoplates were directly inoculated with 150 ml of diluted sample per well, and were incubated at 22 C in the dark for one week. Absorbance at 590 nm was measured with a BMG labtech SPECTROstat Nano absorbance microplate reader (Isogen life Science, The Netherlands). The average well-colour development (AWCD) was calculated as described by (Garland and Mills, 1991) and following this formula (AWCD ¼ P A 590i /31) in which A 590 is the well colour development minus that of the control. The AWCD was calculated for each plate and for each of the six substrate groups: carbohydrates, carboxylic, amines, amino acids phenolic compounds and polymers. The Shannon-Weaver index (H), was calculated as following: H' ¼ P -pi (lnpi) in which pi is the normalized value of well development (Gryta et al., 2014). Normalization of the absorbance data were performed to reduce the variability due to differences in inocula densities. Measured absorbance of each well, at each time point, were divided by the initial absorbance of the corresponding carbon source (Hackett and Griffiths, 1997;Muñiz et al., 2014). All statistical analyses were performed with R (v3.5.1) (R Core Team, 2016).

DNA extraction and PCR amplification
DNA samples were extracted and purified using the PowerSoil kit (QIAGEN Benelux B.V.), following the supplier's protocol. Samples were stored at À80 C until further analysis. DNA concentrations were determined using the Qubit dsDNA HS Assay kit (Thermo Fischer Scientific, Waltham, MA USA) and the Qubit fluorometer (Thermo Fischer Scientific, Waltham, MA USA). Used primers (8 forward primers and 16 reverse primers) (Appendix A) target the V3eV4 region of bacterial 16S rRNA genes and were yield at~550 bp fragment during PCR amplification. PCR products were purified using an Agencourt AMPure XP magnetic bead (Beckman Coulter Nederland B.V.). Agarose gel electrophoresis was used to control the purification step. Samples were diluted to 1 ng/ml prior to pooling and subsequent analysis by the Illumina MiSeq platform.
Alpha (Shannon H, inverse Simpson, Observed) and beta diversity analyses (Unifrac and Bray-Curtis distance) were carried out using the Vegan package (Oksanen et al., 2017) of R (R Core Team, 2016). The Shannon diversity index (H) (Surovtseva et al., 1985) is a measure of diversity combining both abundance and evenness. Further, data were normalized using the variance-stabilizing transformations, using the R packages Phyloseq and DeSeq2 (ver. 1.22.1) (Love et al., 2014) and following examples from (McMurdie and Holmes, 2014). Dissimilarities of community composition were assessed by non-metric multidimensional scaling (NMDS) using the Bray-Curtis dissimilarity on the normalized data. Principal component analysis (PCA) and Unifrac (Lozupone and Knight, 2005) was also used on normalized data to compare the communities by incorporating phylogenetic information. Permutational multivariate analysis of variance (PERMANOVA) were performed using the Adonis function from the R package Vegan (Oksanen, J., et al., 2017). Finally, differential abundance analyses were performed with Phyloseq and DeSeq2. Differential abundance is defined as the difference in mean relative abundances between groups or conditions (McMurdie and Holmes, 2014). Only groups showing a significant difference (p adj < 0.01) with the original inoculum, in each community, were displayed.

Biodegradation in OECD 310 tests
The biodegradation profile in time of NMP, measured by LC-MS/ MS according to the OECD 310 biodegradability guidelines, is shown in Fig. 1A&B. The colours of the lines indicate the exposure time of the inoculum to NMP in the chemostat. The compound was not degraded after 28 d of incubation in the abiotic control, indicating that abiotic degradation did not take place. Test batches incubated with the original activated sludge were able to degrade approximately 10% of the original NMP content after 28 d, indicating that NMP can be considered as recalcitrant to biodegradation in this test. Moreover, the activated sludge community cultivated in chemostat culture without being exposed to NMP showed no sign of significant biodegradation activity, not even after 9 months (Fig. 1B). In contrast, after 1 month of exposure in the chemostat, the microbial community was able to completely degrade the original concentration of NMP within 10 d. We observed the same result after three and nine months of exposure. The inoculum exposed for nine months took only 6.63 ± 1.81 d to reach 90% degradation of NMP (Appendix B). It should be noted that the density of the inoculum, measured in CFU/ml, was similar in every test.
The biodegradation of 4CA in OECD 310 biodegradability tests, using inocula from activated sludge communities exposed to 4CA in chemostat culture, is shown in Fig. 1C and the results using unexposed inocula in Fig. 1D. We can observe that the results from the abiotic control do not correspond to the three-parameter log-logistic model and that the degradation pattern fluctuated over time with high standard deviations. However, tests with activated sludge communities exposed for one month as inocula reached 90% of degradation within 28 d. In contrast, communities exposed for three and nine months could only degrade 45% of the 4CA added before the end of the test, indicating loss of the capacity to degrade 4CA.

Community dynamics
The Illumina MiSeq platform produced 3,488,627 sequences in total and an average of 61,203 reads per samples, that belong to 17,650 operational taxonomic units (OTUs). OTUs with a number of reads lower than 3 were discarded, in order to filter single and artefact OTUs, reducing the total number to 5666. The abundance and taxonomic data were used to produce the different alpha and beta diversity indices display in Fig. 2.
Alpha diversity indices are used to follow the community evolution in each chemostat over time. The observed richness represents the number of unique species observed in our samples, while the Shannon diversity index (H) combines both richness and the evenness of the communities. The Simpson diversity index also includes both richness and evenness of the community. However, we used the inverse Simpson index since this diversity index is not influenced by the sampling effort in contrast to the regular Simpson. The chemostats were inoculated at the same time using the same sample of activated sludge, which is confirmed by similar extents of alpha diversity at the start of culturing. We then observe that their diversity quickly diverges in time, and even fluctuates. Only the chemostats exposed to 4CA and NMP seem to stabilize at the end of the experiment. However, ANOVA and Kruskal Wallis tests, show that these indices do not differ significantly over time.
Based on the Bray Curtis dissimilarity (Fig. 2B), it appears that the cultivated communities are different from each other. The microbial community structure of each chemostat at the first day of the experiment, just after inoculation, showed considerable similarity, as can be observed from the cluster in the higher right-hand corner of the NMDS plot in Fig. 2B. Subsequently, the communities quickly started to diverge. This result is confirmed by the output from the PERMANOVA Adonis test (Appendix C). The three chemostat communities do not share the same centroid and the permutation test for homogeneity of multivariate indicates that the Adonis results are not due to group dispersions. We can therefore conclude that the communities are indeed different. Furthermore, based on the PCA in Fig. 2C, we observe differences over time in each chemostat. As expected, cultivating a microbial community in a chemostat with artificial media changes its composition. These results suggest that the exposure to 4CA or NMP, as well as the cultivating conditions, are influencing the community's evolution.
Some members of these genera were found in our samples. Acinetobacter and Pseudomonas together represented approximately 30e40% of the total community at the start (Fig. 3), while after 1 month of exposure in chemostat 4CA, species from the genus Delftia became more abundant until they represented approximately 40% of the community. Furthermore, these observations are supported by the data in Fig. 4, which displays every OTU in chemostat 4CA, whose abundance significantly (p adj < 0.01) differs from the original inoculum. An OTU with a negative log2-FoldChange (Fig. 4A) indicates that the abundance significantly decreased (Fig. 4B), while a positive one indicates an increase (Fig. 4C). Several OTUs from the genera Acinetobacter and Pseudomonas are amongst the numerous OTUs whose abundance decreased. While one OTU from the genus Delftia increased over time, its contribution to the biodegradation of 4CA, in our experiment, is unknown.
Our data do not enable the identification of NMP-degrading microorganisms. However, the results shown in Figs. 2 and 5 show a shift in the community composition over time. The NMP community is different to the other chemostats and the initial activated sludge used. Fig. 5 shows that some OTUs from the genera Variovorax, Terrimonas, Delftia, Mycobacterium, an uncultured bacterium and an OTU from the order Burkholderiales increased in abundance after one months of exposure. Among these genera, Delftia and Terrimonas increased significantly in abundance in the unexposed chemostat as well (Appendix D). Fig. 5B represents the OTUs that were outcompeted from the community exposed to NMP.

Community physiological profiling
Biolog Ecoplates can provide extensive information regarding the general use of carbon sources by the cultivated communities. We did not perform any kinetic analyses on our data as this can be influenced by the biomass, which was not controlled in this experiment. However, to compare the catabolic diversity among communities, the data were normalized, as described in the materials and methods section, in order to reduce the variability due to differences in inocula densities (Hackett and Griffiths, 1997;Muñiz et al., 2014). The Shannon Weaver index (H) was used to calculate the physiological diversity of each community.
Analysis of the variance (ANOVA), performed on the calculated Shannon-Weaver indexes (Table 2), indicates that the chemostat communities significantly differed from the original activated sludge after 1 month of cultivation (p < 0.001), with the exception  Relative abundance (%) of the known 4CA-degrading genera in chemostats over time. Relative abundance calculated based on the whole community abundance. Left: community exposed to 4CA; Right: community that was not exposed to 4CA. of chemostat blank, that significantly differed after 3 months. However, no significant differences could be measured between the chemostats, with the exception of the chemostat exposed to NMP after 9 months of exposure and the chemostat exposed to 4CA after 1 and 3 months.
Cluster analysis of the AWCD for each substrate group shows different patterns that may be shared between chemostats (Fig. 6). For example, after 3 months of exposure the NMP and 4CA chemostats are statistically not different from each other, based on ANOVA and the post hoc Tukey's test. Overall, carbohydrates were the less well degraded group and activated sludge was the worst performing community, showing the least AWCD development. It seems that the presence or absence of NMP and 4CA did not affect the community's physiological profile in a distinctive way. Exposed and unexposed communities could be significantly different at one time point, but this may no longer be the case at the next time point.

Discussion
Adaptation of microbial communities to chemicals being persistent in laboratory systems is a known and documented process (Itrich et al., 2015;Poursat et al., 2019a). Despite chemostats being suitable to study the adaptive evolution of microorganisms  (Gresham and Hong, 2014), very few studies used them as model system to induce adaptation. Previous investigation showed that communities exposed in chemostat could adapt and influence the results of biodegradability tests (Thouand and Block, 1993). Adapted inocula enhanced the reproducibility of biodegradability tests, decreased the lag time and increased the biodegradation rate of the test chemicals (Thouand and Block, 1993). Based on our results, we can also confirm that exposure of microbial communities in chemostat systems modify their biodegradation abilities. Moreover, it should be noted that the community cultivated in the unexposed chemostat, and only fed with acetate, yeast extract and micronutrients, was unable to degrade NMP and completely lost its ability to degrade 4CA. Loss or changes of biodegradation potential by cultivated activated sludge community, in absence of the test chemical, have been observed in the MITI tests (OECD 301C and 302 C) (V azquez-Rodríguez et al., 2011) and in chemostat systems Block, 1993, Poursat et al., 2019b). Long term culturing in the absence of the test chemical, without sludge recycling, will lead to changes in the community structure, which might result in a loss of function. Moreover, community diversity in culture system is strongly affected by changes in the operating  3.14 ± 0.01 3.13 ± 0.07 3.14 ± 0.04 Chemostat blank 3.07 ± 0.03 3.14 ± 0.04 3.32 ± 0.02 parameters. Configuration of the system has been shown to have the strongest effect on the community diversity, above loading rates and medium composition (Pholchan et al., 2010). 4CA has been reported to show inconsistent behaviour in ready biodegradability testing with variable outcomes and in many cases no degradation at all (Comber and Holt, 2010). Exposure to 4CA has been shown to promote its biodegradation and to greatly reduce the lag phase (Tor€ ang and Nyholm, 2005). Using batch and semicontinuous pre-exposure procedure (SCEP) to expose surface water microbial communities to low concentrations (0.1e100 mg/l) of 4CA resulted in a reduction of the lag-phase from 88 d, for the original river water, to 9 d after a 5 week exposure period (Tor€ ang and Nyholm, 2005). However, the semi-continuous system allowed the biomass to be retained, while in a chemostat, microorganisms with a generation time lower than the dilution rate will be washed out from the system. Growth in chemostats affected the taxonomic and physiological diversity of the communities, without there being a clear distinction in biodegradation capacity between the chemostats. Differential abundance analysis indicates that the loss of microorganism from the genera Pseudomonas (Loidl et al., 1990;Nitisakulkan et al., 2014;Surovtseva et al., 1985;Vangnai and Petchkroh, 2007) and Acinetobacter (Hongsawat and Vangnai, 2011;Vangnai and Petchkroh, 2007), which have been identified in the literature as 4CA degraders, might have been the reason for changes in the biodegradation capacity of the community. Interspecies competition for resources led to selection of a community able to degrade the easily biodegradable carbon sources, such as sodium acetate and peptides and to the low 4CA concentration (Ferenci, 2007).
NMP is used as a starting material for the synthesis of pharmaceuticals (Dorokhova et al., 1974) and the fate of this compound upon emissions to the environment is unknown. Biodegradation of NMP has not been reported in the scientific literature, as this compound is rarely studied in environmental science. No data are available on the occurrence and fate of NMP in wastewater and the environment. The parent compound piperazine has been reported to have a lag phase for biodegradation up to 50 d (Cai et al., 2013). The exposed community could adapt to NMP and completely degrade it, however, factors influencing and regulating the adaptation to this compound remain unclear and need further investigation. The presence of the chemical itself, or a structural analogue, is a crucial condition to induce community adaptation (Alidina et al., 2014;Birch et al., 2017;Oh et al., 2013;Pfaender et al., 1985;Thouand et al., 1996). To fully understand the mechanisms behind the adaptation to NMP, we should investigate the microbial interactions within the community, the catabolic potential of the whole community, and the influence of the environmental factors on the community. A population found in NMP a stable substrate, leading to the development of a degrading population for this ecological niche, and therefore to an adaptation at the community level (van der Meer, 2006). However, these changes observed in the taxonomic diversity could not be linked to the changes of the functional diversity. Functional diversity analysis only points out that the cultivated communities are different from the original activated sludge used as inoculum, as demonstrated by the taxonomic diversity analyses. Abundance in activated sludge of NMPdegrading microorganisms must have been too low to efficiently remove NMP under the stringent conditions of the RBT. The ECHA database on the safety data sheets (SDS) of NMP also point out this hypothesis to explain the absence of biodegradation by activated sludge, while an adapted community could quickly remove the compound (ECHA, 2013). The supplier of the SDS observed that NMP was degraded within a week in an enhanced biodegradation test (OECD 301) by an inoculum exposed for four weeks to 10 and 100 mg/l in a semi continuous activated sludge system (SCAS).
Furthermore, similar results were obtained by using sessile microorganisms collected from glass beads during flow through of river water containing mg/l of NMP (ECHA, 2013). Hence, this paper confirms the previously observed mechanism that was not published in the scientific literature.
Differences in response between the communities exposed to NMP and 4CA show that biodegradation capacity in an open system is the result of a fine balance between the increase of competent species (i.e, adaptation) and loss of these strains. The presence of the test compound, the relative abundance of key degraders, interand intra-species competition and the presence of additional energy sources are some of the different parameters that will influence the capacity of a community to biodegrade substances. 4CA is a well-known environmental pollutant (Zhu et al., 2011), to which microbial communities may already been exposed and subsequently adapted. The environmental concentrations of NMP is unknown but its use as intermediate (ECHA, 2013) suggest that this compound is not prevalent in wastewater, which indicates that microbial communities might not be naturally exposed.
Our results confirm the importance of characterising inocula and their exposure history before using them in biodegradation testing, as this can affect the outcome of the test (Birch et al., 2017;Kim et al., 2017;V azquez-Rodríguez et al., 2011). Including information about the exposure history of the inoculum would enhance the OECD biodegradation guidelines, however, it is challenging to obtain this type of information. Long term exposure in chemostats is one method, with its advantages and disadvantages, that can be used to investigate microbial adaptation and its implication for ready biodegradability testing. Chemostats offer the possibility to adapt a community in a controlled and reproducible environment. However, as described in this study, they can lead to a loss of biodegradation capacity (Poursat et al., 2019b), as they cannot fully mimic environmental conditions, and they can be expensive and difficult to handle in laboratory. Moreover, long term exposure may lead to the development of a highly artificial community, that does not represent the environment. However, it should be noted that results of the regular RBTs, without previous exposure, did not express a high results variability or false negatives, as sometime reported in the literature (Dick et al., 2016;Federle et al., 1997), hence, use of an enhanced test did not reduce the expected test variability or number of false negatives For these reasons, we do not advise to use highly adapted communities from chemostats as standard inocula in ready biodegradation test, as more investigation should be conducted before one can truly understand and therefore include adaptation in the tests. However, other systems with biomass retention may avoid some of these disadvantages and therefore be more suitable. Short-term pre-exposure, which is not allowed under the ECHA guidelines, seems to be a more reliable technique to produce a standardized inoculum for the tests, and deserve to be investigated. As proposed in a previously published work (Poursat et al., 2019a), we would recommend to perform preexposure, as an add-on, in parallel with the standard tests. In one test, the inoculum would be pre-exposed in a closed system, while in the second test, the inoculum will not be pre-exposed, and only processed according to the guideline. However, characterising the inoculum, by using multi-omics approaches (Kowalczyk et al., 2015), is a first step to reduce test outcome variations and false negative results.

Conclusions
This study highlighted the potential of chemostat systems to induce microbial adaptation to specific chemicals. The present study is another example of exposure leading to adaptation and is a first step to understand the community mechanisms and dynamics involved in this process, in order to implement the adaptation of microbial communities to the ready biodegradability testing guidelines. The two test chemicals, 4-chloroaniline and N-methylpiperazine, were respectively completely eliminated and recalcitrant to biodegradation in the OECD 310 biodegradation test, using activated sludge as inoculum. However, after one month of exposure, the community adapted and was able to completely remove n-methylpiperazine from the batch degradation test within 10 days. On the other hand, cultivation in chemostat in the presence of 4-chloroaniline, led to a partial loss of biodegradation capacity for this compound after 3 months of exposure. Results of this study confirmed the importance of characterizing and standardizing the inocula used in ready biodegradability testing. Indeed, exposure to a chemical, in either the laboratory or in the environment, clearly changes the community behaviour in biodegradation testing. Investigating the implications of the chemical exposure on the test outcomes, will ultimately lead to a better prediction of the environmental fate of organic chemicals. Estimated time to reach 10%, 50% and 90% of degradation (DT10, DT50 and DT90) of N-methylpiperazine by each inoculum. Degradation times were estimated using a three-parameters log-logistic model.


extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Psychology

text:
Brain weight predicts yawn duration across domesticated dog breeds

Abstract Previous research shows that yawning enhances intracranial circulation and regulates brain temperature. Consistent with these functional outcomes, yawn duration correlates positively with interspecies variation in brain weight across mammals, with robust relationships documented at both the taxonomic rank of class and the more restricted scale of family (e.g., Felidae). This study provides the first investigation into whether differences in brain weight within a single species, domesticated dogs Canis lupus familiaris, can predict intraspecific variation in yawn duration. Measures of yawn duration were obtained from public videos available online and then paired with previously published brain and body weight data of different dog breeds. The final sample consisted of 272 yawns from 198 dogs across 23 breeds. Consistent with recent studies, we find robust correlations between yawn duration and brain weight across breeds. Moreover, these correlations remain significant after controlling for differences in body weight across breeds. These findings replicate and extend upon past work in this area and provide further support that yawns evolved to serve an important and large-scale neurophysiologic function.

Yawning is characterized by a powerful gaping of the jaw with inspiration, followed by a brief period of peak muscle contraction and a passive closure of the jaw with shorter expiration (Barbizet 1958). While yawning appears to be a stereotyped action pattern across vertebrates (Provine 1986), recent research has identified distinct variations in the expression of this behavior. For example, some primate species display different yawn types that vary in intensity and morphology, which may correspond with transitions in daily repertoires, activity patterns, and internal states (Palagi et al. 2009;Vick and Paukner 2010;Leone et al. 2014). This work supports research suggesting yawns might be multifunctional across different contexts (Baenninger 1997;Provine 2005;Gallup 2011).
Studying differences in the expression of yawning across species could also help elucidate the neurophysiological significance of this highly conserved response. Based on hypotheses stating that the motor action pattern of yawns serves to enhance intracranial circulation (Walusinski 2014) and regulate brain temperature (Gallup andGallup 2007, 2008;Ramirez et al. 2019), the duration or magnitude of this response should correspond to the degree of neurophysiological change. In particular, the brain cooling hypothesis states that the gaping of the jaw and deep inhalation of air during yawning functions to promote thermal homeostasis by altering the rate and temperature of arterial blood traveling to the skull (Gallup and Hack 2011). Support for this hypothesis comes from measuring changes in brain temperature before and after yawns (Shoup-Knox et al. 2010;Shoup-Knox 2011) and naturalistic and experimental studies assessing the relationship between yawning and ambient temperature (e.g., Massen et al. 2014;Eldakar et al. 2015). Although the neural structures necessary for yawning appear to be located in the brain stem (Heusner 1946), recent comparative research shows that yawns produce large-scale neurophysiologic effects as evidenced by sizable cooling at the surface of the skull (Eguibar et al. 2017;. Therefore, according to this hypothesis, animals with larger and more complex brains should display longer and more powerful yawns to achieve these outcomes. In the first study to test this prediction, Gallup et al. (2016) collected yawn duration data across a diverse sample of mammals and linked these with published data on both brain weight and cortical neuron number. The results of this work revealed average yawn duration to be highly variable across 24 mammalian taxa (range: 0.8-6.5 s). Moreover, consistent with neurophysiologic hypotheses of yawning, interspecific differences in yawn duration were robustly correlated with both brain weight (correlation coefficient ¼ 0.911) and cortical neuron number (correlation coefficient ¼ 0.951) across species. Variation in yawn duration was not tied to particular features of body/jaw size, but was correlated positively with encephalization quotient across species.
As a direct follow-up to this work,  investigated whether the relationship between yawn duration and brain size could be identified within a single family of mammal species. Examining this link at a restricted taxonomic scale allowed for the determination of whether variability in yawn duration was related specifically to brain size or to other features pertaining to phylogenetic history and morphological diversity between species. In a sample of 58 animals across 7 wild cat species in the family Felidae, a similarly robust association (correlation coefficient ¼ 0.937) was observed between yawn duration and brain weight . This study also showed that body size was not a significant predictor of yawn duration. These combined findings support the view that the motor action pattern of yawning serves an important function tied to brain size.
Here, we extend upon this line of research by testing whether intraspecific variation in brain weight predicts differences in yawn duration within a single species. Domesticated dogs Canis lupus familiaris were chosen due to large differences in brain size across breeds and the availability of open access yawning data on the Internet. A report by Bronson (1979) that documented the average adult brain and body weights from over 2 dozen dog breeds was used as a basis for the inclusion of breeds within this study. This report also included a scaling function for estimating brain weight based on the average body weight for a given breed, which was recently validated and applied to a much larger sample of domesticated dogs (Horschler et al. 2019). The same methods used in Gallup et al. (2016 were applied to acquire representative yawn durations from each breed using public videos available online. We hypothesized that dog breeds with larger brains would display longer yawns, and that this relationship would remain significant after controlling for body size.

Materials and Methods
Using the list of 26 dog breeds compared in Bronson (1979), a researcher was instructed to find up to 12 adult dogs yawning from each breed by searching videos posted on YouTube and related websites. For each video identified, the researcher noted the time(s) at which the yawn(s) on the clip occurred and the duration of all yawning events using the operational definition provided by Barbizet (1958). Durations were recorded to the nearest 0.01 s using the stopwatch feature of a smartphone. Data collection occurred during Fall 2016, and yielded a total of 299 yawns from 213 dogs across the 26 breeds.
In Fall 2017, an independent rater was provided with the complete list of uniform resource locators (URLs) identified from the original search and instructed to score the duration of each yawn across the sample. In doing so, 4 videos were no longer available and 8 yawns were excluded because they were noted as either being too difficult to distinguish or there was no clear way to determine their beginning or end. In addition, links for 3 dogs were excluded because they were identified by both raters as juveniles. The interrater reliability for yawn duration was high for the remaining sample: intraclass correlation ¼ 0.87.
The average yawn duration for each breed was then calculated by summing all the respective durations from members of that group and then dividing by the number of yawners. To avoid pseudoreplication in cases where the same individual displayed multiple yawns in a given video, the average duration for this animal was used as a single data point when generating the overall average for that breed. The maximum yawn duration was also recorded from each breed. Following the same criterion set within Gallup et al. (2016, analyses were restricted to breeds with yawns from at least 3 individuals. The final sample included a total of 272 yawns from 198 dogs across 23 breeds. Individual yawning data from the following breeds were included in the sample: Beagle (11), Boston Terrier (12), Boxer (12), Bulldog (11), Chihuahua (9), Cocker Spaniel (3), Collie (5), Dachshund (12), Doberman Pinscher (9), German Shepherd (11), Golden Retriever (12), Great Dane (12), Labrador Retriever (12) (12), and Weimaraner (10). The Supplemental Material includes URLs to all videos used in the study. For the average brain and body weights across breeds, Bronson (1979) provided data based on an average sample of 31.69 dogs per breed (range: N ¼ 5-110). This report also included a scaling function for estimating brain weight based on the average body weight for a given breed. Horschler et al. (2019) recently validated this scaling function and applied it to a much larger sample of domesticated dogs (N ¼ 7,397 across 74 breeds). Therefore, where applicable in this study (19/23 breeds), these recent and more representative body weight and estimated brain weight data were used (see Table 1).
Based on output from Kolmogorov-Smirnov normality tests, Pearson correlations were used to assess the relationships between measures of yawn duration and brain and body weights across the 23 breeds. Separate partial correlations were then run controlling for body and brain weight, respectively. Despite our a priori hypotheses, all analyses included more conservative 2-tailed tests with the alpha set at 0.05. Moreover, Benjamini and Hochberg corrections were applied for multiple comparisons, and thus adjusted P-values are provided.

Results
The average yawn duration for all dogs in the sample was 2.04 6 0.59 s, with the breed average ranging from 1.43-2.83 s (Table 1). Consistent with previous research, there were strong positive correlations between measures of yawn duration and brain weight across breeds (average duration: r 23 ¼ 0.819, P < 0.001; maximum duration: r 23 ¼ 0.703, P < 0.001; Figure 1). Breed brain and body weights were highly correlated across the sample (r 23 ¼ 0.957, P < 0.001), and thus yawn duration was also positively correlated with body weight (average duration: r 23 ¼ 0.749, P < 0.001; maximum duration: r 23 ¼ 0.649, P ¼ 0.001). Partial correlations showed that average yawn duration remained significantly correlated with brain weight after controlling for body weight (average duration: r 20 ¼ 0.535, P ¼ 0.016; maximum duration: r 20 ¼ 0.374, P ¼ 0.111). However, yawn duration was not correlated with body weight when controlling for brain weight (average duration: r 20 ¼ À0.222, P ¼ 0.361; maximum duration: r 20 ¼ À0.125, P ¼ 0.579).

Discussion
Recent research has shown that mammalian brain weight is a robust predictor of yawn duration across different taxonomic scales (Gallup et al. 2016, supporting past research suggesting that this ubiquitous motor action pattern serves an important neurological function. This study reveals that yawn duration is also tied to intraspecific variation in brain size among domesticated dogs. In particular, we show a wide range in the duration of this action pattern across breeds, with larger brained dogs tending to have longer yawns. Notably, over two-thirds of the variability in the average yawn duration across breeds can be explained by brain weight differences between these populations. Partial correlations show that this association is not due to differences in body size across breeds, but rather differences in yawning appear to be specifically associated with variation in the brain size of these animals. It remains unknown exactly which neurological features are driving this association, but past research suggests that cortical neuron number is the strongest predictor of interspecies differences in yawn duration (Gallup et al. 2016). Although research on neuron numbers in dogs is currently limited, in a comparison of 2 animals Jardim-Messeder et al. (2017) provide some evidence that larger breeds have more neurons (i.e., Golden Retriever: 627 million; unspecified smaller breed: 429 million).
While we show here that intraspecific variation in brain size appears to be a robust predictor of yawn duration in domesticated dogs, it is unknown how these measures would correlate within a given breed or in a different wild type species. Due to artificial selection, the variation in brain and body weights across dog breeds exceeds that of any other species, making this study a distinctive and atypical within-species comparison. Based on the findings of recent work uncovering breed differences in neuroanatomical organization (Hecht et al. 2019), one could argue that assessing neurological differences across dog breeds becomes similar to interspecies comparisons. Thus, until further research is conducted, we are hesitant to conclude that yawn duration would be predictive of intraspecific variation in brain size in other animals. One intriguing species to investigate would be humans, since there is sizable variation in the Research on humans or another species could also provide a more precise analysis of this relationship, by linking brain volumes acquired through magnetic resonance imaging to electromyography from the masseter and submental muscles to measure yawn duration and intensity (see Ertekin et al. 2015).
It is important to acknowledge limitations to this study with regard to the methods of data collection. Although visual media can provide a powerful tool for investigating and developing insights into animal behavior (Nelson and Fijn 2013), and research utilizing YouTube for this purpose has increased within the past decade (see Measey et al. 2019), one potential issue with human posted videos of this nature is nonrandom sampling (Burn 2014). While we do not view this as a major concern given the strength of the correlations both for average and maximum yawn duration, naturalistic observations could be conducted to confirm the true variability in spontaneous yawn duration for a given breed. The inability to determine the contextual triggers of yawning from online videos also presents a limitation. While there is consistency in that most videos occurred within the context of close dog-human interaction, other potentially important factors known to alter the morphology of yawning in other species (i.e., Leone et al. 2014;Palagi et al. 2009;Vick and Paukner 2010) could not be evaluated. To date, however, no study has reported differential yawn types for dogs.
Overall, this study replicates and extends upon previous results showing a robust link between yawn duration and brain weight, providing further support for predictions derived from the brain cooling hypothesis. However, it is important to acknowledge that the current findings do not necessarily exclude other neurophysiological hypotheses of yawning. We propose that, in addition to examining the connection between yawn duration and brain size in other classes of vertebrates (e.g., avian species), future research could investigate how differences in the duration of yawning events correlate with immediate changes in neurophysiological measures and cognitive processing. Work of this nature could help elucidate the evolutionary significance of this ubiquitous response.