Impact Factor 4.076

The 3rd most cited journal in Microbiology

This article is part of the Research Topic

Microbial food safety along the dairy chain


Front. Microbiol., 22 June 2016 |

Flow Cytometric and 16S Sequencing Methodologies for Monitoring the Physiological Status of the Microbiome in Powdered Infant Formula Production

  • 1Teagasc, Food Research Centre, Fermoy, Ireland
  • 2UCD Centre for Food Safety, Science Centre South, University College Dublin, Dublin, Ireland

The aim of this study was to develop appropriate protocols for flow cytometric (FCM) and 16S rDNA sequencing investigation of the microbiome in a powdered infant formula (PIF) production facility. Twenty swabs were collected from each of the three care zones of a PIF production facility and used for preparing composite samples. For FCM studies, the swabs were washed in 200 mL phosphate buffer saline (PBS). The cells were harvested by three-step centrifugation followed by a single stage filtration. Cells were dispersed in fresh PBS and analyzed with a flow cytometer for membrane integrity, metabolic activity, respiratory activity and Gram characteristics of the microbiome using various fluorophores. The samples were also plated on agar plates to determine the number of culturable cells. For 16S rDNA sequencing studies, the cells were harvested by centrifugation only. Genomic DNA was extracted using a chloroform-based method and used for 16S rDNA sequencing studies. Compared to the dry low and high care zones, the wet medium care zone contained a greater number of viable, culturable, and metabolically active cells. Viable but non-culturable cells were also detected in dry-care zones. In total, 243 genera were detected in the facility of which 42 were found in all three care zones. The greatest diversity in the microbiome was observed in low care. The genera present in low, medium and high care were mostly associated with soil, water, and humans, respectively. The most prevalent genera in low, medium and high care were Pseudomonas, Acinetobacter, and Streptococcus, respectively. The integration of FCM and metagenomic data provided further information on the density of different species in the facility.


Currently, culture-based methods such as agar plates are the most commonly used technique for assessing the microbiological status of food processing environments. In using agar plates, it is only possible to determine the presence of bacteria that the investigator is looking for (using selective agars) and the method gives little information on the physiological status of the cells, except that they are alive if they grow. Some additional information can also be obtained by using different nutrient-content agars and assessing the number of stressed cells. With advances in microbiological techniques, there are a number of different methodologies that can be applied to samples to gain more information about the bacterial population of a processing environment and the physiological state of the bacteria present. Two of these methodologies are 16S rDNA sequencing and flow cytometry (FCM), the former giving information about the type of bacteria present and the latter providing information on the physiological state of the bacteria. Combining the results of these two methods can give valuable information about the microbiological status of a food processing environment.

Though most micro-organisms are cultured using traditional culture-based methods, under stress conditions, many microorganisms, some yet unidentified, are unable to grow on conventional growth media due to lack of effective culture techniques and/or induction of the so-called viable but non-culturable (VBNC) state (see Quigley et al., 2011; Ramamurthy et al., 2014). This has made normal culture techniques ineffective in describing the entire microbiome of complex environments. This problem could be overcome by using 16S rDNA sequencing, which is a culture-independent next generation sequencing method and has been successfully used for describing the composition of the microbiome in depth. The 16S rDNA gene encodes for 16S ribosomal RNA which is universally present in prokaryotic microorganisms (Coenye and Vandamme, 2003). The variations within the 16S rDNA sequences facilitate identification of bacteria at the species level and, therefore, has made 16S rDNA gene sequencing an ideal tool for bacterial taxonomic studies (Neefs et al., 1993). Even though 16S rDNA gene sequencing is widely used to identify inter-species variations, the development of modern high throughput sequencing technologies combined with downstream bioinformatics analysis has made this tool ideal for species identification within complex communities (Logares et al., 2014). It has been shown that 16S rDNA amplicon fragments as short as ~82 base pairs are sufficient for classification at the phylum level (Lazarevic et al., 2009), and with good primer design and analysis methods, fragments of 100–200 base pairs could show the same clustering information as long fragments used in phylogenetic studies (Liu et al., 2007). Compared to other molecular microbiology techniques, 16S rDNA sequencing is not only cheaper in price, but also the interpretation of the resulting data is easier and faster. With some online analysis platforms, such as Illumina Basespace or software such as QIIME and MOTHUR (Schloss et al., 2009; Caporaso et al., 2010) one could easily generate very straightforward information and give an overview about detailed structures and compositions of the target environment microbiome. However, such methodologies will detect gene sequences and will not differentiate live from dead cells.

Flow cytometry is a powerful and rapid technique for simultaneous quantification and multi-parameter analysis of the microbial populations at single cell level (Müller and Nebe-Von-Caron, 2010). Cells are focused and aligned one behind the other in a narrow stream with a diameter close to the diameter of the cells so that single cells can be introduced to the light beam (generally laser). When cells are subjected to the light, they scatter light in all directions, although it is generally detected in two directions: forward scatter (FSC), along the axis of the light source; and side scatter (SSC), perpendicular to the light beam. Data from FSC and SSC scatters are generally used to characterize the morphological state of the cells, as rough indicators of the cell size and granularity, respectively. In addition, the light absorbed by the cells can result in emission of fluorescence (either due to presence of naturally fluorescent compounds or staining with various fluorophores), the intensity of which could be detected by FCM (Shapiro, 2003). Consequently, staining the cells with various fluorophores or fluorescence-conjugated antibodies can be used for understanding a wide range of physiological parameters of the cell (e.g., viability, metabolic and respiratory activities, internal pH, etc.) as well as detection of specific microorganisms at an analysis rate of up to 10,000 cells per second. Furthermore, comparing the viability results obtained with FCM with those of the plate counting can be used to determine the number of VBNC cells in a sample.

The aim of this study was to develop appropriate protocols for flow cytometry and 16S rDNA sequencing investigation of the environmental microbiome in a powdered infant formula (PIF) production facility in the Republic of Ireland.

Materials and Methods


