Combined analysis of microbial metagenomic and metatranscriptomic sequencing data to assess in situ physiological conditions in the premature infant gut.

Microbes alter their transcriptomic profiles in response to the environment. The physiological conditions experienced by a microbial community can thus be inferred using meta-transcriptomic sequencing by comparing transcription levels of specifically chosen genes. However, this analysis requires accurate reference genomes to identify the specific genes from which RNA reads originate. In addition, such an analysis should avoid biases in transcript counts related to differences in organism abundance. In this study we describe an approach to address these difficulties. Sample-specific meta-genomic assembled genomes (MAGs) were used as reference genomes to accurately identify the origin of RNA reads, and transcript ratios of genes with opposite transcription responses were compared to eliminate biases related to differences in organismal abundance, an approach hereafter named the "diametric ratio" method. We used this approach to probe the environmental conditions experienced by Escherichia spp. in the gut of 4 premature infants, 2 of whom developed necrotizing enterocolitis (NEC), a severe inflammatory intestinal disease. We analyzed twenty fecal samples taken from four premature infants (4-6 time points from each infant), and found significantly higher diametric ratios of genes associated with low oxygen levels in samples of infants later diagnosed with NEC than in samples without NEC. We also show this method can be used for examining other physiological conditions, such as exposure to nitric oxide and osmotic pressure. These study results should be treated with caution, due to the presence of confounding factors that might also distinguish between NEC and control infants. Nevertheless, together with benchmarking analyses, we show here that the diametric ratio approach can be applied for evaluating the physiological conditions experienced by microbes in situ. Results from similar studies can be further applied for designing diagnostic methods to detect NEC in its early developmental stages.

Introduction

Physiological conditions within the gut are important metrics to measure when studying gut inflammatory diseases [1,2], yet are notoriously difficult to measure in vivo. Transcriptional profiling provides information on the pool of genes that microbial cells express, and therefore can reveal the physiological conditions experienced by these cells. Microbial transcriptional patterns have been analyzed using many methods, including reverse transcription quantitative PCR [3], microarrays [4], and in recent years, RNA sequencing [5]. Analyzing transcription patterns within microbial communities, i.e. meta-transcriptomics, is challenging because it is necessary to specifically identify the microbial species from which transcripts originate. In addition, it is also important to account for changes in the abundance of different microbial species from which transcripts originate, as changes in transcript abundance can also be related to changes in organismal abundances [6].

A relevant case in which such an approach would be useful is the study of necrotizing enterocolitis (NEC) in premature infants. The hallmark of NEC is inflammation of the small and/or large bowel that can progress rapidly to intestinal necrosis, sepsis, and death [7,8]. Because the onset of the disease is often fulminant, treatment options for severe cases are limited and often futile. Thus, the need for disease biomarkers, to enable early and accurate diagnosis, motivates ongoing research on early stages of NEC development. A recent meta-analysis of 14 DNA-based studies reported that fecal samples from preterm infants, later diagnosed with NEC, contained modest but statistically significant increased abundance of facultative anaerobes from the Proteobacteria phylum and a modest decrease in abundance of strict anaerobes [9]. As others have previously noted, an understudied approach to effectively study the bacterial response to local conditions within the infant gut is to pair taxonomic profiling with functional information [9,10]. This understudied approach can be addressed by applying the transcriptomic analysis approaches described here.

Here we established a new method of measuring transcriptomic data, named the diametric ratio, to measure physiological conditions from microbial community transcriptomic data. We used this approach to study NEC in a pilot cohort of premature infants, to infer about the transcriptional patterns associated with physiological conditions occurring before NEC is diagnosed, specifically targeting genes related to oxygen exposure.

Methods

Study design and sampling

This study made use of a previously analyzed dataset of metagenomic DNA sequencing of premature infant stool samples [10-13]. In this study a subset of these stool samples from four infants, two of which were diagnosed with NEC, were additionally subjected to RNA sequencing. As previously described [12], stool samples for establishing these datasets were collected after perineal stimulation [14], so that fecal samples were collected under direct vision immediately upon evacuation, to minimize changes in the transcription pattern of gut microbes outside the intestine. After collection of the samples, they were stored at -80°C until DNA and RNA extraction. Medical information for these infants is presented in Table 1, and sampling schedule for each infant is available in Table 2. To analyze changes in transcription patterns over time and with a balanced distribution of samples, five samples from both NEC and Control infants from day of life (DOL) 10–19 were grouped over the 1st time block, and those from DOL 20–37 were grouped over the 2nd time block (Table 2).

Table 1

Infant medical information.

Infant	Gestational age (Weeks)	Study	Gender	Delivery	Weight (g)	Feeding	Condition	NEC Diagnosis (DOL)
64	28	NIH2	M	Vaginal	1100	Combination	Control	–
66	28	NIH2	F	Vaginal	1028	Breast	Control	–
69	26	NIH2	M	C-section	637	Combination	NEC	32
71	25	NIH2	M	C-section	754	Combination	NEC	31

Table 2

Infant stool sampling scheme.

Diagnosis	infant	1st time block								2nd time block
NEC	71	✓				✓			✓		✓	✓	✓		*
	69				✓			✓					✓	✓		*
Control	66							✓	✓						✓		✓		✓
	64		✓	✓			✓			✓								✓
		10	11	13	14	15	16	17	19	20	21	25	27	28	31	32	33	36	37
		Days of life

indicate dates stool samples were taken

indicate NEC diagnosis

DNA and RNA sequencing

Procedures for DNA extraction and sequencing were previously described [13]. RNA was extracted from selected stool samples using MOBIO PowerMicrobiome RNA isolation kit. The only modification to the manufacturer’s protocol was that phenol:chloroform:isoamyl alcohol was added to the glass bead tubes prior to the addition of the stool sample. RNAseq libraries were prepared with Illumina's 'TruSeq Stranded RNAseq Sample Prep kit'. Prior to library preparation, eukaryotic rRNA was removed using the 'Ribo-Zero rRNA Removal Kit (Human/Mouse/Rat)', and bacterial rRNA was removed using the 'Ribo-Zero rRNA Removal Kit (Bacteria)' from Illumina. DNA and RNA were sequenced on HiSeq2500, RNA sequencing was performed at the Roy J. Carver Biotechnology Center at the University of Illinois at Urbana-Champaign. DNA sequencing yielded an average of 35,127,007 reads per sample, while RNA sequencing of corresponding samples yielded an average of 13,786,716 reads per sample.

Genome reconstruction

