Massively parallel metabarcoding of droplet co-cultures

Massively parallel metabarcoding of droplet co-cultures | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Massively parallel metabarcoding of droplet co-cultures View ORCID Profile James Y. Tan , Jessica D. Li , Auden Bahr , View ORCID Profile Xiaoxia Nina Lin doi: https://doi.org/10.1101/2025.01.14.633022 James Y. Tan 1 Department of Chemical Engineering, University of Michigan , Ann Arbor, MI Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for James Y. Tan Jessica D. Li 2 Department of Microbiology and Immunology, University of Michigan , Ann Arbor, MI Find this author on Google Scholar Find this author on PubMed Search for this author on this site Auden Bahr 3 Gilbert S. Omenn Department of Computational Medicine and Bioinformatics, University of Michigan , Ann Arbor, MI Find this author on Google Scholar Find this author on PubMed Search for this author on this site Xiaoxia Nina Lin 1 Department of Chemical Engineering, University of Michigan , Ann Arbor, MI 2 Department of Microbiology and Immunology, University of Michigan , Ann Arbor, MI 4 Department of Biomedical Engineering, University of Michigan , Ann Arbor, MI Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Xiaoxia Nina Lin For correspondence: ninalin{at}umich.edu Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract A predictive understanding of microbiomes is necessary to engineer them for human and environmental health. While ‘omics data provide insight into patterns of microbial diversity and function, we lack understanding on how microbial interactions contribute to community function and emergent ecological properties. Synthetic ecology approaches seek to address this by constructing and interrogating defined co-cultures of cultured representatives. This provides phenotypic observations on how interactions change in different contexts and contribute to overall community function. However, this cannot be applied to systems where relevant isolates are not available. Here we developed Cocoa-seq (combinatorial co-cultivation and amplicon sequencing), a microfluidic workflow utilizing the high throughput nature of nanoliter-scale, water-in-oil droplets to generate and profile co-cultures generated by the stochastic co-encapsulation of cells from samples, circumventing laborious axenic isolation. The workflow multiplexes over a thousand 16S amplicon libraries from droplet co-cultures into one, with potential for even higher scalability, to enable metabarcoding of each droplet community. To validate Cocoa-seq, we compared community distributions generated from Cocoa-seq with two benchmark communities, a model synthetic bi-culture and synthetic mock communities of various compositions, to those from other measurements or expectations. Benchmarking showed qualitative agreement with fluorescence microscopy and quantitatively similar compositions and beta-diversity patterns compared to expectations. We also introduced discrete spike-in molecular standards to estimate absolute abundance, but uneven amplification due to stochastic bias prevented accurate quantification. We discuss Cocoa-seq’s limitations and recommend how Cocoa-seq can be employed to identify interactions or functionally critical consortia. Introduction The contribution of microbiomes to host-associated ( 1 – 3 ) and environmental health ( 4 – 6 ) is becoming clearer, and ways to predict and manage them are increasingly necessary. While environmental factors certainly influence overall community composition, microbe-microbe interactions are hypothesized to be the major driving factors, especially at more resolved scales ( 7 ). However, identifying and measuring these interactions to represent them in theoretical models is difficult ( 8 ). The accessibility of culture-independent sequencing has generated large 16S metabarcoding and ‘omics datasets to infer interactions computationally, such as through statistical correlations ( 9 , 10 ), genome-scale models ( 11 , 12 ), or multi-omics ( 13 ). While these approaches are scalable, they need empirical data to validate or improve them, which remains challenging to obtain. Synthetic ecology, which involves studying a microbiome using simplified laboratory co-cultures composed of relevant isolates, is a promising approach for these high-resolution phenotypic and functional data ( 14 ). Synthetic ecology has enabled strides in the design and study of microbial communities and their functions; it has been used to study starch degradation in a six-member soil community ( 15 ), promote growth of a probiotic strain in different combinations of a 19-member soil community ( 16 ), and optimize butyrate production in combinations of a 25-member human gut consortia ( 17 ). Since the total number of community combinations is 2 n − 1 for a community of n members, these studies use multiwell plates, microfluidics, and a liquid handling robot, respectively, to overcome the challenge of combinatorial scalability. A recent meta-study utilizing these datasets generated statistical models which predicted community function from member presence or absence ( 18 ), which is a major stride towards rational microbial community design. While this reductionism is certainly advantageous, the major technical barriers are the bias and labor of laboratory isolation and the low scalability of co-culturing apparatuses compared to what is necessary for a comprehensive investigation of the expansive combinatorial compositional space of a complex microbiome. Microfluidic droplets provide an avenue to employ the reductionism of synthetic ecology while embracing natural complexity by enabling the combinatorial co-cultivation of native cells from a complex sample. Specifically, droplets as uniform water-in-oil emulsions co-encapsulate and co-cultivate random combinations of native cells from samples in localized sub-nanoliter co-cultures ( 19 – 21 ). These droplet co-cultures enable massively parallel co-cultivation which circumvent laborious strain isolation. However, while droplet co-culture generation and incubation is highly scalable, the current bottleneck is the applicability and throughput of droplet co-culture analysis. The most high-throughput droplet-based workflows to date utilize fluorescence ( 21 ) or droplet-color coding ( 16 ) to identify which members are in droplets, which require genetic manipulation or strain isolation, respectively, limiting investigation to cultured representatives. Another approach applies amplicon and shot-gun metagenomic sequencing to droplet co-cultures derived from the human gut microbiome ( 20 ), but is limited to two dozen co-cultures. To improve the scalability of droplet co-culture analysis for combinatorial co-cultivation, we developed Cocoa-seq ( co mbinatorial co culturing and a mplicon seq uencing), a droplet-based workflow for the massively parallel 16S amplicon sequencing of thousands of droplet co-cultures that can be applied to both synthetic consortia and natural microbiomes. Inspired by droplet-enabled workflows for single cell sequencing ( 22 , 23 ), single cell transcriptomics ( 24 – 27 ), and 16S amplicon-based spatial mapping ( 28 ), we utilize hydrogel beads to enable multi-step processing of the droplet co-cultures and barcoding for multiplexed sequencing. To validate Cocoa-seq, we characterized the quality of droplet-derived amplicon libraries and benchmarked the workflow’s ability to accurately measure community composition with a model synthetic co-culture and various mock communities. We also describe attempts to improve Cocoa-seq by incorporating unique molecule standards for absolute quantification and how PCR stochastic bias at the discrete template level prevents this. Finally, we discuss recommendations for the further improvement and the application of Cocoa-seq to study interactions in natural communities. Materials and Methods Lithography SU-8 mastermolds with microfluidic device features were generated with standard SU-8 photolithography procedures in a cleanroom environment. Both PDMS microfluidic devices for the generation of agarose beads and high-throughput sample multiplexing are detailed in Fig S1 and were generated with standard soft lithography methods. Full details are in the Supplementary Information. Microbial cultures For benchmarking, we used six microbial cultures: E. coli K12 BW25113 (ΔmetA mCherry), B. subtilis 168 (Pveg-GFP), B. subtilis 168 (Phyperspank-GFP), Pseudomonas putida KT2440, B. thetaiotaomicron VPI-5482, and L. crispatus ATCC 33820. Each overnight cell culture was washed twice, resuspended in PBS, and quantified by counting with a disposable haemocytomer under phase contrast microscopy. The E. coli and B. subtilis 168 (Phyperspank-GFP) strains were used together as an obligately-syntrophic, two-species model co-culture similar to the one used in Hsu et al. ( 21 ). When co-cultivated in droplets, the E. coli/ B.subtilis co-culture was grown in M9. Three 4-member mock communities were comprised of E. coli , B. subtilis 168 (Pveg-GFP), P. putida , and B. thetaiotaomicron in different ratios. Details on media and culture conditions are in the Supplementary Information. Encapsulation of cells in agarose beads Droplet generation was performed in a modified large oven incubator (VWR Scientific 1535) to maintain a temperature between 37-40°C to keep low-melting agarose suspensions fluid (Fig S2). Cell suspensions were mixed with agarose solution to obtain the desired cell concentrations and 1.5% agarose solution. Syringe pumps were used to flow oil (HFE7500 Novec engineering oil with 2% fluorosurfactant) and aqueous suspensions through PFTE tubing into the microfluidic device at rates of 3 µL/min each. Droplet emulsions were collected in 1.5 mL microcentrifuge tubes. For co-cultivation, the emulsions were transferred to a 37°C incubator for 20-24 hours. Agarose in droplets was set by placing the emulsions in ice for 10-20 minutes. To collect agarose beads from droplets, we washed the emulsion with PBS-wash buffer and 20% (v/v) perfluorooctanol in HFE7500 oil, washed with 1% (v/v) Span-80 in hexane, and performed final washes in PBS-wash. Full details are in the Supplementary Information. Microscope imaging of droplets For quantification of E. coli/ B.subtilis droplet co-culture community abundances, fluorescence images of hundreds of droplets were taken after co-cultivation. Images were processed in ImageJ, and custom MATLAB scripts were used to quantify total fluorescence from both strains in individual droplets to estimate relative abundance. Full details are in the Supplementary Information. Cell lysis in agarose beads Agarose beads were washed with 1 mL 10 mM Tris-HCl three times, suspended in lysis buffer (1 mM DTT, 1 mM Tris-HCl pH 8.0, 2.5 mM EDTA, 100 mM NaCl, 0.8% lysozyme), and incubated in a shaker at 37°C overnight. Beads were then washed in 10 mM Tris-HCl pH 8.0, suspended in digestion buffer (30 mM Tris-HCl pH 8.0, 10 mM EDTA, 0.8% Triton X-100 (v/v), 0.5% SDS, 1 ug/µL proteinase K), and incubated for 30 minutes at 50°C. To deactivate proteinase, the beads were washed in 10 mM Tris-HCl pH 8.0, 10 mM EDTA, 0.1% Tween-20 (v/v), 5 mM phenylmethylsulfonyl fluoride (PMSF) (Sigma-Aldrich 93482), washed in TET buffer (10 mM Tris-HCl pH 8.0, 10 mM EDTA, 0.1% Tween-20 (v/v)) five times, and stored at 4°C until use. To check that genomic DNA was immobilized in the microgels after cell lysis, a subset of agarose beads was stained with SYBR Green I and visualized with fluorescence microscopy. Full details are in the Supplementary Information. Barcoded hydrogel bead processing Barcode hydrogel beads were purchased as a unit of 1 million suspended in TET buffer (RAN Biotechnology, CustomSeqReady-1M, InDrop variation, ca. 147k). 16S 515f primers were extended onto the barcode beads by complementary extension oligos based off a protocol from Zilionis et al. ( 25 ). Beads were washed extensively with multiple buffers to stop extension and remove extension oligos. Any unelongated oligonucleotides on the bead were removed by ExoI clean-up (Thermo Fisher Scientific, EN0581). Final beads are washed and stored at 4°C. After the extension protocol, barcode beads were checked to ensure proper 515f extension through the hybridization of complementary fluorescent oligo probes and visualization with fluorescence microscopy. Full details are in the Supplementary Information. 16S multiplexing 515f-extended barcode beads, gDNA-agarose gels at a 50% (v/v) suspension, and 2x PCR mastermix solution containing 816r reverse primer, SuperFi II buffer, dNTPs, Platinum SuperFi II DNA polymerase, Pluronic F-68, and bovine serum albumin (BSA) were flowed into the microfluidic device with syringe pumps at rates of 0.45-0.5 µL/min, 0.8-1 µL/min, and 3 µL/min, respectively. Oil (Droplet Generation Oil for EvaGreen, Biorad) was flowed at 5 µL/min. Droplets of around 150 µm diameter were generated. For each sample, 20-40 µL of droplets and 20 µL of fresh Biorad oil were added into PCR tubes and exposed to 365-nm UV-light for 10 min utes on ice. 50 µL of mineral oil was placed on top of the droplets to prevent evaporation, and droplets were thermocycled (95 °C for 2 min; 98 °C for 30 sec; and 30 cycles of 98 °C for 10 seconds, 60 °C for 10 seconds, and 72 °C for 30 sec) with a 2 °C/s ramp rate between all steps with no heated lid. Droplets underwent clean-up, involving droplet breaking with 20% perfluorooctanol in oil, primer degradation with ExoI, and Ampure bead purification. Full details for droplet multiplexing are in the Supplementary Information. Unique molecular identifier (UMI) benchmarking To quantify the amplification of individual 16S template molecules, we designed an oligonucleotide composed of a modified 16S V4 sequence from Thermus thermophilus ATCC 33923 with seven nucleotides directly after the forward primer region replaced with seven degenerate nucleotides for use as unique molecular identifiers. The 16S standard was purchased from GenScript as a ssDNA stock, resuspended in molecular-grade water, and quantified with an in-house digital droplet polymerase chain reaction (ddPCR) to determine the concentration of functional 16S standard molecules. The stock was diluted and co-flowed into droplets in the droplet barcoding device and workflow for a desired initial λ of 50 molecules/droplet. Full details for oligonucleotide design and in-house ddPCR protocol are in the Supplementary Information. Illumina library preparation and sequencing Illumina adaptors and indices were attached to each multiplexed library via PCR with NEBNext Q5 HotStart HiFi polymerase. The thermocycler program was as follows: 98 °C for 30 seconds; 5-7 cycles of 98 °C for 10 seconds, 68 °C for 20 seconds, and 65 °C for 30 seconds; 65 °C for 2 minutes; and a 12 °C hold. PCR product was verified with gel electrophoresis and purified with Ampure XP bead clean-up. Libraries were sent to the University of Michigan Advanced Genomics Core for quality control with Agilent TapeStation and sequencing with the Illumina NextSeq 2000 with P1 300 cycle chemistry. Details for library preparation PCR and sequencing details are in the Supplementary Information. Bioinformatic processing A custom demultiplexing Python script (run on v.3.10.4) based on a previous study ( 25 ) was modified and used to demultiplex reads and trim barcodes from the forward reads. Because of the large number of barcode families (groups of reads with the same barcode) with a very small number of total reads, we used the expected number of barcodes (Note S1) to filter on the barcode families with the highest number of reads. Using vsearch (v.2.13.7-linux-x86-64), operational taxonomic units (OTUs) were constructed for reads within a 95% similarity threshold to expected 16S v4 sequences for known isolates with the usearch_global command ( 29 ). To run sample rarefaction and beta-diversity analyses, we used mothur (v.1.44.1) ( 30 ). To compare sequencing reads to cell abundances, we normalized abundance by 16S copy number for each strain, retrieved from rrnDB ( 31 ). Sequence read processing scripts were run on the Great Lakes HPC Cluster (University of Michigan). Stochastic amplification simulation To simulate the effect of stochastic bias for a low number of initial molecules, we programmed a simple Python script with a model in which individual template molecules were represented uniquely as different characters in an array. The probability of success of amplification in one round of PCR was represented by the Bernoulli distribution in which success of a molecule is given by η amp (amplification efficiency) and failure is given by 1 – η amp . If a molecule was replicated, its unique character string was added to the array. This was propagated for a designated number of cycles for each molecule. At the end, to simulate sequencing, we subsampled the number of molecules to 1000 molecules and the final abundances of molecules were compared to the initial. Results Development of a workflow for droplet-enabled high-throughput co-culturing and 16S multiplexing The Cocoa-seq workflow utilizes the stochasticity and throughput of droplet co-encapsulation to both generate and multiplex thousands of random co-cultures in 16S amplicon sequencing libraries ( Fig. 1 ). We enable the high-throughput, multi-step processing of individual droplet co-cultures with the compartmentalization of individual co-cultures in individual hydrogel beads. Download figure Open in new tab Fig. 1. Experimental overview of Cocoa-seq for highly multiplexed 16S amplicon sequencing of droplet co-cultures. (a) Schematic for generation of co-culture agarose beads which hold cells and gDNA of co-cultures generated from random droplet co-encapsulation. (b) Microscope images of agarose beads with gDNA from lysed droplet co-culture cells generated at the end of (a) (left is bright field image, right is SYBR stain overlay). Scale bar is 80 µm. (c) Schematic for multiplexed 16S library generation using co-culture beads from (a) and barcoded hydrogel beads (BHBs) functionalized with bead-unique barcode oligonucleotides. (d) Droplets generated from the co-encapsulation of co-culture beads and BHBs. Polyacrylamide beads are darker in outline and slightly smaller than agarose beads. Scale bar is 150 µm. (e) Droplets after bead dissolution and droplet PCR. The agarose beads have melted while polyacrylamide beads remain. Scale bar is 150 µm. Droplets are not the same subset as those in d. (f) After processing droplets and sequencing, bioinformatic demultiplexing of introduced barcodes from BHBs enables identification of reads from each individual droplet. The first half of the workflow involves generating these co-culture beads ( Fig. 1a ), and the second half individually barcodes them for highly multiplexed sequencing libraries ( Fig. 1b ). In the first half of the workflow, cells are co-encapsulated and co-cultivated in droplets, immobilized in gelled agarose beads, and lysed to expose genomic DNA (gDNA) in beads for downstream amplicon sequencing library preparation ( Fig. 1a ). First, initial cells are suspended in growth media with molten agarose and randomly encapsulated following Poisson statistics. The average number of encapsulated cells per droplet (λ) is determined by the concentration of the cellular suspension and the size of the droplet generated ( 32 ), which is set by flow rates and microfluidic channel dimensions. Generated droplet co-cultures are co-cultivated under appropriate temperatures and other conditions. After co-cultivation, co-cultures are immobilized in hydrogel beads by lowering the temperature and setting the agarose. The droplet emulsion is then broken, and co-culture beads are washed to remove residual oil and resuspended in an aqueous suspension. The co-culture beads are subject to enzymatic and chemical cellular lysis, which leaves high molecular weight gDNA in the agarose matrix ( Fig. 1b ). This was an effective cell lysis method for most strains tested, except one, suggesting it may be insufficient for cell types that require more intense physical lysis protocols (Fig. S3). The second half of the workflow prepares the previously processed gDNA-containing co-culture beads by pairing them with barcoded hydrogel beads (BHBs) in another droplet co-encapsulation step and barcoding 16S bacterial marker genes with droplet PCR ( Fig. 1c ). These BHBs are polyacrylamide hydrogel beads functionalized with oligonucleotide sequences, which act as 16S V4 primers with bead-specific unique barcode overhangs in downstream PCR. A microfluidic device is used to co-encapsulate co-culture beads and barcoded hydrogel beads (BHBs) with PCR reagents. This device accommodates close-packed flows of two types of hydrogel beads with different mechanical properties to enable bead pairing (Note S2) and is demonstrated in the Supplementary Video. Pairing of these two types of beads during co-encapsulation is only partially efficient due to microfluidic constraints ( Fig. 1d , Fig. S4), but a large fraction of droplets contains both. After bead co-encapsulation in droplets, barcoded oligos are photocleaved from hydrogel beads by 365-nm UV photoexposure, and agarose beads are melted to release gDNA. Droplet PCR then generates droplet-specific barcoded 16S amplicons. The droplets remain thermostable ( Fig. 1e ) until intentionally broken. Illumina adaptors and indices are attached to the multiplexed amplicons by bulk PCR, and library amplicon products are confirmed by size (Fig. S5) before sequencing. After sequencing, demultiplexing of the barcoded 16S amplicon reads by barcodes enables the profiling of individual droplet co-cultures ( Fig. 1f ). The BHBs which enable co-culture library multiplexing for Cocoa-seq adapts designs from a previous single-cell workflow ( 25 ) and requires modification when purchased from vendors ( Fig. 2 ). The oligonucleotides functionalized onto the BHBs have three major components: a photocleavable spacer, a dual-index barcode region, and an extension adaptor. First, the photocleavable spacer is sensitive to a specific UV-wavelength; when cleaved, it detaches the oligo from the polyacrylamide bead. Secondly, the dual-index barcode scheme utilizes split-and-pool combinatorial synthesis ( 25 ), which enables rapid generation of hundreds of thousands of barcodes. Lastly, Cocoa-seq utilizes the extension adaptor to incorporate 515f 16S V4 forward primers onto BHB oligos by isothermal polymerase, and the attachment of the 515f primer sequence is confirmed and quantified with probes (Fig. S6). All three components allow for these oligos to act as barcoded forward primers for the bacterial 16S marker gene along with the appropriate reverse primer in PCR. The exact oligonucleotide sequence and molecular workflow for multiplexing with BHBs is specified in Fig. S7. BHBs can also be generated in-house according to relevant protocols ( 25 ) with appropriate modifications, such as replacing polyT regions originally intended for single cell eukaryotic RNA sequencing. Download figure Open in new tab Fig. 2. Molecular biology of barcoded hydrogel beads (BHBs) functionalized with barcode oligonucleotides used for multiplexing. Oligonucleotides provide unique bead-specific barcodes through combinatorial dual indices. After purchase, oligonucleotides on BHBs are extended in bulk with the 515f primer for amplifying the 16S bacterial universal marker gene. Full oligonucleotides, when cleaved off beads in droplets with a specific wavelength of UV light, act as primers during droplet PCR, resulting in barcoded 16S marker gene amplicons. Benchmarking shows barcode filtering is required to remove sequence noise from droplet libraries To assess the quality of the sequencing libraries generated by Cocoa-seq, we benchmarked with a model bi-culture and a set of synthetic mock communities. To demonstrate Cocoa-seq’s ability to accommodate and assess co-cultivation, we utilized an obligately syntrophic two-species co-culture, comprising of an engineered E. coli methionine auxotroph (mCherry-labeled) and a native B. subtilis tryptophan auxotroph (GFP-labeled). We generated and co-cultivated droplet co-cultures with an initial λ of 5 cells per each strain per droplet to provide both strains in each droplet for co-cultivation, which after overnight incubation showed increased total growth ( Fig. 3a ). Additionally, to demonstrate Cocoa-seq’s ability to accurately represent more complex co-cultures across compositional space, we created mock communities comprised of E. coli , B. subtilis , P. putida , and B. thetaiotaomicron present in different rank abundances and generated droplet libraries by randomly encapsulating community members in droplets. The three mock communities were: the “even” library comprised of a 1:1:1:1 λ ratio of each member; the “Ec-biased” library comprised of mostly E. coli and the other members in decreasing rank abundance; and the “Bt-biased” library with opposite rank abundance ( Fig. 3b ). During droplet generation, we generated droplets at a λ of 100 cells/droplet with these mock communities, which led to diverse community compositions within each droplet library due to Poisson statistics. Download figure Open in new tab Fig. 3. Bi-culture and mock community benchmarking libraries used to validate Cocoa-seq and distribution of barcode family (groups of reads sharing the same barcode) sizes from the bi-culture droplet library. (a) (Top) The model bi-culture is comprised of two cross-feeding species: B. subtilis 168, which is a GFP-labeled (green) native auxotroph for tryptophan, and E. coli BW25113, which is an mCherry-labeled (red) metA knockout mutant. (Bottom) Fluorescence microscopy with brightfield overlays before and after co-cultivation for 20 hours. Initial λ for each species was 5 cells/droplet. Imaged droplets are not the same between timepoints. (b) Three synthetic mock communities were comprised of four members ( E. coli , B. subtilis , P. putida , and B. thetaiotaomicron ) present in different rank abundances with the λ (average cells/droplet) specified in the labelled ratios. Lower rank abundance members occur less than once per droplet. (c) Histogram of the number of reads per barcode family from the bi-culture droplet library with a bin size of 10. The subset histogram shows a more resolved view (bin size of 1) of the smallest barcode families. (d) Rarefaction curve of the cumulative number of reads contributed from sampling barcoding families. We sequenced these benchmarking libraries generated with Cocoa-seq and assessed library quality throughout bioinformatic processing ( Table 1 ). After initial Illumina index demultiplexing (Table S1), individual droplet libraries were demultiplexed from each sequencing library, and sequences associated with ambiguous or low-quality barcodes were filtered out, resulting in 81-83% of total reads with meaningful and differentiable barcodes. There were more barcode families (groups of reads with the same barcode) than expected based on the number of droplets processed during droplet barcoding PCR, echoing results from a previous method utilizing a similar microfluidic workflow ( 28 ). For example, for the biculture library, which was sequenced the deepest, we observed 109 459 total barcode families, which is much higher than expected at around 1 500. Looking at the histogram of reads per barcode family in the bi-culture library shows most barcode families contain very few amplicons ( Fig. 3c ), and this is also observed in the other benchmarking libraries (Fig. S8). A cutoff for the number of reads per barcode family to be considered a genuine barcoding event was set conservatively to match the expected number of co-cultures barcoded (Note S2), keeping anywhere between 50%-80% of the filtered reads. We suspect that “noisy” reads below the barcode family threshold were generated from failed co-encapsulation events (e.g. droplets with only barcode beads), combinatorial barcode index hopping, barcode oligo synthesis error, or amplification error, and not from contamination since a majority of these reads were identical to reference sequences (between 91.3-96.4%). Overall, the total retention rate of reads through bioinformatic processing was between 40% and 65%. Additionally, the sequencing depth of our libraries was intensive enough to capture available droplet library, as demonstrated by a rarefaction curve of total reads from barcode families for the biculture library ( Fig. 3d ) as well as the other libraries (Fig. S8). View this table: View inline View popup Download powerpoint Table 1. Retention of reads through bioinformatic processing for each benchmark library. We removed reads that did not match expected barcode sequences, belonged to barcode families with too few total reads, and did not match expected 16S V4 sequences within 95 percent identity. The minimum reads/barcode threshold was set based on the number of expected barcodes from the volume of droplets thermocycled during Cocoa-seq (Note S1). we estimated 1500 communities derived from each mock community and 2500 for the bi-culture and made sure to set the threshold to keep less than the estimated number. Cocoa-seq shows qualitative agreement with fluorescence microscopy-based measurements in determining droplet bi-culture composition We used the model bi-culture to compare community composition distributions between Cocoa-seq and fluorescence microscopy. While fluorescence microscopy is semi-quantitative due to diffraction at the droplet-oil interface, fluorescence variation between cells and strains, focusing on a two-dimensional cross-section, and cell clumping, we accounted for the first two limitations by normalizing the fluorescence per cell with image processing ( Fig. 4a ). From this we calculated total fluorescence for both strains to estimate relative abundance across hundreds of droplets. The abundance distributions from both methods are qualitatively comparable: the community equilibrium state is heavily biased towards B. subtilis with significant variation between individual droplet co-culture outcomes ( Fig. 4b ) similar to other droplet co-culture studies ( 21 , 33 ) and is likely a result of cell-to-cell variability in viability and phenotype. However, the bias towards B. subtilis is dramatically pronounced in the sequencing dataset compared to fluorescence. Given the limitations in fluorescence measurements, we suspect that fluorescence estimates are most likely underestimating B. subtilis abundance. In particular the extensive growth and clumping of B. subtilis , visible as large green fluorescent aggregates, will not be quantitatively represented in a two-dimensional cross-section microscopy image. Download figure Open in new tab Fig. 4. Community composition distributions from fluorescence microscopy and from Cocoa-seq for droplet co-cultures of the syntrophic bi-culture qualitatively agree. (a) Example processing of the fluorescence image from Fig. 3a to detect droplets and normalize fluorescence signal per cell. (b) Histograms of B. subtilis relative abundance in individual droplets as measured by total fluorescence compared to Cocoa-seq. Each histogram represents a different set of 500 individual droplets from the same droplet pool for comparison. However, it is uncertain if Cocoa-seq introduces some degree of systemic sequencing bias, particularly for low-rank members. Cocoa-seq replicates expected beta-diversity patterns of droplet mock communities and captures low abundance members While the previous benchmarking demonstrated Cocoa-seq on droplet co-cultures, due to the lack of a ground truth comparison and the bi-culture’s simplicity, we sought to more rigorously evaluate the quantitative capability of Cocoa-seq with the previously-described 4-member synthetic mock communities. Because there is no cultivation, droplet community composition distributions are predictable from initial mock community composition before droplet encapsulation. To generate these, we simulated droplet co-encapsulation and sequencing in silico based on Poisson statistics and known 16S rRNA copy numbers for each strain. We first compared the distribution of droplet community compositions from both in silico expectations and Cocoa-seq with principal component analysis (PCA) ( Fig. 5a ). The droplet communities from Cocoa-seq cluster distinctly by mock communities of origin as expected but have a higher abundance of B. subtilis and a larger spread than the in silico libraries, which can be partially explained by cell clumping during droplet encapsulation (Fig. S9a) and the large cells of the B. subtilis strain used for mock communities, which may have contributed to overrepresentation in visual cell counting (Fig. S9b). Despite these deviations, the average relative abundances for each member for each mock community library are quite similar between Cocoa-seq and expectation, suggesting minimal overall bias (Table S2). Additionally, due to the rank abundance composition of the Ec-biased and Bt-biased mock communities, we assessed Cocoa-seq’s ability to detect rarer members which are only present in a subset of droplet communities. The prevalence of members in Cocoa-seq datasets across the 3 different mock communities, matches expectation well, with slight deviations from expectation which we suspect are a result of cell counting and clumping independent of Cocoa-seq capability ( Fig. 5b ). This highlights the ability of Cocoa-seq to capture rare members given the appropriate sequencing depth. Download figure Open in new tab Fig. 5. Community composition distributions from Cocoa-seq for droplet subsets of mock communities resemble those from expectations. (a) Principal component analysis (PCA) of community compositions in individual droplets derived from a mock community, as expected from simulations based on Poisson statistics and observations with Cocoa-seq (100 instances for each sample). The tetrahedron in 3D (with the third principal component coming out of the page) represents all possible compositions of the four-species ecosystem, with each vertex corresponding to a monoculture (gray circle). Point colors designate the mock community of origin. (b) Prevalence of each member across droplets, with circles representing observed prevalence from Cocoa-seq and crosses for expected prevalence from Poisson statistics. PCR introduces significant bias for single molecule resolution, preventing the use of spike-in molecular standards for absolute quantification We hypothesized that we could introduce spike-in standards during the droplet library generation step of Cocoa-seq to infer absolute abundances of droplet communities, enabling the differentiation of total biomass between droplet co-cultures ( 34 ) (Fig. S10). We proposed to include 16S V4 oligonucleotides with 7-mer random unique molecular identifiers (UMI) ( 35 ) (for a total of 4 7 possible unique molecules) at a λ of 50 molecules/droplet to be included during barcoding PCR and subsequent sequencing. The exact number of spike-in molecules for each droplet would be quantified for each droplet using the total observed UMIs and used to back calculate the absolute abundances of other members. However, when we quantified the amplification of each of these individual spike-in molecules in droplets after Cocoa-seq, there was a high degree of amplification efficiency between templates despite each template molecule being represented once in the droplet before PCR (Fig. S11a). Due to this significant variation, determining the total number of spike-in molecules per droplet was difficult due to the inability to determine a reasonable threshold to differentiate signal from noise. We hypothesized this was a result of stochastic bias at the single molecule level ( 36 , 37 ). To test this, we simulated PCR with a probabilistic model in which the amplification of an individual template molecule during a PCR cycle was successful based on a binary probability distribution dependent on an amplification efficiency parameter. Iteration of this stochastic model in individual templates over multiple cycles produced a distribution of amplification for individual initial templates as expected (Fig. S11b). However, it but did not generate a distribution with most templates having been lowly amplified, as observed in our experiments. There is also stochasticity introduced by subsampling during sequencing and polymerase amplification error that we did not consider in this simple model ( 38 ). We explored if the stochastic bias affects the relative abundance profiling of communities by running this model on in silico mock communities (Fig. S10c). This stochastic bias is low with higher initial template amounts, but the effect increases as initial template decreases. While the simulations are not definitive, they suggest that PCR can introduce stochastic bias, especially when amplifying a small number of discrete templates. Discussion The advancement of microbiome science requires the development of widely applicable, ultra-high-throughput tools to generate and analyze co-cultures directly derived from the microbiome. By bridging the throughput of droplet cultivation with massively multiplexed sequencing, Cocoa-seq represents an initial effort in translating the coveted reductionist approach of synthetic ecology to complex natural microbiomes. In this work, we multiplexed over a thousand individual droplets co-culture per sequencing library, but the microfluidics throughput is scalable and can be applied to tens of thousands of co-cultures per condition as needed. While benchmarking Cocoa-seq is difficult due to the stochasticity of droplet encapsulation, we benchmarked with multiple metrics and various communities to demonstrate that Cocoa-seq reliably reproduces community distributions and is sensitive enough to detect rare members to be used in answering biological questions. However, Cocoa-seq has limitations, and we provide specific recommendations on its usage. The major limitation of Cocoa-seq is its inability to compare absolute biomass between droplets. To address this limitation, we attempted to implement spike-in standards to provide quantitative metrics per droplet. However, the inherent bias introduced by PCR for single discrete template amplification prevented its use in normalizing read abundance for absolute quantification, which remains a significant technical challenge for Cocoa-seq. Nevertheless, as it is currently presented, Cocoa-seq can be used to explore biological questions with careful and comprehensive experimental design. For example, cell or nucleic acid staining of co-culture beads and subsequent flow-activated cell sorting (FACS) ( 39 ) could enable the selection of co-cultures with high biomass, which then can be multiplexed and sequenced with Cocoa-seq to determine which consortia exhibit facilitative mutualistic interactions in specific contexts. Alternatively, a fluorescent reporter cell which senses conditions of interest ( 16 , 40 ) can be incorporated in droplet co-cultures and likewise used to FACS-screen and study consortia of interest. Additionally, relative abundance alone can be an informative metric with proper experimental design if a compositional shift can be distinguished as a proxy for a specific interaction such as inhibition. However, one specific application that we recommend caution is the analysis of low-cell droplet communities. The droplet communities tested in our benchmarking were relatively high biomass, ranging from a hundred to more than a thousand cells per droplet, which roughly corresponds to at least 10 8 cells/mL. Our investigation of stochastic PCR bias suggests there may be increased risk as total template number decreases and if discrete molecule level resolution is important, such as for studying specific amplicon sequence variants ( 41 , 42 ). Further improvements to Cocoa-seq, such as eliminating the use of agarose, would enable its application to more microbial systems. The current workflow assumes agarose is inert and does not affect the metabolic state of the community. This is not universally true, especially in marine or aquatic systems where agarose is a natural carbon source. In addition, with the current workflow, cultivation temperatures need to be above 37°C to prevent the low-melting temperature agarose from setting and constraining cells during cultivation in droplets, which matters if cells require physical contact for interaction ( 43 ). Environmental microbes also exhibit a very wide array of optimal and relevant growth temperatures and current inability to accommodate this diversity is another limitation. Ultra-low melting point agarose is a potentially feasible alternative since it can remain liquid at room temperature. Alternatively, polyacrylamide could replace agarose as the hydrogel of choice, but its monomer toxicity is understudied ( 44 ) and it may need to be added to droplet co-cultures after cultivation by microfluidic picoinjection ( 45 ). More recently, commercially available semi-permeable capsules enable multi-step processing of microbial cultures ( 46 ) and adapting Cocoa-seq to incorporate these may enable more universal workflows. While we designed Cocoa-seq to investigate complex natural communities, it can also aid traditional synthetic ecology. Sequencing-based synthetic ecology can accommodate a higher number of laboratory strains and analysis of larger consortia. While approaches to high-throughput experimental design currently rely on linear, objective-based progressions through the design space ( 17 , 47 ), the stochastic nature of Cocoa-seq facilitates search through a broader nonlinear design space. The results from an approach like Cocoa-seq can then inform more precise, linear experimental design to determine which initial points in the design space to start with, which may increase the speed of optimization and prevent trapping in local optima. Finally, we want to acknowledge that a major challenge for widespread utilization of microfluidic technology like Cocoa-seq is the requirement of related technical knowhow. To help mitigate the issue, in addition to the Supplementary Information, which details more specifics, we have also posted a protocols.io live document with more nuanced experimental specifications and opportunity for inquiry. We also recommend those who are interested in applying Cocoa-seq to leverage commercially available droplet platforms and/or collaborate with groups that utilize microfluidics at their own or neighboring institutions. Competing Interests The authors declare no conflicts of interest. Data Availability Statement Raw reads are available on NCBI Sequence Read Archive (SRA) (PRJNA1208117). Bioinformatic scripts for read processing and analysis are available on https://github.com/jamesyitan/Cocoa-seq-paper . Microscopy images, including fluorescence images of the bi-culture droplet communities and ddPCR droplets used for analyses are also available in relevant directories within the Github repository. Protocols.io link is available at https://www.protocols.io/view/combinatorial-co-culturing-and-amplicon-sequencing-n2bvj3nxblk5/v1 . Acknowledgements The authors would like to acknowledge helpful conversations with Dr. Freeman Lan in the design of the microfluidic and molecular workflow as well as Dr. Yu-yu Cheng and Dr. Ophelia Venturelli for providing strains. Integrated Training in Microbial Sciences (ITiMS) at the University of Michigan (Burroughs Wellcome) and National Science Foundation (Award 2120909) were critical in providing funding. This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. Footnotes Study funding: National Science Foundation (Award 2120909) (awarded to XNL) https://github.com/jamesyitan/Cocoa-seq-paper https://www.protocols.io/view/combinatorial-co-culturing-and-amplicon-sequencing-n2bvj3nxblk5/v1 References 1. ↵ Koropatkin NM , Cameron EA , Martens EC . How glycan metabolism shapes the human gut microbiota . Nat Rev Microbiol . 2012 ; 10 . 2. Qin J , Li R , Raes J , Arumugam M , Burgdorf KS , Manichanh C , et al. A human gut microbial gene catalogue established by metagenomic sequencing . Nature . 2010 ; 464 ( 7285 ): 59 – 65 . OpenUrl CrossRef PubMed Web of Science 3. ↵ Aldunate M , Srbinovski D , Hearps AC , Latham CF , Ramsland PA , Gugasyan R , et al. Antimicrobial and immune modulatory effects of lactic acid and short chain fatty acids produced by vaginal microbiota associated with eubiosis and bacterial vaginosis . Front Physiol . 2015 ; 6 ( JUN ): 1 – 23 . OpenUrl CrossRef PubMed 4. ↵ Cavicchioli R , Ripple WJ , Timmis KN , Azam F , Bakken LR , Baylis M , et al. Scientists’ warning to humanity: microorganisms and climate change . Nat Rev Microbiol . 2019 ; 17 ( 9 ): 569 – 86 . OpenUrl CrossRef PubMed 5. Singh BK , Bardgett RD , Smith P , Reay DS . Microorganisms and climate change: Terrestrial feedbacks and mitigation options . Nat Rev Microbiol . 2010 ; 8 ( 11 ): 779 – 90 . OpenUrl CrossRef PubMed Web of Science 6. ↵ Michalak AM , Anderson EJ , Beletsky D , Boland S , Bosch NS , Bridgeman TB , et al. Record-setting algal bloom in Lake Erie caused by agricultural and meteorological trends consistent with expected future conditions . Proc Natl Acad Sci U S A . 2013 ; 110 ( 16 ): 6448 – 52 . OpenUrl Abstract / FREE Full Text 7. ↵ Cordero OX , Datta MS . Microbial interactions and community assembly at microscales . Curr Opin Microbiol . 2016 Jun ; 31 : 227 – 34 . OpenUrl CrossRef PubMed 8. ↵ Widder S , Allen RJ , Pfeiffer T , Curtis TP , Wiuf C , Sloan WT , et al. Challenges in microbial ecology: building predictive understanding of community function and dynamics . ISME J . 2016 ;(November 2015): 1 – 12 . 9. ↵ Freilich S , Kreimer A , Meilijson I , Gophna U , Sharan R , Ruppin E . The large-scale organization of the bacterial network of ecological co-occurrence interactions . Nucleic Acids Res . 2010 Feb 27; 38 ( 12 ): 3857 – 68 . OpenUrl CrossRef PubMed Web of Science 10. ↵ Friedman J , Alm EJ . Inferring Correlation Networks from Genomic Survey Data . PLoS Comput Biol . 2012 ; 8 ( 9 ): 1002687 . OpenUrl CrossRef 11. ↵ Henry CS , Dejongh M , Best AA , Frybarger PM , Linsay B , Stevens RL . High-throughput generation, optimization and analysis of genome-scale metabolic models . Nat Biotechnol . 2010 ; 28 ( 9 ): 977 – 82 . OpenUrl CrossRef PubMed Web of Science 12. ↵ Stolyar S , Van Dien S , Hillesland KL , Pinel N , Lie TJ , Leigh JA , et al. Metabolic modeling of a mutualistic microbial community . Mol Syst Biol . 2007 ; 3 : 1 – 14 . OpenUrl CrossRef 13. ↵ Goyal A , Wang T , Dubinkina V , Maslov S . Ecology-guided prediction of cross-feeding interactions in the human gut microbiome . Nat Commun . 2021 Feb 26; 12 ( 1 ): 1 – 10 . OpenUrl CrossRef PubMed 14. ↵ Vrancken G , Gregory AC , Huys GRB , Faust K , Raes J . Synthetic ecology of the human gut microbiota . Nat Rev Microbiol . 2019 ; 17 ( 12 ): 754 – 63 . OpenUrl CrossRef PubMed 15. ↵ Sanchez-Gorostiaga A , Bajić D , Osborne ML , Poyatos JF , Sanchez A . High-order interactions distort the functional landscape of microbial consortia . Vol. 17 , PLoS Biology . 2019 . 1 – 34 p. OpenUrl CrossRef 16. ↵ Kehe J , Kulesa A , Ortiz A , Ackerman CM , Thakku SG , Sellers D , et al. Massively parallel screening of synthetic microbial communities . Proc Natl Acad Sci U S A . 2019 ; 116 ( 26 ): 12804 – 9 . OpenUrl Abstract / FREE Full Text 17. ↵ Clark RL , Connors BM , Stevenson DM , Hromada SE , Hamilton JJ , Amador-Noguez D , et al. Design of synthetic human gut microbiome assembly and butyrate production . Nat Commun . 2021 May 31; 12 ( 1 ): 1 – 16 . OpenUrl CrossRef PubMed 18. ↵ Skwara A , Gowda K , Yousef M , Diaz-Colunga J , Raman AS , Sanchez A , et al. Statistically learning the functional landscape of microbial communities . Nat Ecol Evol . 2023 Oct 2; 7 ( 11 ): 1823 – 33 . OpenUrl CrossRef PubMed 19. ↵ Park J , Kerner A , Burns MA , Lin XN . Microdroplet-Enabled Highly Parallel Co-Cultivation of Microbial Communities . Herrera-Estrella A , editor. PLoS One . 2011 Feb 25; 6 ( 2 ): e17019 . OpenUrl CrossRef PubMed 20. ↵ Tan JY , Wang S , Dick GJ , Young VB , Sherman DH , Burns MA , et al. Co-cultivation of microbial sub-communities in microfluidic droplets facilitates high-resolution genomic dissection of microbial ‘dark matter.’ Integr Biol . 2020 Nov 18; 12 ( 11 ): 263 – 74 . OpenUrl CrossRef 21. ↵ Hsu RH , Clark RL , Tan JW , Ahn JC , Gupta S , Romero PA , et al. Microbial Interaction Network Inference in Microfluidic Droplets . Cell Syst . 2019 ; 9 ( 3 ): 229 – 242 .e4. OpenUrl CrossRef PubMed 22. ↵ Lan F , Demaree B , Ahmed N , Abate AR . Single-cell genome sequencing at ultra-high-throughput with microfluidic droplet barcoding . Nat Biotechnol . 2017 Jul 29; 35 ( 7 ): 640 – 6 . OpenUrl CrossRef PubMed 23. ↵ Zheng W , Zhao S , Yin Y , Zhang H , Needham DM , Evans ED , et al. High-throughput, single-microbe genomics with strain resolution, applied to a human gut microbiome . Science (80-) . 2022 Jun 3; 376 ( 6597 ). 24. ↵ De Rop F V. , Ismail JN , Bravo González-Blas C , Hulselmans GJ , Flerin CC , Janssens J , et al. Hydrop enables droplet-based single-cell ATAC-seq and single-cell RNA-seq using dissolvable hydrogel beads . Elife . 2022 Feb 23; 11 : 1 – 26 . OpenUrl CrossRef PubMed 25. ↵ Zilionis R , Nainys J , Veres A , Savova V , Zemmour D , Klein AM , et al. Single-cell barcoding and sequencing using droplet microfluidics . Nat Protoc . 2017 Dec 8; 12 ( 1 ): 44 – 73 . OpenUrl CrossRef PubMed 26. Klein AM , Mazutis L , Akartuna I , Tallapragada N , Veres A , Li V , et al. Droplet Barcoding for Single-Cell Transcriptomics Applied to Embryonic Stem Cells . Cell . 2015 May ; 161 ( 5 ): 1187 – 201 . OpenUrl CrossRef PubMed 27. ↵ Macosko EZ , Basu A , Satija R , Nemesh J , Shekhar K , Goldman M , et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets . Cell . 2015 ; 161 ( 5 ): 1202 – 14 . OpenUrl CrossRef PubMed 28. ↵ Sheth RU , Li M , Jiang W , Sims PA , Leong KW , Wang HH . Spatial metagenomic characterization of microbial biogeography in the gut . Nat Biotechnol . 2019 ; 37 ( 8 ): 877 – 83 . OpenUrl CrossRef PubMed 29. ↵ Rognes T , Flouri T , Nichols B , Quince C , Mahé F . VSEARCH: A versatile open source tool for metagenomics . PeerJ . 2016 Oct 18; 2016 ( 10 ): e2584 . OpenUrl CrossRef 30. ↵ Schloss PD , Westcott SL , Ryabin T , Hall JR , Hartmann M , Hollister EB , et al. Introducing mothur: Open-source, platform-independent, community-supported software for describing and comparing microbial communities . Appl Environ Microbiol . 2009 ; 75 ( 23 ): 7537 – 41 . OpenUrl Abstract / FREE Full Text 31. ↵ Stoddard SF , Smith BJ , Hein R , Roller BRK , Schmidt TM . rrnDB: improved tools for interpreting rRNA gene abundance in bacteria and archaea and a new foundation for future development . Nucleic Acids Res . 2015 Jan 1; 43 ( Database issue ): D593 . OpenUrl CrossRef PubMed 32. ↵ Mazutis L , Gilbert J , Ung WL , Weitz DA , Griffiths AD , Heyman JA . Single-cell analysis and sorting using droplet-based microfluidics . Nat Protoc . 2013 ; 8 ( 5 ): 870 – 91 . OpenUrl CrossRef PubMed 33. ↵ Tan JY , Saleski TE , Lin XN . The effect of droplet size on syntrophic dynamics in droplet-enabled microbial co-cultivation. Papadimitriou K, editor . PLoS One . 2022 Mar 31 ; 17 ( 3 ): e0266282 . OpenUrl CrossRef PubMed 34. ↵ Stämmler F , Gläsner J , Hiergeist A , Holler E , Weber D , Oefner PJ , et al. Adjusting microbiome profiles for differences in microbial load by spike-in bacteria . Microbiome . 2016 Jun 21; 4 ( 1 ): 1 – 13 . OpenUrl CrossRef PubMed 35. ↵ Kivioja T , Vähärautio A , Karlsson K , Bonke M , Enge M , Linnarsson S , et al. Counting absolute numbers of molecules using unique molecular identifiers . Nat Methods . 2012 Nov 20; 9 ( 1 ): 72 – 4 . OpenUrl CrossRef PubMed Web of Science 36. ↵ Best K , Oakes T , Heather JM , Shawe-Taylor J , Chain B . Computational analysis of stochastic heterogeneity in PCR amplification efficiency revealed by single molecule barcoding . Sci Rep . 2015 Oct 13; 5 ( 1 ): 1 – 13 . OpenUrl CrossRef PubMed 37. ↵ Kou R , Lam H , Duan H , Ye L , Jongkam N , Chen W , et al. Benefits and challenges with applying unique molecular identifiers in next generation sequencing to detect low frequency mutations . PLoS One . 2016 ; 11 ( 1 ). 38. ↵ Blomquist T , Crawford EL , Yeo J , Zhang X , Willey JC . Control for stochastic sampling variation and qualitative sequencing error in next generation sequencing . Biomol Detect Quantif . 2015 Sep 1; 5 : 30 – 7 . OpenUrl CrossRef PubMed 39. ↵ Yanakieva D , Elter A , Bratsch J , Friedrich K , Becker S , Kolmar H . FACS-Based Functional Protein Screening via Microfluidic Co-encapsulation of Yeast Secretor and Mammalian Reporter Cells . Sci Rep . 2020 Jun 23; 10 ( 1 ): 1 – 13 . OpenUrl CrossRef PubMed 40. ↵ Saleski TE , Kerner AR , Chung MT , Jackman CM , Khasbaatar A , Kurabayashi K , et al. Syntrophic co-culture amplification of production phenotype for high-throughput screening of microbial strain libraries . Metab Eng . 2019 ; 54 : 232 – 43 . OpenUrl CrossRef PubMed 41. ↵ Callahan BJ , McMurdie PJ , Rosen MJ , Han AW , Johnson AJA , Holmes SP . DADA2: High-resolution sample inference from Illumina amplicon data . Nat Methods . 2016 ; 13 ( 7 ): 581 – 3 . OpenUrl CrossRef PubMed 42. ↵ Callahan BJ , McMurdie PJ , Holmes SP . Exact sequence variants should replace operational taxonomic units in marker-gene data analysis . ISME J . 2017 Dec 1; 11 ( 12 ): 2639 – 43 . OpenUrl CrossRef PubMed 43. ↵ Pande S , Shitut S , Freund L , Westermann M , Bertels F , Colesie C , et al. Metabolic cross-feeding via intercellular nanotubes among bacteria . Nat Commun . 2015 May 23; 6 ( 1 ): 6238 . OpenUrl CrossRef PubMed 44. ↵ Petka K , Tarko T , Duda-Chodak A . Is Acrylamide as Harmful as We Think? A New Look at the Impact of Acrylamide on the Viability of Beneficial Intestinal Bacteria of the Genus Lactobacillus . Nutrients . 2020 Apr 21; 12 ( 4 ): 1157 . OpenUrl CrossRef PubMed 45. ↵ Eastburn DJ , Sciambi A , Abate AR . Picoinjection Enables Digital Detection of RNA with Droplet RT-PCR . PLoS One . 2013 Apr 26; 8 ( 4 ): e62961 . OpenUrl CrossRef PubMed 46. ↵ Leonaviciene G , Leonavicius K , Meskys R , Mazutis L . Multi-step processing of single cells using semi-permeable capsules . Lab Chip . 2020 ; 20 ( 21 ): 4052 – 62 . OpenUrl CrossRef PubMed 47. ↵ Dama AC , Kim KS , Leyva DM , Lunkes AP , Schmid NS , Jijakli K , et al. BacterAI maps microbial metabolism without prior knowledge . Nat Microbiol . 2023 ; 8 ( 6 ): 1018 – 25 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted January 14, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Massively parallel metabarcoding of droplet co-cultures Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Massively parallel metabarcoding of droplet co-cultures James Y. Tan , Jessica D. Li , Auden Bahr , Xiaoxia Nina Lin bioRxiv 2025.01.14.633022; doi: https://doi.org/10.1101/2025.01.14.633022 Share This Article: Copy Citation Tools Massively parallel metabarcoding of droplet co-cultures James Y. Tan , Jessica D. Li , Auden Bahr , Xiaoxia Nina Lin bioRxiv 2025.01.14.633022; doi: https://doi.org/10.1101/2025.01.14.633022 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Microbiology Subject Areas All Articles Animal Behavior and Cognition (7636) Biochemistry (17703) Bioengineering (13897) Bioinformatics (41960) Biophysics (21458) Cancer Biology (18597) Cell Biology (25524) Clinical Trials (138) Developmental Biology (13381) Ecology (19905) Epidemiology (2067) Evolutionary Biology (24325) Genetics (15612) Genomics (22512) Immunology (17738) Microbiology (40421) Molecular Biology (17187) Neuroscience (88627) Paleontology (667) Pathology (2834) Pharmacology and Toxicology (4825) Physiology (7645) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)