The PIF production unit consisted of three care zones of low, medium and high care. From each zone, twenty swab samples were collected, representing the critical production, storage and packaging sites within each PIF production zone (see Supplementary Table 1). Sponges pre-moistened with neutralizing buffer (Labplas Inc., Sainte-Julie, Canada) were used for environmental swabbing. Each hydrated sponge was used for swiping a single sampling zone of 50 × 50 cm. The zigzag wiping procedure for surface sampling was performed as described by Nicolau and Bolocan (2014). In total, the twenty samples taken represented 5 m2 of that zone and were placed into five bags, each consisting of four sponges. To each bag, 40 mL of phosphate buffered saline (PBS; Sigma, St Louis, USA) was added not earlier than 30 min post-sampling. The purpose of the delay was to prevent the dilution of the neutralizing buffer, allowing the effective neutralization of the possible chlorine and quaternary ammonium compound residues in the sample. The bags were sealed with the tabs provided, kept on ice, transported to the laboratory and processed within 24 h. The sampling procedure and the overall protocol used in this study is schematically represented in Figure 1.


Figure 1. Schematic representation of the steps involved in sampling and sample preparation for flow cytometric and 16S rDNA sequencing studies.

Sample Preparation

Sample Preparation for Flow Cytometric and Plating Studies

The bags were massaged vigorously by hand for 1 min and subsequently stroked in a lab blender (BagMixer® 400 P, Interscience, Saint Nom la Bretèche, France) at a fixed speed of 8 strokes/s for 30 s to release the cells into the PBS dispersion medium. The sponges were then forcefully squeezed by hand and the resultant suspension was transferred into a sterile beaker, representing three composite samples, one from each zone. Clarification of the suspension and separation of cells from debris were achieved by three-step centrifugation followed by a single-step filtration. The cell suspension from each composite sample was first sub-sampled into six polypropylene conical skirted centrifuge tubes (Sarstedt, Wexford, Ireland) and centrifuged at 300 × g for 5 min at 4°C using a Sorvall Legend RT centrifuge (Thermo Electron Corporation, Waltham, USA) to remove the large dust particles and debris. The supernatant was subsequently transferred to a clean centrifuge tube and centrifuged at 4000 × g for 20 min in order to harvest the cells. Finally, the supernatant was transferred again to a clean tube and centrifuged at 4000 × g for 45 min to harvest the possible remaining cells. For the first two steps of centrifugation, the acceleration and breaking were set at 5 and 7, respectively, while for the last step they were set at 9 and 1, respectively. The pellets from steps one, two and three of centrifugation were dispersed in 0.5, 5, and 0.5 mL of 4°C PBS, respectively. The dispersed pellets of the second, third and first stages of centrifugation were then filtered sequentially through sterile Whatman™ Grade 113V prepleated qualitative filter papers (GE Healthcare, Little Chalfont, UK) with pore size of ~30 μm. This was done to ensure that the final sample was free from large particles which could block the fluidics system of the flow cytometer. The filtrate (final sample volume of ~32 mL) was then transferred into a clean centrifuge tube and kept at 4°C before analysis within 24 h.

Sample Preparation for 16S Sequencing Study

For 16S sequencing studies, the sample preparation was similar to that mentioned above with minor modifications. For each centrifugation, 50 ml of the suspension was used. Following the first centrifugation step (300 × g for 5 min at 4°C) and the removal of large dust particles and debris, the supernatant was transferred to a new 50 mL centrifuge tube and centrifuged at 10,000 × g for 10 min at 4°C to harvest the cells. For both steps of the centrifugation, the acceleration and breaking were set at 9. The supernatant was discarded and pellet was resuspended in 1 ml sterile PBS and transferred to a new 1.5 ml microcentrifuge tube. The two centrifuge steps were repeated until all the suspension in the sampling bag was centrifuged. The cell suspensions were stored at −80°C until analyzed.

Flow Cytometric Study


A BD FACSCanto II flow cytometer (BD Bioscience) equipped with green (488 nm air-cooled solid state; 20 mW laser output) and red (633 nm HeNe; 17 mW power output) lasers was used in this study. The fluorescence detector/filters relevant to this study were FITC (FL1; 530 ± 30 nm), PerCP-Cy5.5 (FL3; > 670 nm), APC (FL5; 660 ± 20 nm), and APC-Cy7 (FL6; 780 ± 60 nm). FSC and SSC detectors were also used for determining the light scatter parameters of forward scatter and side scatter, respectively. The instrument was cleaned before and after use and its performance was validated according to the manufacturer's instructions. In order to differentiate between the background noise and the signal (particles of interest), the threshold channel number was set at 200 on FSC. The number of background events (noise) detected by the instrument upon analysis of deionized water (deH2O) at a fast flow rate (≈ 80 μL/min) was less than 4 events/min (approximately less than 1 noise particle per 20 μL deH2O) (data not shown). For each parameter, the height, width and area (integral) of the voltage pulses of each event was measured and recorded.

Determination of the Flow Rate

The exact flow rate of the instrument was determined on the day of the experiment. Ten flow tubes (Sarstedt) were filled with 1 mL of deH2O. Each sample was acquired on the flow cytometer for 30 to 600 s at low, medium or high flow rate settings. The weights of the tubes were determined before and after analysis with the help of an analytical balance with accuracy of ± 0.001 g (Denver instruments, Göttingen, Germany). Assuming the density of deH2O to be 1000 kg/m3, a calibration curve was generated by plotting time vs. volume, to calculate the flow rate (μL/s) of the instrument. Therefore, by knowing the time required for recording a certain number of cells within a sample, it was possible to determine the volume, hence the flow rate of the instrument using Equation 1 (see Supplementary Figure 1).

Flow rate (μLmin)= Number of beads counted ×Sample volume (1,000 μL)Acquisition time (1 min)×number of beads added (5,110 beads)    (1)

Identification of the Cells and Gating Strategy

Depending on the parameter and/or the type of fluorophore used in the study, the Photomultiplier tube (PMT) voltage for each parameter was adjusted so that the cells could be displayed in the center of the investigating plot. The samples were diluted to ensure that the flow rate was between 800 and 1200 events/s at medium speed setting (~45 μl/min). After displaying the data in a density plot of FSC vs. SSC and identifying the particles of interest, the latter was defined by creating a gating region (P1) around it, based on the light scatter properties of the particles (Plot a[1] and b[1]; Figure 2).