DNA reads were trimmed using Sickle (v1.33) (https://github.com/najoshi/sickle), and assembled into scaffolds with IDBA_UD (v1.1.1) [15]. Scaffolds were binned into genomes using DAS tool (v1.0), which uses a combination of established binning algorithms [16]. Reconstructed genomes from each infant were de-replicated according to 99% average nucleotide identity using dRep (v2.3.2) [17]. One genome of each identified bacterial species was manually chosen from each infant gut microbiome for downstream RNA analysis. Scaffolds from MAG’s are accessible on ggkbase interface (https://ggkbase.berkeley.edu/), with a ‘scaffold to bin’ file available as supporting information files that can be used to reconstruct all selected genome bins used in this study (S3–S6 Files). In addition, Table A in S1 File shows accessions for the different Escherichia spp. genes used in this study.

Phylogenetic profiling

Taxonomic classification was done according to ‘shared affiliation of predicted proteins’ procedure, which compares each predicted protein on genome scaffolds, assembled from metagenomic sequences, to UniProt database, as described previously [11]. When more than 50% of predicted proteins on a scaffold shared the same taxonomic affiliation, which can be on any taxonomic level (from species to kingdom), this scaffold was classified according to the shared taxonomic affiliation [11].

Gene identification and annotation

Open reading frames (ORFs) were predicted using Prodigal (v2.6.3) [18] with the option to run in metagenome mode selected. Sequences of predicted ORFs were annotated using Hidden Markov Models (HMM) [19]. Annotations of the set of genes later analyzed in transcript ratio analyses were also confirmed by aligning against UniProtKB and UNIREF100 databases.

Calculating diametric ratios (DR)

Before analysis RNA reads were trimmed using Cutadapt [20], these RNA reads are available on the short-read archive (SRA) on NCBI, Bio-project ID PRJNA505710. To calculate Diametric Ratios (DR), RNA reads were mapped to the nucleotide sequence of all open reading frames, identified by Prodigal, in each scaffold of the de-replicated genomes sets, of each infant, using Bowtie2.

Further filtering of RNA reads, mapped to each gene sequences, was performed using the script ‘mapped.py’ (https://github.com/christophertbrown/bioscripts/blob/master/ctbBio/mapped.py). This script filters out any reads that map with more than one mismatch to the reference genome.

Transcript abundance was measured as coverage depth, by calculating the number of RNA reads per gene length. To calculate RNA reads coverage depth on each gene sequence we used the script ‘calculate_coverage.py’ (https://github.com/christophertbrown/bioscripts/blob/master/ctbBio/calculate_coverage.py). mapped.py and calculate_coverage.py are part of ctbbio version 0.45.

Gene expression levels were measured using diametric transcript ratios (DR), calculated according to Eq (1): q

Where α is the transcript (RNA) abundance of a certain gene/genes (average abundances if abundances of several genes are considered, as denoted with the overbar) in a specific genome (G), and β is the transcript (RNA) abundance of a different gene/genes in the same genome (G) with an opposite transcriptional response to a surrounding physiological condition (E).

For example, if a certain surrounding physiological condition (E; e.g. oxygen level) increases transcription of gene α in genome G than for the same physiological condition (E) transcription of gene β in genome G decreases (Table 3). Because of the opposite transcriptional responses of examined genes these transcript ratios are referred to as diametric ratios (DR). Coupled genes for DR must have a transcriptional response that is opposite to a given physiological condition (such as norVW vs norR; Table 3) or respond to opposing physiological condition (such as ompC vs ompF Table 3). Calculating transcriptional responses this way allows measuring shifts in transcriptional responses while minimizing biases associated with changes in genome abundance. As noted, transcripts of the two sets of genes (α and β) originate from the same genome (G). Thus, given that the values of α and β are composed from transcript abundance per genome multiplied by genome abundance, the DR formula factors out genome abundance and the DR are comparable between samples.

Table 3

Genes examined in this study and the factors controlling their transcription.

Coupled genes for DR	Genes	Factors controlling transcription§	Variable in Eq (1)	Reference
1)	cydAB	Micro-aerobic↑	α	[21]
1)	cyoABCD	Aerobic↑	β	[21]
2)	arcA	Anaerobic ↑	α	[22]
2)	fnr	Oxygen in-depended	β	[23]
3)	nrdDG	Anaerobic and Micro-aerobic↑	α	[24]
3)	nrdAB	Aerobic↑	β	[24]
4)	norVW	Nitric oxide ↑	α	[25]
4)	norR	Nitric oxide ↓	β	[25]
5)	ompC	High osmolarity ↑	α	[26]
5)	ompF	Low osmolarity ↑	β	[26]

§ Arrow direction indicates whether the controlling factor up or down regulate transcription of the gene.

In each sample of our DR analysis we discarded genes that didn’t have any RNA reads mapped to them after the stringent filtering, to avoid α or β values of 0 leading to extreme DR values of 1 or 0.

Mapping to cytochrome oxidases from infant MAG’s compared to cytochrome oxidases from KEGG

In order to compare between mapping of RNA reads to genes from assembled genomes from premature infant’s gut microbiome to mapping of RNA reads to genes from available databases, E. coli (K-12 MG1655) cytochrome oxidases sequences were downloaded from KEGG (Kyoto Encyclopedia of Genes and Genomes) database. RNA reads were mapped to these genes and filtered as described previously for genes from assembled genomes. Diametric ratios were created for those mappings.

Statistical analyses

Differences in transcript ratios were visualized with boxplots constructed with R software for statistical computing (v3.5.1), using ggplot2 package (v3.1.1) [27].

Differences in DR of examined microbial genes were compared within time blocks between NEC and control infants. The first-time block was from 10th day of life till 19th day of life and a second time block was from 20th day of life till 36th day of life. In addition, differences in the DR of examined microbial genes were compared between all NEC and all control infants. Comparisons were done with Welch's t-test. To avoid type 1 error due to multiple comparisons, p values were adjusted with Bonferroni correction.

Results

Establishment of the diametric ratio as a quantitative metric

In this study, a new approach was examined to circumvent potential biases associated with meta-transcriptomic analysis. This approach involved accurate mapping of RNA reads to metagenomic assembled genomes (MAG’s) from the same pool of samples from which RNA reads were retrieved, to be more confident to the genomic origin of the transcripts. Next, transcript abundance ratios of genes that are known to have opposing transcriptional responses to a specific environmental exposure were calculated. To avoid biases related to changes in genome abundances, ratios were calculated only using transcripts belonging to the same genome within the same sample. This approach was designated as Diametric Ratios (DR) analysis, due to the expected opposite expression patterns of chosen genes. The approaches described here are hypothesis-driven, and are distinct from other recent approaches described for meta-transcriptomic observations, which aim to achieve a wide-spread comprehensive view of the changes in relative abundances of gene families and pathways [6,28].

To assess this approach, we analyzed 20 meta-transcriptomic datasets paired with MAG’s from gut microbiomes of 4 premature infants. Two of those premature infants were eventually diagnosed with NEC.