Figure 2. Gating strategy used in this study. Cells were acquired (A) before and (B) after staining with SYTO 62 dye. Based on Boolean logic, the events recorded within P1 (20,000 events) were passed through a series of gates (P2-P5) as shown in plots a(2)-a(5) (for unstained cells) and b(2)-b(5) (in the case of stained cells) to determine the number of noise particles (particles in gate P6 of plot a[6]; i.e., 79 noise particles) as well as cells of interest (particles in gate P6 of plot b[6] minus those shown in gate P6 of plot a[6], i.e., 17212–79 = 17133 cells).

The defined gated population of P1 contained not only the presumed cells but also acellular particulates. Therefore, in order to differentiate between the two and detect the cells of interest, samples were stained with ~162 nM SYTO62 (Molecular Probes, Eugene, USA). This dye is a cell-permeant nucleic acid stain which is capable of staining most live and permeabilized (i.e., dead) bacteria. Its maximum excitation and emission wavelengths (λmax) are 652 and 676 nm, respectively. Considering the fact that SYTO62 was the only fluorophore used in this study that could primarily be excited by the red-laser, this made it possible to make exclusive use of the red-laser and its two detectors, FL5 and FL6 for detection of SYTO62-positive particles, hence the cells of interest. The principle behind the multi-gating strategy used in this study and using both FL5 and FL6 for detecting SYTO62 positive particles was based on the one reported by Buzatu et al. (2014). As the λmax emission of SYTO62 is 676 nm, it is primarily detected by FL5; however, due to its broad emission spectrum (620–800 nm), its fluorescence could also be detected by FL6 detector (with relative fluorescence intensity of less than 24% compared to the λmax emission wavelength). The gating strategy used in this study is described in details in Supplementary Data 2.

Determination of Cell Density and Signal to Noise (S/N) Ratio

The cell density in the sample was calculated using the following equation:

Cell denisty (Cellscm2)=[(P6(Stained)-P6(Unstained)) (Cells)]×Sample volume (mL)×1000 (μLmL)×60 (smin)Flow rate (μLmin)*×Acquisition time (s)×Area sampled (cm2)×dilution factor    (2)

P6(stained) and P6(unstained) refer to the number of events within gate P6 for stained and unstained samples, respectively. By knowing the number of events for P6(unstained) (i.e., noise particles), it was also possible to determine the S/N ratio by using the following equation:

Signal to noise ratio (SN)=Number of signal events (cells) within P6Number of noise events within P6    (3)

Physiological Studies and Staining Strategy

Considering the heterogeneity of the microbiome in environmental samples and the variations between the stainability of different microorganisms with various dyes, it was decided to take a holistic multi-staining approach for studying the physiological status of the cells by staining the samples with a wide range of fluorescent dyes. Immediately prior to staining, 250 μL of diluted sample was transferred to 12 × 25 mm flow tubes and supplemented with 20 μL of filter-sterilized 100 mM EDTA (Sigma, Wicklow, Ireland) and 20 μL of 0.1% (v/v in deH2O) polyoxyethylene sorbitan monolaurate (Tween® 20) (Sigma) to improve the stainability of the cells.

The staining protocol for all the fluorochromes used in this study is shown in Table 1. In summary, samples were first stained with 10 μL of 5 μM working solution of SYTO62 to differentiate between the cells of interest and acellular particulates, as previously described. Samples were then vortexed for 1–2 s and incubated at room temperature (18 to 22°C) in darkness for 30 min before analysis. Upon identification of the cells of interest, SYTO62-stained sub-samples were also stained with the following combination of dyes to determine the physiological status of the cells: (a) Propidium iodide (PI) and SYTOX Green Dead Stain to study the membrane integrity; (b) PI in combination with SYTO9 or thiazole orange (TO) to determine the membrane integrity and viability; (c) PI with DiBAC4(3) (Bi-oxonol or BOX) to investigate the membrane potential and viability; (d) Fluorescein diacetate (FDA) and its derivatives [5-(and-6)-carboxyfluorescein diacetate (cFDA) and 5(6)-carboxyfluorescein diacetate N-succinimidyl ester (cFDA-SE)] as indicators of esterase activity; (e) 5-Cyano-2,3-di-(p-tolyl)tetrazolium chloride (CTC) for semi-quantitative analysis of the respiratory activity; and (f) Hexidium iodide (HI) in combination with SYTO9 for Gram staining. All samples were protected from light during the staining process. The protocols used for preparing the fluorophore stock and working solutions as well as the rationale behind using each one is described in detail in Supplementary Data 1.


Table 1. The staining protocol used in this study.

Color Compensation

Live and heat-killed samples of one Gram positive (Lactobacillus rhamnosus GG; LGG) and one Gram negative (Escherichia coli; EC) strain were used as compensation controls. The color compensation was performed for each fluorophore and its primary detector by first plotting a contour plot of the primary detector (the one used for measuring the fluorescence intensity of the fluorophore) vs. non-primary detector. The values for each fluorescence parameter on the plots were transformed bi-exponentially using the BD FACSDiva software version 6.1.3 (BD Biosciences). After gating the negative (non-fluorescent) and positive (fluorescent) cells, the median fluorescence intensity (FI) of the observed populations in primary and non-primary detectors was measured. When the FI of the positive cells in the non-primary detectors was larger than the one for negative cells, a percentage of the FI of the primary detector was subtracted from the affected non-primary FI in order to remove the spillover. The approximate values required for color compensation were calculated using the following equation:

Spillover correction (%)=(median FI of positive control-median FI of positive control)non-primary detector(median FI of positive control-median FI of positive control)primary detector    (4)

The color compensation was performed using the calculated values and verified by visualizing the effects of the applied values on the median FI of the non-primary detectors. If necessary, using the calculated value as a guide, an arbitrary percentage of the FI of the primary detector was subtracted from the affected non-primary FI until the median FI of the latter was relatively the same for both the positive and negative controls. The information on the control samples and the color compensation values used in this study could are shown in Supplementary Tables 2, 3, respectively.

Analysis and Display of FCM Data

The data were acquired and analyzed using BD FACSDiva software version 6.1.3. For FSC vs. SSC contour plots, the scale of the axis for each parameters was logarithmic. When samples were dual-stained (including SYTO62), histograms were used for displaying the FL1 (in the case of FDA, cFDA, cFDA-SE) or FL3 (for CTC). For histograms, the scales for the x-axis (fluorescence channel number) were transformed bi-exponentially (logical x-axis) while the y-axis (count) was displayed on a linear scale beginning at zero. In the case of triple-staining (including SYTO62), the channel numbers (x-axes) for all the fluorescence parameters were transformed bi-exponentially and plotted in two-dimensional quantile (probability) contour plots. All contour plots displayed the events with 99% probability as well as the outliers. The gates (rectangles, polygons, quadrants and vertical or horizontal markers) for each sub-population were set manually.

Plating Study

Growth Media

The growth media used in this study consisted of M9 minimal salt media agar (M9), nutrient agar (NA) (Oxoid, Altrincham, UK) and Brain Heart Infusion (BHI) agar (Sigma). M9, NA, and BHI were used as low, medium and rich nutrient media, respectively, and were used to investigate the effects of the nutrient content of the growth medium on recovery of cells from the processing environment samples. Unless stated otherwise, all media and diluents were prepared according to the manufacturer's instructions using distilled water and autoclaved at 121°C for 15 min. After autoclaving, the salt and agar solutions were allowed to cool to 50°C. The salt broth was then supplemented with 20 mL of 20% (w/v) D-(+)-glucose (Fisher Scientific, Loughborough, UK), 2 mL of 1 M MgSO4 (Fisher Scientific), 0.2 mL of 0.5 M CaCl2 (Reagecon, Shannon, Ireland), and 0.1 mL of 0.5% (w/v) thiamine hydrochloride (Sigma-Aldrich). All the supplements for M9 agar were filter sterilized using 0.22 μm syringe filters (Sarstedt). M9 minimal agar was finally prepared by mixing 500 mL of the supplemented M9 salt broth with 500 mL of 3% (w/v) sterile agar solution.

Recovery of Viable but Non-Culturable Cells (VBNC)

Prior to pouring as agar plates, M9, NA and BHI were also occasionally supplemented with the reactive oxygen species (ROS) scavengers, catalase (Sigma), and sodium pyruvate (SP) (Sigma-Aldrich). To investigate the possible resuscitation of injured and stressed VBNC cells, the supplementation of the agar plates (≈ 20 mL) with 2000 units of catalase per plate (100 μL stock solution) or 0.3% of SP (200 μL stock solution) was performed around 10 min prior to plating the sample. Respectively, 20,000 units/mL and 300 mg/mL stock solutions of catalase and SP were prepared by vigorous vortexing of the compounds in distilled water followed by filter sterilization with 0.22 μm filters. Stock solutions were prepared on the day of experiment and stored at 4°C until use.

Total Aerobic Count (TAC)

In order to determine the total aerobic count (TAC), samples were first decimally serially diluted in maximum recovery diluent (MRD; Fluka), allowed to stand for 30 min after which 100 μL of each dilution was plated on the aforementioned growth media. The plates were incubated either at the room temperature (≈ 21°C), 30°C or 37°C for 48 h. The colonies were examined and counted after 24 h and 48 h incubation.

16S rDNA Sequencing Study

Genomic DNA Extraction and Quality Check

In this study, the processing environment samples were collected from a PIF production site which had strict sanitary standards; therefore, they contained far less cells compared with environment samples collected from the natural environment such as water and soil. As a result, commercial kits with filter column genomic DNA extraction was not used, to avoid filter column clogging and genomic DNA loss during the binding-washing step. Genomic DNA was extracted using a chloroform-based method. Samples (see Section Sample Preparation for 16S Sequencing Study) were centrifuged at 8000 × g, at room temperature for 2 min, the supernatant was discarded and the pellet was resuspended in 1 ml sterile PBS. This washing step was repeated three times after which the pellet was resuspend in 300 μl DNase/RNase free H2O. The resuspended pellet was boiled at 100°C on a heating block for 5 min. After boiling, it was vortexed for ~5 s to ensure the disruption of cell walls and then centrifuged at 8000 × g, for 5 min at room temperature. Subsequently, the supernatant was transferred to a new microcentrifuge tube, mixed with chloroform (Sigma) in a 1:1 ratio and vortexed for 5 s to ensure thorough mixing of the aqueous and chloroform phases. Finally, the mixture was centrifuged at 13000 × g, 4°C for 10 min and ~75% of the aqueous (upper) phase that contained genomic DNA was transferred to a new tube. Nanodrop® Spectrophotomer ND-1000 (1 μl genomic DNA) and Qubit® 2.0 Fluorometer (1 μl genomic DNA; Thermo Fisher Scientific) were used to check genomic DNA concentration and quality, and the genomic DNA was stored at −20°C for further use.

Total RNA Extraction

As with DNA extraction, no commercial kit was used to avoid filter column clogging and total RNA loss during the binding-washing step. Samples were centrifuged at 8000 × g, at room temperature for 2 min, the supernatant was discarded and the pellet was resuspended and washed with 1 ml sterile PBS. The centrifuge and wash steps were repeated twice and the pellet was resuspend in 1 ml sterile PBS, transferred to a 50 ml centrifuge tube and sterile PBS was used to adjust the final volume of the cell suspension to 4 ml. To this suspension, 1.6 mL of ice cold phenol-ethanol solution (95% ethanol and 5% acidic phenol; pH 4.3) was added and the tube was incubated on ice for at 30–120 min to stabilize the RNA and prevent degradation (Tedin and Bläsi, 1996). The mixture was centrifuged at 3300 × g at 4°C for 10 min and most of the supernatant was discarded. The pellet was resuspended with the remaining supernatant in the tube and transferred to a 1.5 ml tube. The tube was centrifuged at 18,000 × g, 4°C for 1 min, the supernatant was discarded and the pellet was resuspended in 1 ml of TRIzol (Ambion, Foster City, USA). The mixture was transferred to a 2 ml heavy phase lock tube (5 Prime) and supplemented with 400 μl chloroform (Sigma-Aldrich). The tube was immediately gently inverted for 10 s (no vortexing) and incubated at room temperature for 2–5 min. The mixture was centrifuged at 16,000 × g, at room temperature for 15 min and the aqueous phase was transferred to a new microcentrifuge tube, supplemented with 450 μl isopropanol (Sigma-Aldrich) and mixed immediately. The mixture was incubated at room temperature for 20 min and stored at −20°C overnight to provide higher yields. After storage, the mixture was centrifuged at 16,000 × g at room temperature for 30 min, the supernatant was discarded, the pellet was washed with 350 μl of 70–75% ethanol and centrifuged at 16,000 × g at room temperature for 10 min. After air drying the pellet, 25 μl of pre-heated (65°C) DNA/RNase-free H2O was added and the tube was incubated on a thermomixer at 900 rpm, 65°C for 5 min. During this time, the tubes were vortexed briefly 2–3 times to improve the dispersion of the pellet. The liquid in the tube contained total RNA. DNA removal was carried out using Ambion TURBO DNA-free™ kit (Thermo Fisher Scientific, USA). An Agilent 2100 Bioanalyzer (1 μl total RNA) and a Nanodrop® Spectrophotomer ND-1000 (1 μl total RNA) were used to check total RNA quantity and quality, and the sample was stored at −80°C for further use.