We focused our analysis to Escherichia spp. as this genus was ubiquitous in almost all of the analyzed samples (as representatives of this genus occurred in most of the samples, allowing statistical analysis between NEC and control samples), while other species or genera were not adequately ubiquitous for conducting comparisons between NEC and control samples (Fig A in S1 File). More than 85% of the predicted proteins on the same scaffolds of each of Escherichia spp. genomes had shared affiliation with Escherichia genus. On the species taxonomic level more than 84% of the predicted proteins on the same scaffolds had shared affiliation with either E. coli or E. vulneris (Table B in S1 File). Relative abundances of different bacterial genera were calculated by the ratio of genome coverage to sum of coverages of all genomes in a sample (Fig A in S1 File). Coverage depth and relative abundances data of other MAG’s in each sample is also available in S3–S6 Files. Each sample included about 15 MAG’s (max-24, min-10), with the genus Escherichia occurring in most samples. It is important to note, that RNA reads from each sample were mapped to the Escherichia sp. genome found in the same sample, as not all of the infants had E.coli species. Subsequent analyses were carried out on the genus level of Escherichia. Therefore, calculated DRs for all Escherichia spp. found in the different samples were gathered and compared between NEC infants and control infants. As Escherichia spp. are facultative anaerobes, they can adjust their life style according to available oxygen, which was well represented in its transcription pattern [3]. Thus, by examining Escherichia spp. we addressed a longstanding idea that inadequate oxygen tension in the intestine, due to reduced regional blood flow (ischemia), may be a key contributor to NEC development [29,30]. In vivo patterns of microbial gene expression in the infant gut may be an indicator of insufficient oxygen supply to the intestine and the progression of NEC. Furthermore, diametric ratios of different sets of genes related to other physiological conditions, such as exposure to different nitric oxide (NO) and osmolarity levels, were examined as well.

We examined transcripts abundances of cydAB and cyoABCD, genes encoding cytochrome oxidases with high and low affinity to oxygen (respectively; see below ‘cydAB cyoABCD diametric ratio’ section). Transcript abundances of each of these two sets of genes (either cydAB or cyoABCD, without any normalization accounting for genome abundance) was extremely variable within all-time blocks of either NEC or control infants. This high variation was exemplified through the high standard deviation compared to the average transcript abundance, yielding high coefficients of variance (between 93–174; Table 4). However, in NEC samples there was a trend showing differences between abundances of cydAB genes compared to cyoABCD genes while in Control infants no differences were found, inspiring the idea of diametric ratios. After DR calculation, much lower variation was found, yielding much lower coefficients of variance (20 and 24 for NEC and control infants, respectively; Table 4).

Table 4

Comparison between transcript abundances and diametric ratio (DR) of Escherichia spp. cytochrome oxidases across all time blocks in NEC and control infants.

Infants	Genes/DR	Average	Standard Deviation	Coefficient of Variance	F-TEST§	T-TEST§
NEC	cydAB	233	217	93.0	1.53E-12	0.008
	cyoABCD	71.7	125	174	2.03E-09	0.117
	cydAB-cyoABCD DR	0.82	0.17	20.3	0.226	2.11E-05
Control	cydAB	2.18	2.73	125
	cyoABCD	3.09	4.40	142
	cydAB-cyoABCD DR	0.43	0.10	23.9

§ p-values

cydAB cyoABCD diametric ratio

To examine the association between microbiome oxygen exposure and NEC development in the gut of premature infants, we first examined genes encoding for cytochrome oxidases. These protein complexes are a part of the electron transport chain that pass electrons to O2 during aerobic respiration. We constructed DR from transcript abundances of bacterial cytochrome oxidase genes that have different affinities for oxygen: cytochrome bd oxidase (cydAB) with high affinity for oxygen and cytochrome o oxidase (cyoABCD) with low affinity for oxygen [31,32]. Consistent with these biochemical predictions, a previous study showed that under microaerophilic conditions there was higher expression of cydAB whereas under aerobic conditions expression of cyoABCD genes increases [21]. Our results showed that Escherichia spp. had higher cydAB to cyoABCD transcript ratio in the gut of NEC infants compared to control infants (Fig 1A), and these differences were significant within both time blocks (p < 0.05) and also across all time blocks (p < 0.01).

Fig 1

Transcriptional response to oxygen by Escherichia spp. in the gut of NEC and control premature infants.

(A) Diametric ratios were compared between NEC and control infants in each time block (short lines above) and across all time points (longer lines). Distributions of diametric ratios were compared using Welch’s t-test with Bonferroni correction. Asterisks and double asterisks (*, **) represent p < 0.05 and p <0.01, respectively. (B) Diametric ratios of cydAB and cyoABCD transcript abundances of RNA reads mapped to gene sequences of Escherichia spp. genomes found in infants’ gut. Filtering of reads with more than 1 miss matches was applied. (C) Diametric ratios of cydAB and cyoABCD transcript abundances of RNA reads mapped to gene sequences of E. coli (K-12 MG1655) downloaded from KEGG (Kyoto Encyclopedia of Genes and Genomes) database. Filtering of reads with more than 1 miss matches was applied. (D) Diametric ratios of cydAB and cyoABCD transcript abundances of RNA reads mapped to gene sequences of Escherichia spp. genomes found in infants’ gut. No filtering of reads with miss matches was applied.

To evaluate whether mapping RNA reads to genes from sample-specific MAG’s could improve results sensitivity compared to mapping to genes retrieved from KEGG databases, we examined results of cydAB and cyoABCD diametric ratios between the two mapping approaches. We found that mapping RNA reads to genes from Escherichia spp. from each infant microbiome MAG’s had more significantly distinguishable DR between NEC and control infants then DR found for RNA reads mapped to genes retrieved from KEGG database (Fig 1A and Fig 1B, respectively). Further evaluation of differences between MAG’s and E.coli K12 gene sequences through alignments of cydA genes showed variable number of mismatches (S8 File), ranging from 1 to 19 in E.coli species and 269 mismatches with E.vulneris. This signal was enhanced when filtering out reads that had more than 1 base mismatch (Fig 1A compared to Fig 1C). After filtering, mapping to KEGG database genes was even less assured as there were less data points due to filtering out of reads that were inaccurately mapped (Fig 1B compared to Fig 1D). These results highlighted the necessity for having MAG’s retrieved from the same sample set as RNA reads were retrieved, as well as accurate mapping of RNA reads.

fnr arcA and nrdGD nrdAB diametric ratios

To further strengthen our observations that Escherichia spp. in the gut of NEC infants were exposed to lower oxygen levels, we analyzed another set of genes that are differentially transcribed in response to oxygen levels, fnr and arcA (Fumarate and nitrate reductase, and Aerobic respiration control protein, respectively). fnr transcription is consistent in both aerobic and anaerobic conditions [22], and the regulatory mode is through changes in FNR protein conformation at different oxygen levels [33]. However, the transcription of arcA, which is regulated by fnr, increases when low oxygen conditions prevail [22].