16S rDNA Sequencing and Bioinformatics Analysis

16S rDNA sequencing was carried out using the Illumina MiSeq platform. Such sequencing is amplicon-based, targeting the variable region V3-V4 of the 16S rDNA gene and generates PCR products with a length of ~460 bp. Paired-end, 300 × 2 bp sequencing was carried out to cover the whole PCR product from the two opposite ends. For library construction, the Illumina official guide for 16S rDNA sequencing library preparation was used as reference. Bioinformatics analysis was carried out using the 16S Metagenomics app v1.0.1 provided in the Illumina BaseSpace online platform, including reads quality control, reads alignment, assembly and annotations. For taxonomic classification of the 16S rDNA reads, the Ribosomal Database Project (RDP) Classifier was used for classification (Wang et al., 2007) and an Illumina-curated version of the GreenGenes taxonomic database was used as reference database (DeSantis et al., 2006). Both applications were inbuilt in BaseSpace.

Results and Discussion

Flow Cytometric and Plating Studies

Table 2 shows the FCM data obtained for three samples collected from low, medium and high care zones. The greatest number of cells (regardless of their physiological state) was detected in the low care zone, followed by the medium and high care zones (p < 0.0001). This was expected considering the implementation of stricter levels of hygiene and working practices by the PIF manufacturer in the latter two zones. However, the reduction in the total cell count, did not necessarily translate to a concurrent reduction in the number of viable cells as determined by FCM. By knowing the total cell count and the percentage of viable cells (based on the exclusion of PI viability dye), it was possible to calculate the density of viable cells per cm2 of the sampling zone. The number of viable cells per cm2 in medium care was nearly 3 times greater than that detected in both low and high care zones. This was probably due to greater level of humidity, hence greater access of microorganisms to available water in this zone. The lack of significance between the FCM viability results obtained for the two dry zones (low and high care) could be considered as further evidence of the primary role of humidity in improved viability of the cells in the wet medium care zone.


Table 2. Flow cytometric (FCM) total cell and viable cell density.

The results of the flow cytometry study were in agreement with those of the plate counting technique in the sense that the cells from medium care exhibited significantly greater culturability as well as viability compared to those from the other two care zones (Figure 3). In both dry zones, compared to the FCM, the plate counting significantly under-estimated the number of viable cells, while in the medium care, significantly lower number of cells were detected by FCM, compared to the plate counting. The under-estimation of viable cells by plate counting due to the presence of stressed and starved VBNC cells is well-documented (Oliver, 2005) and was, therefore, expected. On the other hand, the under-estimation of viable cells by FCM in medium care may be due to the presence of ultramicrobacteria and ultramicrocells (e.g., bacteria of the genera Flavobacterium, Bacteroides, and Chryseobacterium) with sizes smaller than the detection limit of the flow cytometer used in this study (< 0.5 μm). This could have resulted in the cells being considered as background noise. This could also mean that the number of FCM viable cells, hence VBNC cells in dry care zones of low and high could have been significantly higher than determined in this study. Moreover, the possible presence of very high concentration of surfactants, detergents and washing solutions in that zone, could have rendered some of the cells un-stainable with SYTO62, hence not detectable based on the proposed protocol (Vives-Rego et al., 1999).


Figure 3. Total viable count (log10 cells/cm2) based on the flow cytometry (FCM) and plate counting (CFU) techniques. For FCM, a cell was considered as viable if it could not be stained with either of PI or SYTOX Green Dead Stain, while for CFU, it formed a colony on nutrient agar plate at 37°C, 48 h post-inoculation. The FCM values are the same as those calculated in row N of Table 2. All the values are the mean ± SD of two technical replicates. Capital letters are used for comparing the FCM or CFU data between each zone, while lower-case letters compare the FCM and CFU values within each zone. Columns marked with similar capital or lower-case letters are not statistically significantly different based on unpaired two-tailed Student's t-test; p > 0.05) (Data from sampling in October 2015).

Table 3 shows the FCM data regarding the physiological status as well as Gram characteristics of the cells in each care zone. With regard to membrane integrity, membrane polarization, and metabolic activity, relatively similar results were obtained. The greatest mean percentage of cells with intact and polarized membranes and metabolic activity was observed in medium care, followed by high and low care zones. For instance, with the exception of SYTO9/PI staining, the percentage of cells with intact membranes was significantly lower in low care compared to the other two zones, regardless of the fluorophores used. The discrepancies observed between the results obtained for each fluorophore combination could be due to the difference between their staining mechanism (Netuschil et al., 2014). Although, the greatest percentage of cells with respiratory activity was found in the medium care sample, the values obtained were significantly higher than the percentage of viable cells. The reason behind this observation is not clear, however, it could be due to the residual activity of the electron transfer chain in cells with depolarized and compromised cells. Moreover, the difference between the stainability of Gram negative and Gram positive cells could have played a role in the discrepancies observed (Holm and Jespersen, 2003).


Table 3. The physiological status and Gram characteristics of the microbiome.

Changing the incubation temperature and or the supplementation of the growth media had a significant effect on the recovery of the cells on a care-zone basis (Figure 4 and Supplementary Table 4). For instance, decreasing the incubation temperature from 37 to 30°C or room temperature, had a highly significant effect on the total aerobic count for samples collected from the low care zone (p < 0.0001). This could indicate the presence of a significant number of psychrotrophic bacteria in low care samples (Hantsis-Zacharov and Halpern, 2007). The change in nutrient content of the media did not make a statistically significant difference in the recovery of the stressed cells. On the other hand, supplementation of the growth media with reactive oxygen species scavengers seemed to improve the recovery of the cells. Catalase and SP have previously been shown to be effective in recovery and resuscitation of injured and VBNC cells (Mizunoe et al., 2000; Bang et al., 2007). Catalase was effective in recovery of the cells collected from both medium and high care zones, while SP only improved the recovery in medium care samples. In the case of low care samples, the results were similar to those observed in the medium care zone, however, they did not reach the level of statistical significance.


Figure 4. Recovery of stressed cells using catalase and sodium pyruvate. Samples were spread plated in duplicate on nutrient agar (NA), M9 Salt agar or Brain Heart Infusion agar (BHI) solid growth media (with or without catalase or SP) and incubated at either of 21, 30, or 37°C. For each temperature/sample combination (e.g., low care sample at 37°C) shown, the data are the mean ± SD (technical duplicate plating of a single sample) of the plate count obtained from all the plates that were subjected to the indicated parameter, regardless of the effects of others (n = 18). For instance, the plate count for the low care sample at 37°C (3.17 ± 0.17) is the mean ± SD of the plate counts for all the eighteen plates that were incubated at 37°C, regardless of the growth medium or supplementation. See Supplementary Table 4 for further information on the three-way interaction between temperature, media and supplemnetation variables for each sample (Data from sampling in May 2015). (**: p < 0.01; ****: p < 0.0001).

16S rDNA Sequencing Study

The mean DNA content in low, medium and high care zones was 278, 168, and 53 pg/cm2, respectively. Considering the vast difference between the DNA content of different bacterial species, it was not possible to establish a direct correlation between the DNA content of the sample and the cell count. Nonetheless, the greatest DNA content was found in the low care sample, followed by medium and high care zone samples, which closely resembled the results for cell counts obtained using from the FCM total count, as previously described. In addition, the mean total RNA content for samples of low, medium and high care zones was 2, 29, and 1 pg/cm2, respectively. By making the presumption that the presence RNA in the cell is an indicator of protein synthesis, hence a degree of cellular vitality, the RNA content of the cells was in agreement with both the FCM viable count and plate counts.

16S sequencing provided valuable information on the type of bacteria present in each care zone and the percentage contribution of each genera to the entire microbiome. This technique has been used successfully to characterize the microbiome in various environments such as soil, water, food and hospitals (Gomez-Alvarez et al., 2012; Oberauner et al., 2013; Rampelotto et al., 2013; de Boer et al., 2015). By integrating these data with that of the total FCM cell count, it was possible to calculate the number of cells belonging to each genera in each care zone. 16S sequencing of the samples revealed the presence of 243 bacterial genera (with more than 0.05% distribution) in the microbiome of the PIF production facility (Figure 5). The combination of 16S rRNA-targeted oligonucleotide probes with flow cytometry for analyzing mixed microbial populations has been reported (Amann et al., 1990), although there are no reports of its use to characterize the microbial population of a food production environment. The results pointed to a striking similarity between the type of bacteria present in different care zones. For instance, 42 out of 243 genera were common to all three care zone, contributing to nearly 70% of the microbiome. Similarly, 58 genera were common between low care and the other two care zone. On the other hand, although a third of the identified genera were unique to low care, they only made up 7.3% of the microbiome (total 1180 cells/cm2). Similarly, 19 and 39 genera were unique to medium and high care zones, making up 1.0 and 2.4% of the microbiome, respectively.


Figure 5. Schematic representation of the number of genera identified in each care zone. The size of each circle and the overlap areas are propotionate to the number of genera idnetified in that zone. The values outside the parantheses show the number of genera idenfied within the zone/area. The values within the parentheses show the total cells/cm2 and the percentage contribution of those genera to the overall microbiome of the PIF production unit. The total cell counts exclude the unpecified genera and/or species (2682, 540, and 819 cells/cm2 of unspecified bacteria in low, medium and high care zones, respectively). This Venn diagram was generated using the BioVenn software by Hulsen et al. (2008).

Looking at the top thirty genera, in terms of % occurrence, in the PIF production unit provided a better picture of the type of microorganisms associated with each care zone (Table 4). For the complete list, readers are referred to Supplementary Table 5. Twenty out of 30 genera which were predominantly present in low care zone, are mainly associated with soil and the general environment, which included species belonging to Pseudomonas, Spirosoma, and Sphingomonas genera. On the other hand, those predominantly present in wet care zone such as Acinetobacter, Chryseobacterium, and Paucibacter are mainly associated with water and sewage, as well as soil and other general environment sources. In contrast, the greatest number of human and milk-associated genera such as Streptococcus, Lactococcus, Corynebacterium, Lactobacillus, and Kocuria were found in the high care zone. Washing the sponges of the low care zone in PBS resulted in the formation of very dark gray cell suspension with a significant soil and debris sedimentation. Unlike the other two care zones, the majority of drains and sampling points (mainly drains) within medium care were wet. Based on the current data, it was not possible to definitively determine the primary reason behind the greater prevalence of human-associated microorganisms in high care zone. However, it is plausible that, while strict segregation of the high care zone led to a substantial reduction in the entry of soil and drain-associated microorganisms from the other two care zones, human-associated microorganisms within that zone still contributed to the microbiome of high care zone.


Table 4. The top 30 genera in the PIF production unit (all care zones).