A third set of genes that their transcription is associated to oxygen levels, nrd, encodes ribonucleotide reductase, which catalyzes the enzymatic reduction of ribonucleotides to deoxyribonucleotides. Escherichia spp. have two sets of nrd genes that are differently expressed depending on prevailing oxygen levels. Under aerobic conditions transcription of nrdAB is upregulated, while under anaerobic condition transcription of nrdDG is upregulated [24].

According to the cydAB and cyoABCD diametric ratios we hypothesized that diametric ratios of arcA and fnr genes and nrdDG and nrdAB genes transcribed by Escherichia spp. in the gut of NEC infants would be higher compared to control infants. Indeed, results of these diametric ratios confirmed the results of cydAB and cyoABCD diametric ratios, showing that higher ratios of both arcA and fnr genes nrdDG and nrdAB genes of Escherichia spp. in the gut of the NEC premature infants were significantly higher, across all time blocks (p < 0.05 and p < 0.01, respectively; Fig 2A and 2B) and specifically in the 1st time block (p < 0.01 and p < 0.05, respectively; Fig 2A and 2B). These results further suggested that Escherichia spp. in the gut of the NEC premature infants were exposed to lower oxygen levels.

Fig 2

Transcriptional response to oxygen, Nitric oxide and osmotic conditions by Escherichia spp. in the gut of NEC and control premature infants.

(A) Diametric ratios were compared between NEC and control infants in each time block (short lines above) and across all time points (longer lines). Distributions of diametric ratios were compared using Welch’s t-test with Bonferroni correction. Asterisks and double asterisks (*, **) represent p < 0.05 and p <0.01, respectively.(A) Diametric ratios of arcA and fnr transcript abundances of RNA reads mapped to sequences of Escherichia spp. genomes found in infants’ gut. (B) Diametric ratios of nrdDG and nrdAB transcript abundances of RNA reads mapped to sequences of Escherichia spp. genomes found in infants’ gut. (C) Diametric ratios of ompC and ompF transcript abundances of RNA reads mapped to gene sequences of Escherichia spp. genomes found in infants’ gut. (D) Diametric ratios of norVW and norV transcript abundances of RNA reads mapped to gene sequences of Escherichia spp. genomes found in infants’ gut.

ompC ompF and norVW norR diametric ratios

The next set of genes we examined encode the Outer Membrane Proteins, ompC and ompF. These genes were previously shown to be transcribed under different osmotic conditions, with high abundance of ompC being associated with inflammatory bowel diseases [34]. High expression of ompC, which has small porin size, occurs during high osmotic conditions, while ompF, which has large porin size, occurs during low osmotic conditions [26]. These genes are reciprocally regulated by ompR, depending on its phosphorylation state [26]. Interestingly, significantly higher diametric ratios of ompC to opmF were found to be transcribed by Escherichia spp. in the gut of NEC premature infants compared to control premature infants across all time blocks (p > 0.05), specifically in the 1st time block (p < 0.05; Fig 2C).

The last set of genes we examined were the norVW genes, coding for NO detoxifying enzyme Nitric Oxide Reductase, and their oppositely transcribed regulating gene norR [25]. Nitric oxide binds to constitutively expressed NorR, which up-regulates the transcription of norVW and down-regulates norR expression [25]. Thus, to assess microbial transcriptional response to NO prior to NEC diagnosis, we measured the ratio of transcript abundances for norVW and norR. Higher diametric ratios of NO detoxifying enzymes compared to norR transcribed by Escherichia spp. were found in the guts of the control infants compared to NEC infants, in both time blocks and across time blocks (p < 0.01, for all cases; Fig 2D). Interestingly, the ratio of norVW to norR transcript abundances was higher at earlier compared to later time points, for both NEC and control infants (Welch’s t test, p = 0.009). This may reflect the response of the preterm infant gut to initial microbial colonization.

Discussion

Using diametric ratios and RNAseq mapping to MAG’s to infer physiological conditions in the gut of premature infants

Here we demonstrate that calculating diametric ratios of genes with opposite transcriptional responses to ambient physiological conditions is an approach that can be effectively used for analyzing meta-transcriptomic data (Table 4; Fig 1). In addition, we show that mapping to sample specific MAG’s provides the most clear and significant signal. Results based on three sets of genes (cytochrome oxidase genes, aerobic/anaerobic regulation genes and ribonucleotide reductase genes) suggest that Escherichia spp. in the gut of two premature infants that developed NEC were exposed to lower oxygen levels than Escherichia spp. in the gut of two premature infants without NEC. Together, these data suggest that hypoxic conditions may exist in the gut prior to NEC development.

Previous animal model studies and analyses of early microbial colonization of premature infant guts showed that at early gestational age the gut milieu was more aerobic [35-38]. During NEC, however, hypoxic conditions were observed in the gut tissue of many patients [39]. Furthermore, histologic examination of removed dead intestine tissue of NEC patients demonstrated coagulation necrosis, evidence for ischemic injury [40], in either small or large intestine. Yet, hypoxia is highly debated as a primary controlling factor of NEC [29,41-43]. Recent theories on NEC development point to the role of gut tissue immaturity, impairing intestinal microcirculation and oxygen delivery [44], which might explain the observed transcriptional response to lower O2 levels in the gut microbiomes of NEC infants (Fig 1A; Fig 2A and 2B). It should be noted also that these circumstances are distinct from those that occur in mature gut systems, where anaerobic conditions are linked to a healthy condition and aerobicity is linked to inflammation [45].

Another potential indicator for progression of inflammatory response in the infant gut is exposure of the gut microbiome to nitric oxide (NO). Increased expression of inducible nitric oxide synthase (iNOS) by host’s gut epithelial cells often occurs as a part of the inflammatory response, and recent studies have suggested that TLR4-mediated iNOS expression is a key element of NEC progression [46]. Thus, diametric ratios for nitric oxide reductase can help infer about exposure of Escherichia spp. to different nitric oxide levels in the gut of premature infants. Counter-intuitively to progression of an inflammatory response in the gut of NEC infants, results of this study indicate down regulation of bacterial genes for NO detoxification in the gut of NEC infants (Fig 2D), indicating lower NO levels in the gut lumen. This might be explained by low oxygen supply to the gut lumen, as suggested by results found by previous examined genes (Fig 1A; Fig 2A and 2B), reducing host epithelial cells iNOS activity to produce the antimicrobial agent NO, as these enzymes need oxygen to produce NO [47,48]. Alternatively, these results might also be explained by the report that norVW transcription decreases with combined oxidative and nitrosative stresses in contrast to nitrosative stress alone [49], as occurs during inflammation [50]. Inflammatory response inducing combined oxidative and nitrosative stresses also stimulates higher cydAB transcription (Fig 1A; [51,52]). An additional physiological condition that can be associated with inflammatory response is altered osmotic conditions [53].