A closer look at the number of cells of the different species of the top three genera could be used as a good indicator of the possible transition of the cells between different zones (Figure 6). For instance, with regard to the Acinetobacter, the top five species of this genus in low care and the top three in low care were also among the top five species of this genus in medium care. Similar results were also obtained for Streptococcus spp. where S. vestibularis, S. bovis, and S. fryi were the top three species of this genus in all three zones. Furthermore, the rate of change in the number of cells for each species in one zone, closely resembled the change observed in the other two zones. It is important to note that the aim of this study was to compare the microbiome of different care zones and therefore, sponges from different locations of a specific care zone were placed in a single bag. Consequently, this could have contributed to the variability observed between the results for each care zone. Further studies are needed to determine the microbiome of each sampling point.


Figure 6. The number of cells for each species of (A) Acinetobacter, (B) Streptococcus and (C) Pseudomonas was determined by multiplying the percentage distribution of each species in each care zone (as determined by 16S rDNA sequencing) by the total cell count for the corresponding care zone (as determined by FCM) (Data from sampling in May 2015).

16S sequencing also provided information on the pathogenic strains present in each care zone. Table 5 shows the percentage distribution of 18 pathogenic species and 1 pathogenic genus in three different care zones. According to the official report of FAO/WHO (2004), these species are divided into three categories “based on the strength of evidence of a causal association between their presence in PIF and illness in infants.” Class A includes Cronobacter spp. and Salmonella enterica for which clear evidence of causality exist. No Class A microorganisms were detected in either of the three care zones. On the other hand, Class B (i.e., causality plausible, but not yet demonstrated) and Class C organisms (i.e., causality less plausible, or not yet demonstrated) were detected at both genus and species level in all three zones.


Table 5. Pathogenic species identified in different care zones.

In summary, the results showed that the physical segregation of a production unit into different care zones has a positive impact on reducing the microbial load within a PIF production unit. However, the reduction in total cell count did not lead to a reduction in either the total viable count or the human associated pathogenic bacterial species. Therefore, better control measures such as stricter monitoring of staff and personal hygiene policies might be necessary to achieve a significant reduction in the human-associated microorganisms in high care. The results also demonstrated that combining the FCM and 16S rDNA sequencing data could be used successfully for hygiene monitoring in PIF production units.

Author Contributions

All authors contributed to the design of the sampling protocol, AA and YC collected samples for FCM and 16S rDNA sequencing studies, respectively, prepared the samples, analyzed the data and drafted the manuscript. AA integrated the FCM and sequencing data, wrote the statistical analysis plan and analyzed the integrated data. SS contributed to the development of sample preparation protocol and analyzed the sequencing data. SF and KJ initiated the study, obtained the funding, designed and supervised the collaborative project, monitored the data collection and contributed to the completion of the manuscript. AA and YC made an equal contribution (co first authors).


This work was funded by the Irish Department of Agriculture, Food and the Marine, Ireland, under the Food Institutional Research Measure (FIRM) initiative (project Smart-PIF; reference number 13/F/423).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


The authors would like to acknowledge the co-operation of the management and staff in the PIF production facility and the provision of access for AA and YC.

Supplementary Material

The Supplementary Material for this article can be found online at:


Amann, R. I., Binder, B. J., Olson, R. J., Chisholm, S. W., Devereux, R., and Stahl, D. A. (1990). Combination of 16S rRNA-targeted oligonucleotide probes with flow cytometry for analyzing mixed microbial populations. Appl. Environ. Microbiol. 56, 1919–1925.

PubMed Abstract | Google Scholar

Bang, W., Drake, M. A., and Jaykus, L. A. (2007). Recovery and detection of Vibrio vulnificus during cold storage. Food Microbiol. 24, 664–670. doi: 10.1016/

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzatu, D. A., Moskal, T. J., Williams, A. J., Cooper, W. M., Mattes, W. B., and Wilkes, J. G. (2014). An integrated flow cytometry-based system for real-time, high sensitivity bacterial detection and identification. PLoS ONE 9:e94254. doi: 10.1371/journal.pone.0094254

PubMed Abstract | CrossRef Full Text | Google Scholar

Caporaso, J. G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F. D., Costello, E. K., et al. (2010). QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336. doi: 10.1038/nmeth.f.303

PubMed Abstract | CrossRef Full Text | Google Scholar

Coenye, T., and Vandamme, P. (2003). Intragenomic heterogeneity between multiple 16S ribosomal RNA operons in sequenced bacterial genomes. FEMS Microbiol. Lett. 228, 45–49. doi: 10.1016/S0378-1097(03)00717-1

PubMed Abstract | CrossRef Full Text | Google Scholar

de Boer, P., Caspers, M., Sanders, J. W., Kemperman, R., Wijman, J., Lommerse, G., et al. (2015). Amplicon sequencing for the quantification of spoilage microbiota in complex foods including bacterial spores. Microbiome 3:1. doi: 10.1186/s40168-015-0096-3

PubMed Abstract | CrossRef Full Text | Google Scholar

DeSantis, T. Z., Hugenholtz, P., Larsen, N., Rojas, M., Brodie, E. L., Keller, K., et al. (2006). Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl. Environ. Microbiol. 72, 5069–5072. doi: 10.1128/A.E.M.03006-05

PubMed Abstract | CrossRef Full Text | Google Scholar

FAO/WHO (2004). Enterobacter Sakazakii and Salmonella in Powdered Infant Formula. [online] Available online at: (Accessed 11 March 2016).

Gomez-Alvarez, V., Revetta, R. P., and Santo Domingo, J. W. (2012). Metagenomic analyses of drinking water receiving different disinfection treatments. Appl. Environ. Microbiol. 78, 6095–6102. doi: 10.1128/AEM.01018-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Hantsis-Zacharov, E., and Halpern, M. (2007). Culturable psychrotrophic bacterial communities in raw milk and their proteolytic and lipolytic traits. Appl. Environ. Microbiol. 73, 7162–7168. doi: 10.1128/AEM.00866-07

PubMed Abstract | CrossRef Full Text | Google Scholar

Holm, C., and Jespersen, L. (2003). A flow-cytometric gram-staining technique for milk-associated bacteria. Appl. Environ. Microbiol. 69, 2857–2863. doi: 10.1128/AEM.69.5.2857-2863.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Hulsen, T., de Vlieg, J., and Alkema, W. (2008). BioVenn — a web application for the comparison and visualization of biological lists using area-proportional Venn diagrams. BMC Genomics 9:488. doi: 10.1186/1471-2164-9-488

PubMed Abstract | CrossRef Full Text | Google Scholar

Lazarevic, V., Whiteson, K., Huse, S., Hernandez, D., Farinelli, L., Østerås, M., et al. (2009). Metagenomic study of the oral microbiota by Illumina high-throughput sequencing. J. Microbiol. Meth. 79, 266–271. doi: 10.1016/j.mimet.2009.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Lozupone, C., Hamady, M., Bushman, F. D., and Knight, R. (2007). Short pyrosequencing reads suffice for accurate microbial community analysis. Nucleic Acids Res. 35:e120. doi: 10.1093/nar/gkm541

PubMed Abstract | CrossRef Full Text | Google Scholar

Logares, R., Sunagawa, S., Salazar, G., Cornejo-Castillo, F. M., Ferrera, I., Sarmento, H., et al. (2014). Metagenomic 16S rDNA Illumina tags are a powerful alternative to amplicon sequencing to explore diversity and structure of microbial communities. Environ. Microbiol. 16, 2659–2671. doi: 10.1111/1462-2920.12250

PubMed Abstract | CrossRef Full Text | Google Scholar

Mizunoe, Y., Wai, S. N., Ishikawa, T., Takade, A., and Yoshida, S. I. (2000). Resuscitation of viable but nonculturable cells of Vibrio parahaemolyticus induced at low temperature under starvation. FEMS Microbiol. Lett. 186, 115–120. doi: 10.1111/j.1574-6968.2000.tb09091.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Müller, S., and Nebe-Von-Caron, G. (2010). Functional single-cell analyses: flow cytometry and cell sorting of microbial populations and communities. FEMS Microbiol. Rev. 34, 554–587. doi: 10.1111/j.1574-6976.2010.00214.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Neefs, J.-M., Van de Peer, Y., De Rijk, P., Chapelle, S., and De Wachter, R. (1993). Compilation of small ribosomal subunit RNA structures. Nucleic Acids Res. 21, 3025–3049. doi: 10.1093/nar/21.13.3025

PubMed Abstract | CrossRef Full Text | Google Scholar

Netuschil, L., Auschill, T. M., Sculean, A., and Arweiler, N. B. (2014). Confusion over live/dead stainings for the detection of vital microorganisms in oral biofilms-which stain is suitable? BMC Oral Health 14:2. doi: 10.1186/1472-6831-14-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicolau, A. I., and Bolocan, A. S. (2014). “Sampling the processing environment for Listeria” in Listeria Monocytogenes Methods and Protocols, eds K. Jordan, E. M. Fox, and M. Wagner (New York, NY: Humana Press), 3–14.

Oberauner, L., Zachow, C., Lackner, S., Högenauer, C., Smolle, K.-H., and Berg, G. (2013). The ignored diversity: complex bacterial communities in intensive care units revealed by 16S pyrosequencing. Sci. Rep. 3:1413. doi: 10.1038/srep01413

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliver, J. D. (2005). The viable but nonculturable state in bacteria. J. Microbiol. 43, 93–100.

PubMed Abstract | Google Scholar

Quigley, L., O'Sullivan, O., Beresford, T. P., Ross, R. P., Fitzgerald, G. F., and Cotter, P. D. (2011). Molecular approaches to analysing the microbial composition of raw milk and raw milk cheese. Int. J. Food Microbiol. 150, 81–94. doi: 10.1016/j.ijfoodmicro.2011.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramamurthy, T., Ghosh, A., Pazhani, G. P., and Shinoda, S. (2014). Current Perspectives on Viable but Non-Culturable (VBNC) Pathogenic Bacteria. Front. Public Health 2:103. doi: 10.3389/fpubh.2014.00103

PubMed Abstract | CrossRef Full Text | Google Scholar

Rampelotto, P. H., de Siqueira Ferreira, A., Barboza, A. D. M., and Roesch, L. F. W. (2013). Changes in diversity, abundance, and structure of soil bacterial communities in Brazilian Savanna under different land use systems. Microb. Ecol. 66, 593–607. doi: 10.1007/s00248-013-0235-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Schloss, P. D., Westcott, S. L., Ryabin, T., Hall, J. R., Hartmann, M., Hollister, E. B., et al. (2009). Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541. doi: 10.1128/AEM.01541-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Shapiro, H. M. (2003). Practical Flow Cytometry, 4th Edn. Hoboken, NJ: John Wiley & Sons.

Google Scholar

Tedin, K., and Bläsi, U. (1996). The RNA chain elongation rate of the lambda late mRNA is unaffected by high levels of ppGpp in the absence of amino acid starvation. J. Biol. Chem. 271, 17675–17686. doi: 10.1074/jbc.271.30.17675

PubMed Abstract | CrossRef Full Text | Google Scholar

Vives-Rego, J., Guindulain, T., Vazquez-Dominguez, E., Gasol, J. M., Lopez-Amoros, R., Vaque, D., et al. (1999). Assessment of the effects of nutrients and pollutants on coastal bacterioplankton by flow cytometry and SYTO-13 staining. Microbios 98, 71–85.

Google Scholar

Wang, Q., Garrity, G. M., Tiedje, J. M., and Cole, J. R. (2007). Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261–5267. doi: 10.1128/AEM.00062-07

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: powdered infant formula (PIF), systems microbiology, flow cytometry, 16S sequencing, viable but non-culturable (VBNC), microbial stress response, environmental sampling, microbial physiology

Citation: Anvarian AHP, Cao Y, Srikumar S, Fanning S and Jordan K (2016) Flow Cytometric and 16S Sequencing Methodologies for Monitoring the Physiological Status of the Microbiome in Powdered Infant Formula Production. Front. Microbiol. 7:968. doi: 10.3389/fmicb.2016.00968

Received: 11 April 2016; Accepted: 03 June 2016;
Published: 22 June 2016.

Edited by:

Juan Aguirre, Universidad de Chile, Chile

Reviewed by:

Nuno F. Azevedo, University of Porto, Portugal
Julio Parra-Flores, Universidad del Bío-Bío, Chile

Copyright © 2016 Anvarian, Cao, Srikumar, Fanning and Jordan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Kieran Jordan,

These authors have contributed equally to this work and co-first authors.