Results of diametric ratios between transcripts of outer membrane protein genes indicate that Escherichia spp. in the gut of NEC premature infants might be exposed to high osmotic conditions (Fig 2C). Consistent with previous observation of high expression of ompC by E. coli during high osmolarity levels and increased adherence to host gut epithelial cell through the development of Crohn’s disease [34]. To the best of our knowledge, little is known about gut lumen osmolarity and the development of NEC. Osmolality of feeding formula and its association with NEC development has been studied, but no clear connection was found [54].

Confounding factors of this study dataset

It should be noted, however, that this data set contains confounding factors limiting the interpretation of Escherichia spp. transcriptional differences between NEC and control infants solely as a result of NEC development. First of all, the limited number of infants (two NEC and two control) examined in this study impaired our ability to draw conclusions. Low sample number also restricted our ability to define cutoffs to optimize DR analysis, as discarding more data points would have limited statistical analysis on the time block level. Secondly, other factors, such as gestational age, mode of delivery and birth weight, also differed between NEC cases and controls (Table 1). These factors were previously shown to affect the development of microbial communities [36,55], suggesting that occurrence of different physiological conditions bring about such alterations in infant gut microbiomes.

Nevertheless, transcription results shown here were opposite than expected according to some of those factors, such as gestational age. According to previous studies, younger gestational age infants might have more aerobic conditions as their microbial communities are associated with a more facultative anaerobic life style compared to older gestational age infants associated with an obligate anaerobic life style [36,38,56]. Whereas our results show that lower oxygen exposure occurred at earlier gestational age, in infants that develop NEC (Table 1; Fig 1A; Fig 2A and 2B). Although, larger sized and higher time point resolution gut microbiome meta-transcriptomic studies would be fundamental to confirm shifts in aerobicity state in the intestine of premature infants. In addition, many of the studies examining how these factors affect microbial communities were done on full term infants and data on gut microbiome of preterm infants is still very scarce. A recent paper examining gut microbiome in preterm infants showed that cesarean or vaginal birth mode did not significantly affect microbial communities [57], unlike full-term infant where delivery mode is a major factor shaping infant gut microbiome [55,58].

Although it is hard to conclude whether observed differences were associated with NEC development or other factors, the methodological approach described here still shows a clear and significant signal distinguishing between NEC and control premature infants. A larger study is needed to verify these results and confirm that the gut microbiome in early stages of NEC development senses and responds to different physiological conditions compared to microbiomes in guts of premature infant where NEC is not developing. Further experiments can verify this approach, either by experiments where stool samples are inoculated into artificial media with varying physiological condition (e.g. oxygen levels or NO) or experiments with animal models inducing wanted physiological conditions in-vivo.

Concluding remarks on examined approaches

The approach described here can add new insights into gut microbiome analysis. It is important to note that this approach relies upon prior knowledge on the transcriptional responses of different genes to physiological conditions. It is, thus, essential to do a preliminary literature survey on gene expression pathways of dominant species in examined samples to decide on genomes and genes from which to calculate diametric ratios. In addition, a significant factor affecting the accuracy of DR measurements is proper read mapping. Mapping RNA reads to sample-specific MAG’s is shown here to give more significantly distinguishable signal compared to mapping to genes retrieved from an outside database (Fig 1). The mapping exercise shown here might potentially describe the extreme end of strain heterogeneity, as K-12 E. coli is a laboratory strain that can be different enough from gut E. coli to result in inquorate RNA read mapping. Choosing genomes from databases originating from gut microbiome databases might be more adequate for RNA mapping than KEGG database, as long as chosen genomes are most similar to those genomes for which RNA was sequenced. Genome databases are valuable for taxonomic identification, as done here for confident identification of MAG’s in our samples. However, based on propositions put forth in this study, further research might further enforce the idea that MAG’s enable more accurate mapping and promote better quantification of RNA reads than genomes from databases. It is also important to note that potentially fewer metagenomes are required per individual than meta-transcriptomes, as it was previously shown that specific strains can remain stable within an individual human host [12,59].

In conclusion, we show here how meta-transcriptomic data combined with sample-specific MAG’s can be applied effectively to probe the physiological conditions that gut microbial communities experience by comparing diametric ratios. Further application of this approach can bring new insights on microbe-host interactions within the GI tract systems, and potentially help identify biomarkers for early detection gut diseases, such as NEC, onset and progression.

Fig A. Relative abundances of different Bacterial genera; Table A. Accession phrases for Escherichia spp. genes in on ggkbase interface; Table B. Percent of shared affiliation of the predicted proteins.

(DOCX)

Click here for additional data file.

Genomic relative abundances and coverage depths.

(XLSX)

Click here for additional data file.

Scaffold to bin file, infant 64.

(TSV)

Click here for additional data file.

Scaffold to bin file, infant 66.

(TSV)

Click here for additional data file.

Scaffold to bin file, infant 69.

(TSV)

Click here for additional data file.

Scaffold to bin file, infant 71.

(TSV)

Click here for additional data file.

Diametric ratios.

(XLSX)

Click here for additional data file.

CydA-Alignments.

(TXT)

Click here for additional data file.

3 Sep 2019

PONE-D-19-22108

Combined analysis of microbial metagenomic and metatranscriptomic sequencing data to assess in situ physiological conditions in the premature infant gut

PLOS ONE

Dear Dr. Sher,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Your manuscript has been evaluated by two reviewers. As you can see in their comments, they have numerous comments concerning clarity of the study and its presentation. Hence, submitting a revision might be challenging. However, since I agree with both reviewers that the study is interesting, I give you the opportunity to submit a revised version of the manuscript. Please make sure that all concerns raised by the reviewers are well addressed in the revision.

We would appreciate receiving your revised manuscript by 20 October 2019. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Erwin G Zoetendal, PhD

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

1. Thank you for including your competing interests statement; "The authors have declared that no competing interests exist."

We note that one or more of the authors are employed by a commercial company: Enview, inc.

Please provide an amended Funding Statement declaring this commercial affiliation, as well as a statement regarding the Role of Funders in your study. If the funding organization did not play a role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript and only provided financial support in the form of authors' salaries and/or research materials, please review your statements relating to the author contributions, and ensure you have specifically and accurately indicated the role(s) that these authors had in your study. You can update author roles in the Author Contributions section of the online submission form.

Please also include the following statement within your amended Funding Statement.

“The funder provided support in the form of salaries for authors [insert relevant initials], but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The specific roles of these authors are articulated in the ‘author contributions’ section.”

If your commercial affiliation did play a role in your study, please state and explain this role within your updated Funding Statement.

2. Please also provide an updated Competing Interests Statement declaring this commercial affiliation along with any other relevant declarations relating to employment, consultancy, patents, products in development, or marketed products, etc.

Within your Competing Interests Statement, please confirm that this commercial affiliation does not alter your adherence to all PLOS ONE policies on sharing data and materials by including the following statement: "This does not alter our adherence to  PLOS ONE policies on sharing data and materials.” (as detailed online in our guide for authors http://journals.plos.org/plosone/s/competing-interests) . If this adherence statement is not accurate and  there are restrictions on sharing of data and/or materials, please state these. Please note that we cannot proceed with consideration of your article until this information has been declared.

Please include both an updated Funding Statement and Competing Interests Statement in your cover letter. We will change the online submission form on your behalf.

Please know it is PLOS ONE policy for corresponding authors to declare, on behalf of all authors, all potential competing interests for the purposes of transparency. PLOS defines a competing interest as anything that interferes with, or could reasonably be perceived as interfering with, the full and objective presentation, peer review, editorial decision-making, or publication of research or non-research articles submitted to one of the journals. Competing interests can be financial or non-financial, professional, or personal. Competing interests can arise in relationship to an organization or another person. Please follow this link to our website for more details on competing interests: http://journals.plos.org/plosone/s/competing-interests

3. Please include captions for your Supporting Information files at the end of your manuscript, and update any in-text citations to match accordingly. Please see our Supporting Information guidelines for more information: http://journals.plos.org/plosone/s/supporting-information.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: No

Reviewer #2: I Don't Know

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: No

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: Yonathan Sher et al describe a method to study physiological conditions based on bacterial gene expression, applied on feces derived from preterm infants that do or don’t develop NEC. While the approach in itself is interesting and could be applied in multiple settings, the method is not described in enough detail to be repeated by others, it is not validated under known ‘physiological conditions’ and I doubt its applicability in diagnostic setting as proposed by the authors. Since NEC most often occurs after 2 to 4 weeks of life, time is a limiting factor when using such methodology for early NEC diagnosis. See my detailed comments below.

[General]

- Many grammatical and spelling errors throughout the manuscript

- Is oxygen deprivation in the gut of preterm infants a local issue, or does it affect the complete intestine? Does it affect both the small and large intestine? What was the situation for the two included NEC infants? If affecting the small intestine in a local manner, are feces samples a good representative?

[Introduction]

Line 36: remove dash for in-vivo.

Line 64: ‘can be used to the study’ change to ‘can be used to study’.

[Methods]

Used software is without version numbers and parameters/settings and can therefore not be reproduced by others.

Lines 69-71: I encourage a more detailed description of how fecal samples were processed and stored to indicate compatibility for RNA sequencing. In addition, the manuscripts focusses on oxygen exposure, which may be challenging taking the method of fecal sample collection into consideration.

Lines 88: Change ‘which use’ to ‘which uses’.

Lines: 92-95: Does this mean these scaffolds are not yet made available? Please make them available with the next submission.

Line 99: Change ‘phyla level’ to ‘phylum level’.

Line 103: Write HMMs fully, since it is the first and only time the word is used.

Line 109-110: Please clarify what is filtered for by mapped.py. How does this complement to the mapping performed using Bowtie2?

Line 112: Use of "coverage" is ambiguous: breadth or depth? (see also https://www.danielecook.com/calculate-depth-and-breadth-of-coverage-from-a-bam-file/)

Line 116-120: This equation only takes transcript abundances into account. I wonder how this corrects for genome abundance. Did the authors considered to correct for genome abundance by using their metagenomics data? Or by other means (by correcting for total reads or total mapped reads (https://www.biostars.org/p/273537/). Also, please clarify how opposite regulating genes are identified. What does opposite regulating genes mean exactly and how can one find them for their function of interest?

[Results]

Write the section in past tense

Line 140: Change to ‘In this study new approaches were examined to…’

Line 142: Change to ‘from the same pool of samples from which RNA reads were retrieved…’

Line 143: Change ‘confidante’ to ‘confident’

Lines 140-150: Many grammatical errors, please correct this.

Line 153: The term NEC was already introduced, so the abbreviation can be used.

Lines 155-164: Move to introduction section.

Figure 1: Change the order of A,B,C and D, since now the authors first refer to figure 1D, and later on to A B and C.

Line 177: The authors describe a significant difference in cydAB to cyoABCD ratio across all time points. I wonder how statistical analysis on the separate time points has been performed, since collection of more than 1 sample on each time point was rare, as shown in table 2, and statistical analysis can therefore not be performed. Maybe the authors meant to say something different, please clarify.

Line 184: Authors refer to figure 1CD, while to understand the differences between using KEGG database of infant MAGs, one should compare A with C.

Lines 182-184: Although not significant, the same trend can be observed when using genes from the KEGG database. Could this be a matter of small sample size?

Lines 187-188: It could also be a matter of picking the right E. coli strain from the database. Why was E. coli K12 selected? This is a laboratory strain and may not reflect functional properties (and therefore genes) that are needed to thrive in the human gut. A better option may be to use public available whole genomes from clinical isolates from the NICU.

Line 193: Add full name of the genes

Line 194: Change FNR to fnr

Line 202: Change to ‘hypothesized’.

[Discussion]

Line 298-299: Did the authors consider to mimic low and high oxygen exposure in an in vitro setting or animal model to confirm their method? Would rather do that to begin with than studying a bigger group of infants. Currently, the method is not validated, but directly applied to clinical samples during which too many aspects are based on assumptions.

Reviewer #2: This paper describes an interesting approach to make more (or better) use of meta-transcriptomics data. Even though a priori knowledge is needed, thus limiting it to exploration of already established pathways and transcriptional relationships.

This idea to make use of the metagenomic derived MAG's poses an interesting approach. However, the described methodology is, in the current manuscript text, not always clear enough to follow.

General remarks:

- The tone of the abstract is rather strong. The main text is better nuanced, especially due to the section on the limitations (i.e. the Confounding factors section that describes the sample size, subject "baseline" differences, etc). Strongly suggest to tone down the abstract to prevent overselling.

- The definition and formula of the diametric ratio remain a bit vague. Though it is understandably depending on each gene, it would help to literally state the coupled genes (Table 3 does not provide direct linking). Furthermore, actual descriptive statistics on the transcript abundances of the used examples would greatly help following the method (so the A's and B's of each of the example gene sets, as the boxplots only show the final outcome ratio). Also, any given external condition/stimulus would likely affect several genes, what to do with other genes that are known to be affected by oxygen or infllammation status?

- Although reasoning to focus on E.coli is understandable, it does make the reader wonder what happened to all the other MAGs. How many were there per sample?

- In line with the previous point how well are the MAGs assigned to a taxonomic unit? How confident is the assignment to E.coli? In the Materials and Methods section the phylogenetic assignment appears to been done at the very high phylum level. How deep were the samples sequenced (both transcriptome and metagenome)? Especially since the rather short read HiSeq platform was used, a reasonable depth would be needed. These technical details, that are needed to assess the quality of the reported analysis, are missing or should be reported directly into the main text.

- Though conceptually it could make sense, the conclusion that mapping RNA reads to sample-specific MAGs is more significantly distinguishable compared to mapping to genes retrieved from an outside database is hard to follow in the text. From the results it is not particularly clear if the mentioned E.coli focus means only one (or more) found MAG(s) that was/were assigned to E.coli was/were used or not. Assuming that this is a (set of) MAG(s) identified in these samples, then the comparison appears to be made against KEGG references genes, but these only originate from one strain (K-12 MG1655)? Is the MAG assignment considered accurate to strain level (see earlier point)? If not, then why not compare to various, if not all, E.coli genes in the online repositories? Especially for E.coli there is a huge genetic variation between all the known strains. Comparing to a single set of genes can therefore not prove this is a better method.

Text is mostly well readable though some grammar checks may bring some improvements - some specific examples:

- line 196 "higher when there are more anaerobic conditions [24].": Did you mean "higher when the condition is more anerobic [24]". Perhaps it is even better to state "less/lower oxigen".

- line 236 'provide" --> provides ?

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

18 Oct 2019

We have addressed each concern raised by reviewers, in the 'response to reviewers letter', specifically addressing comments concerning clarity of the study and its presentation.

Submitted filename: Response to Reviewers_v2.docx

Click here for additional data file.

6 Nov 2019

PONE-D-19-22108R1

Combined analysis of microbial metagenomic and metatranscriptomic sequencing data to assess in situ physiological conditions in the premature infant gut

PLOS ONE

Dear Dr. Sher,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Your manuscript was evaluated by the same reviewers. Although they noticed some improvements, both of them were not convinced about the validity of the approach used, to which I agree. This is a major concern with the study and I was temped to reject the manuscript for publication in PLOS ONE. However, I decided a major revision as reviewer 2 pointed out that the approach could maybe be useful for others. Hence, I give your an opportunity to submit a revision. I recommend you to address all comments raised by the reviewers appropriately in a revised version. Make sure you convince us of the validity of the approach used or alternatively discuss in detail the drawbacks of it and tune down the conclusions.

We would appreciate receiving your revised manuscript by 10 December 2019. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Erwin G Zoetendal, PhD

Academic Editor

PLOS ONE

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: (No Response)

Reviewer #2: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: N/A

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

Reviewer #2: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: Yonathan Sher et al made changes to the manuscript accordingly. I, however, remain hesitant about the validity of the presented approach and drawn conclusions, because of unclarities on various levels (sample type, data processing, DR criteria, data interpretation and no validation). See detailed comments below.

[General]

Is a minimum coverage depth required for the DR to be ‘accurate’? If so, what is suggested? If not, wouldn’t that be required?

[Abstract]

Line 21-22: Organism abundance is actually taken out of the equation. Authors don’t ‘correct’ for organism abundance, because the DR is not affected by organism abundance. I suggest to rephrase this to make it more clear to the readers.

[Introduction]

Line 40-41: Suggestion to change to “Physiological conditions within the gut are important metrics to measure when studying gut inflammatory diseases [1,2], yet are notoriously difficult to measure in vivo.

Line 66-67: Based on the analysis described herein, the conclusion that the method is widely applicable cannot be drawn.

[Methods]

Line 73: Suggestion to say ‘As previously described…’ instead of ‘as noted previously…’

Line 78-79: Clarify on which basis the time blocks were stratified. Sample number-driven (so each block has 5 samples) or ‘biology’-driven, or other?

Table 1: Did any of these infants receive breathing support? Continuous positive airway pressure can result in an increased oxygen supply to the GI tract.

Line 110: This is just a link to the website, how can I find these particular scaffolds? Please add a accession number or something similar.

Line 102: In this section, it is unclear if it is about how RNA OR DNA sequences are processed. I suppose genome reconstruction is based on metagenomics data?, but in the first sentence (Line 103) RNA is mentioned as well, as well as in the later part of this section (Lines 108-111).

Line 113: Same here, I assume metagenomics data is used for phylogenetic profiling, but it is important to mention it somewhere nevertheless.

Line 125: Wouldn’t you expect to obtain more MAG’s from feces? I understand these are samples from preterms in early life, and therefore the microbiota is still minimal and developing. However, looking at the relative abundance profiles I would expect to identify more species than 15. Is it a methodological issue?

Line 125: Could the fact that most MAGs were Escherichia be a result of availability (redundancy) in the database?

Line 133: Change to “To calculate Diametric Ratios (DR), RNA reads were…”

Table 3: The table shows different ‘relationships’ between the coupled genes for DR. Some are opposing in a specific condition (e.g. number 4 in response to NO), while some conditions are opposing which results in transcription of different genes (e.g. number 5). Please clarify on which terms genes are considered opposite transcriptional responses. With other words, which criteria need to be fulfilled for genes to be identified as ‘opposite transcriptional responses’. For example, why are arcA and cyoABCD no opposing genes, as one goes up in anaerobic conditions, and the other goes up in aerobic conditions (which is a similar difference in condition as number 5 for osmolarity).

Line 173-174: Please add the version of R and the ggplot2 package.

Line 175-178: I remain unsure about which comparisons are made. Control versus NEC at time block one; Control versus NEC at time block two; and control versus NEC of all samples? I am confused because the authors write “within time block”, which I interpret as comparing samples within a time block with each other, instead of comparing all samples in a time block between control and NEC.

[Results]

Line 181: ‘a new approach was examined to..’

Line 182-184: How can this be an approach by itself? An approach to achieve what exactly? Which potential bias are circumvented by this? In the introduction, the authors focus on one new method: the diametric ratio. To my opinion, mapping of RNA reads to MAG’s is just part of the process to calculate the DR. However, two methods for mapping are examined: mapping to MAGs and mapping to the genome of E.coli K-12 from database. I suggest to focus on the DR approach and the comparison of the two mapping methods.

193: ‘to assess this approach’

Line 196: When is a genus considered adequately ubiquitous?

Line 203: The fecal microbiota does not represent the microbiota lining the mucosal layer, which may be more affected by reduced regional blood flow than luminal microbiota. In addition, the gut is not oxygen free in early life, which is important to take into consideration, particularly for the early timepoint samples (~first month, probably longer in preterm infants).

Line 211-212: ‘Transcript abundances of each of these two sets of genes (either cydAB or

212 cyoABCD, without any normalization accounting for genome abundance) was extremely variable’. Variable between what? Between samples, between NEC/control?

Line 214: Change to (between 93-174) – 93 instead of 92.

Table 4: Please add the DR for all gene sets as supplementary data.

Table 4: Average coverage of the genes in control infants was only 2-3. Getting back to my earlier comment: “Is a minimum coverage depth required for the DR to be ‘accurate’?”

Line 233-236: Escherichia spp were of much lower abundance in the control than NEC infants (Fig S1). May it be that Escherichia of low abundance express a different set of genes than dominant Escherichia, simply as a result of their abundance? Overall metabolic activity may, independent of oxygen levels, already be different. It looks like the gut environment of preterms developing NEC favors conditions for Escherichia growth and replication. Could a more competitive environment for Escherichia (like in control infants) alter expression of cydAB cyoABCD genes, since other activities may be of more importance for survival? Can the authors check for activity status of Escherichia, maybe by checking the expression of housekeeping genes, or by RNA coverage over the entire genome? Long story short: How sure are the authors that the difference in DR represents differences in oxygen levels at the epithelial lining of the gut, and not overall deprived activity in case of low abundance?

Line 258-260: Since this must be a result of RNA mapping efficiency, can the authors give an impression of the difference in gene sequence (of the specific genes of interest) between the MAGs and E coli K12?

[Discussion]

Line 330-331: On which criteria is the conclusion of ‘effectively used’ drawn? When is a method effective?

Line 412: affecting

Line 417-419: Would it be appropriate to screen the databases for Escherichia spp that result in the highest transcript abundance, and use this/these for analysis? In addition, would it be appropriate if the ‘best pick’ from the database is difference for each sample?

Line 419-422: Rebuild sentence.

Line 429: Suggestion to remove “and potentially help identify biomarkers for early detection of NEC progression in premature infants.”

Reviewer #2: The authors have greatly increased readability/transparency of the analysis. Whether or not this approach is practically suitable in preterm NEC setting is up to debate, but the approach to make use of metatranscriptomics in the described method is interesing and now much clearer for other parties to repeat.

Weaknesses and assumptions are now much more transparent.

Minor comment(s):

- Lines 107-108 "One genome of each 107 identified bacterial species was manually chosen from each infant gut microbiome for downstream RNA analysis" --> I understand the manual choice here, but please clarify where this information (which bacterial species chosen for which smaples) can be found. Online? SI Table 1?

- In the Materials and Methods sections some "results" are presented among the methodology (such as the sentences in line 118-121; line 124-125). These would fit better in the Result section.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

10 Dec 2019

We have fully responded to of the reviewer and editor comments, as addressed in response letter.

Submitted filename: 2nd revision Response to Review Comments.docx

Click here for additional data file.

8 Jan 2020

PONE-D-19-22108R2

Combined analysis of microbial metagenomic and metatranscriptomic sequencing data to assess in situ physiological conditions in the premature infant gut

PLOS ONE

Dear Dr. Sher,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

The reviewers of the previous versions indicated the manuscript was majorly improved, but still had a couple of minor issues that need to be addressed. Please address these.

We would appreciate receiving your revised manuscript by 20 February 2020. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Erwin G Zoetendal, PhD

Academic Editor

PLOS ONE

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: (No Response)

Reviewer #2: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: Yonathan Sher et al made changes to the manuscript accordingly, including clarification of applied data processing methods and DR criteria. I remain hesitant about whether the analysis can be repeated by another party (due to the arbitrary nature of the approach) and about its applicability in clinical setting. However, the proposed method may inspire other researchers to think in this direction, and to be aware of potential biases, when analysing their -omics data.

Minor comments:

- Check grammar and spelling once more, carefully, throughout. E.g. lines 197-204: Mixed use of present and past tense (was examined, involves, were retrieved etc), same accounts for lines 214-218 (focused, was, occur, were not).

- Based on the analysis described herein, the conclusion that the method is widely applicable

cannot be drawn (the sentence stating this has already been removed in the introduction section, but not yet in the abstract).

- I suggest, for the conclusion, to focus on the approach (like lines 443-446), and not so much on the preterms/NEC. E.g. In conclusion, we show here how meta-transcriptomic data combined with sample-specific MAG’s can be applied effectively to probe the physiological conditions that gut microbial communities experience by comparing diametric ratios. Further application of this approach can bring new insights on microbe-host interactions within the GI tract systems, and potentially help identify biomarkers for disease onset/progression.

Reviewer #2: The Manuscript has been substantially improved. Validation of the approach remains a weaker point but this is much more addressed and the manuscript invites others to further explore and validate.

Some minor items to consider:

DR definition:

- The definition is based on the "opposite transcriptional response" of genes to a given environmental condition. Should this be strictly kept at transcription? Or could the "opposite" term also be "different"? Because the fnr/arcA example describes that arcA react to the condition, but the fnr transcription is (expected to be) constant so fnr transcription does not rely on the environmental condition (it's translated protein FNR does react to the condition). All other gene sets are obviously of a much more "opposite" nature.

Methodology - suggestion to enhance repeatability of procedure by others

- in the Results (line 208-210) it is stated that RNA reads were mapped on the MAGs per sample, but later this was combined per genus level. To ensure others can repeat your approach you should state how this "combining" was done. Are mapped results/data further "aggregated" based on the gene name/annotation of the different Escherichia species/MAGs?

Perhaps this is evident from the SI files, but it would make it conceptually easier to understand for a larger audience.

Suggestion to strengthen discussion in "Confounding factors" section

- the paragraph with lines 402-412 states that the presented findings are opposite to published findings related to gestational age. Though I agree that this intuitively might seem to be the case: the figures of the DRs do not seem to show a time development in either Control or NEC group, except for norVW/norV DR but this indeed seems to go in the opposite direction. However, as the breathing support is not (properly) recorded and this might have an impact. Moreover, DR findings are "relative" and no true "later gestational age / confirmed anaerobic" infant fecal sample is analyzed and such a validation could give more insight.

Some minor language items I noticed:

- line 40 "in vivo" should be in italics

- I believe "spp." should not be in italics (many times throughout the text and figure legends)

- suggestion to add the word "to" (or "with") in line 405 after "compared"

- line 407: "occur" should be "occurs"?

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

21 Jan 2020

Thank you for indicating the major improvements of the manuscript.

We fully addressed all the minor issues pointed by the reviewers.

Submitted filename: 3rd revision Response to Review Comments_c.docx

Click here for additional data file.

10 Feb 2020

Combined analysis of microbial metagenomic and metatranscriptomic sequencing data to assess in situ physiological conditions in the premature infant gut

PONE-D-19-22108R3

Dear Dr. Sher,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

I saw that (an)aerobic (for microbes) and (an)oxic (for conditions) may have been mixed up occassionally, such as for example in line 359. Make sure these terms are correctly used when checking the proof.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Erwin G Zoetendal, PhD

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

21 Feb 2020

PONE-D-19-22108R3

Combined analysis of microbial metagenomic and metatranscriptomic sequencing data to assess in situ physiological conditions in the premature infant gut

Dear Dr. Sher:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Erwin G Zoetendal

Academic Editor

PLOS ONE