Lettuce (Lactuca sativa L.) is one of the most common crops in the world, and is also widely consumed in South Korea. Ever since its introduction into Korea at the 6-7th century AD, lettuce has remained as one of the major crops. Nowadays, South Korea produces more than 200,000 tons of lettuce in over 8,000 ha of land, and production and consumption are steadily rising [1]. Lettuce is not only rich in vitamins and potassium, but also contains dietary fiber and minerals, making it a very nutritious food source [2]. However, the consumption of raw lettuce has been closely linked to food poisoning. Although the route of contamination in lettuce has not been established, there were reports of hemolytic-uremic syndrome and bloody diarrhea caused by enterohemorrhagic Escherichia coli serotype O157:H7 contamination [3, 4]. In addition, Aeromonas spp. and Yersinia enterocolitica, which can cause food poisoning, were discovered in packaged lettuce [5]. Furthermore, it has been shown that foodborne pathogens such as Bacillus cereus, Listeria monocytogenes, Salmonella enterica serotype Typhimurium, Staphylococcus aureus, and Campylobacter jejuni cause food poisoning through contaminated lettuce [6-8].

Culture-dependent methods have been traditionally used to detect and identify foodborne pathogens, but there are also obvious limitations [9]. This type of analysis is time-consuming and laborious. In addition, all taxa cannot be identified from total bacteria in the samples. Therefore, 16S rRNA gene-based sequencing has been widely used to explore microbial diversity, including foodborne pathogens. Next-generation sequencing of DNA extracted from environmental samples and metagenomic databases packed with sequence information from diverse habitats are useful for this non-cultivation approach [10, 11]. Although microbial diversity within several plants has been reported [12], few studies have examined the microbial community and potential pathogens in raw food, including lettuce.

Microorganisms live in various parts of plants and interact with them. For example, the root microbiota play a very important role in the growth and health of the host, and more than 10⁵ cell/g bacteria can colonize the surface of plants [13]. Interactions between the plant microbiota and the host continue to be revealed by 16S rRNA gene-based analysis [10]. It has been reported that plant energy metabolism is also mediated by the phyllosphere microbiota [14]. Another study revealed that the phyllosphere microbiota are affected by cross-talks with the rhizosphere microbiota [15]. Therefore, understanding the correlation between the phyllosphere microbiota in lettuce is important to prevent food poisoning caused by contaminated vegetables [16].

In this study, we conducted 16S rRNA gene-based analysis to investigate the microbiota of lettuce collected in South Korea at different seasons. We obtained lettuce from five sampling sites (n = 3 at each site) known for their highest lettuce production in South Korea, and examined their microbiota. Samples were collected at each site between April and July to determine the effect of temperature on the microbiota. 16S rRNA gene sequencing of the microbiota on lettuce was performed, and region-and season-specific differences were examined. In addition, core bacteria in the microbiota were selected, and their correlations were calculated and analyzed through network generation [17]. The results of the present study should extend our understanding of food poisoning by raw lettuce and help us better manage fresh vegetables.

Sample Preparation

We collected 30 samples (100 g of each bundle) from five different sites in South Korea between April (relatively low temperature, monthly maximum average temperature of 21.5 ± 0.6°C, monthly minimum average temperature of 8.0 ± 0.8°C) and July (relatively high temperature, monthly maximum average temperature of 30.6 ± 0.4°C, monthly minimum average temperature of 21.7 ± 0.7°C) of 2016 (Fig. 1). According to information from the National Statistical Office, four sites in the Gyeonggi-do area (31,263 tons in 2015) and one site in the Jeon-buk area (12,645 tons in 2015) were selected on the basis of lettuce yield. Three samples were randomly harvested at each site using sterilized gloves and plastic bags. Samples were packed with ice in a box and transported to the laboratory. Analysis was performed on the same day. Samples were measured out (25 g) into filter bags (FILTRA-BAG; Labplas, Canada), designed to simplify taking an aliquot and has 1,840 holes per square inch with a pore size of approximately 330 microns, and were mixed with 225 ml of Buffered Peptone Water (OXOID, UK). Using a BagMixer 400 W (Interscience, France), the mixtures were homogenized to detach bacterial cells from the plant. The mixtures were filtered and centrifuged at 10,000 ×g for 10 min at 4°C. After the supernatant was discarded, the residual pellet was washed with 5 ml of TES buffer (0.1 M NaCl, 10 mM Tris-HCl, pH 8.0, 1 mM ethylenediaminetetraacetic acid (EDTA)) to remove impurities such as plant tissues. The pellet was repeatedly washed and centrifuged two more times; pellets were stored at -80°C prior to metagenomic DNA extraction.

Figure 1. Map of the sampling sites. All samples were obtained from five different sites in April and July 2016. Maximum and minimum of average temperatures are noted at the bottom.

Extraction Procedure of Metagenomic DNA

Metagenomic DNA was extracted via a slightly modified DNA extraction method described in a previous study [18]. The pellet was suspended in 500μl of PVP/CTAB buffer (1%polyvinylpyrrolidone, cetyltrimethylammonium bromide) with 50 μl of lysozyme solution (100 mg/ml) to remove polyphenols that cause enzymatic browning in lysates. Each pellet was incubated at 37°C for 1 h and frozen at -80°C for 10 min. For thawing, the mixture was incubated at 37°C for 10 min. The DNA pellet was then treated with a proteinase K mixture (140 μl of 0.5 M EDTA, 20 μl of 20 mg/ml Proteinase K, 40 μl of 10% sodium dodecyl sulfate) and incubated at 56°C for 1 h. This was followed by centrifugation at 21,206 ×g for 1 min, after which the supernatant was transferred to a microtube and treated with 5 M NaCl and 80 μl of CTAB/NaCl solution. The supernatant solution was thoroughly mixed with one volume of phenol/chloroform/isoamyl alcohol (25:24:1 (v/v/v)), and was centrifuged at 21,206 ×g for 5 min at room temperature (RT). The upper phase was transferred to a new tube and treated with the same volume of chloroform. After re-centrifuging the upper phase at 21,206 ×g for 5 min at RT, RNase A (100 mg/ml) was added and the mixture was incubated at 37°C for 1 h. Phenol-chloroform extraction was conducted again to remove the RNase A. The DNA mixture was then treated with 10% volume of 3 M sodium acetate (pH 5.0) and two volumes of ice-cold 100% ethanol, and was centrifuged at 21,055 ×g for 20 min at 4°C. The supernatant was discarded and ice-cold 70% ethanol was added for washing. The resultant solution was centrifuged at 21,055 ×g for 5 min at 4°C; the supernatant was discarded. After drying out the pellet at RT, DNA was resuspended with 50 μl of TE buffer and incubated at 55°C for 1 h. Genomic DNA was quantified using Optizen Nano Q (Mecasys, Korea) and was stored at -20°C before amplification.

PCR for Specific Bacterial Detection

Primers used in this study for pathogen detection are listed in Table 1 [19-24]. Amplification of the lettuce samples was carried out in a total volume of 50 μl as follows: 1 μl of extracted DNA template (32.3-119.5 ng), 1 μl of (0.2 μM) primers, 5 μl of dNTP (0.2 mM), 5 μl of 10× nTaq Buffer, 1 μl of nTaq DNA Polymerase (Enzynomics, Korea), and 36 μl of distilled water. The PCR cycling parameters were as follows: 95°C for 5 min, and then 30 cycles of denaturation at 95°C for 50 sec, annealing at 55-61°C for 20-60 sec, extension at 72°C for 2 min, and final extension at 72°C for 6 min. Aliquots (20 μl) of the PCR products were analyzed by electrophoresis on 1.2% and 2% (w/v) agarose gels (SeaKem LE Agarose; Takara, Japan).

Table 1 . PCR primers used in this study..

Target bacteria	Target gene (name)	Primer sequence (5’ → 3’)	Amplicon length (bp)	Reference
Bacillus spp. (B. anthracis, B. cereus, B. mycoides, B. pseudomycoides, B. thuringiensis and B. weihenstephanensis)	nheA	F : TAAGGAGGGGCAAACAGAAG R : TGAATGCGAAGAGCTGCTTC	500	[19]
	nheB	F : CAAGCTCCAGTTCATGCGG R : GATCCCATTGTGTACCATTG	770	[19]
	nheC	F : ACATCCTTTTGCAGCAGAAC R : CCACCAGCAATGACCATATC	582	[19]
	motB	F : CGCCTCGTTGGATGACG R : GATATACATTCACTTGACTAATACCG	280	[20]
	Hemolysin	F : CTGTAGCGAATCGTACGTATC R : TACTGCTCCAGCCACATTAC	185	[21]
Enterococcus casseliflavus	vanC2/C3	F : CTCCTACGATTCTCTTG R : CGAGCAAGACCTTTAAG	439	[22]
Klebsiella pneumoniae	rcsA	F : GGATATCTGACCAGTCGG R : GGGTTTTGCGTAATGATCTG	176	[23]
Pseudomonas aeruginosa	V2/V8	F : GGGGGATCTTCGGACCTCA R : TCCTTAGAGTGCCCACCCG	956	[24]
Total bacteria	16S rRNA V5/V6 (799Fmod6/1114R)	F : CMGGATTAGATACCCKGGT R : GGGTTGCGCTCGTTGC	382	[25, 26]

MiSeq Platform for Analysis

To conduct next-generation sequencing via Illumina MiSeq sequencing (Illumina, USA), 16S rRNA genes were amplified using V5-V6 region-specific primers (Table 1) [25, 26], with pre-adapter and sequencing primers (overhang adapter) designed to minimize the amplification of chloroplast and mitochondria 16S rRNAs attached [25]. The PCR products were purified using the MEGAquick-spin Plus kit (iNtRON, Korea). Index PCR was performed using the Illumina Nextera XT Index kit, and the library was purified with AMPure XP beads (Beckman Coulter, USA) according to the manufacturer’s instructions. The size and quality of the library were validated using the Agilent Bioanalyzer 1000 chip and the KAPA qPCR kit. Equal amounts of libraries from all samples were pooled, and paired-end (2 × 300 bp) sequencing was conducted using the MiSeq system (Illumina, USA) at Macrogen (Korea), based on the manufacturer’s instructions.

Bacterial Load Quantification Using Quantitative Real-Time Polymerase Chain Reaction

Bacterial loads obtained from lettuce were analyzed according to previous descriptions with slight modification [27]. We quantified the 16S rRNA gene in lettuce by using real-time PCR to determine total bacterial loads. Real-time PCR was conducted using a primer set that targeted the same region as that used to conduct MiSeq sequencing with the CFX Connect Optics Module (Bio-Rad, USA). Reaction mixtures (20 μl) containing 10 μl of SSoAdvanced Universal SYBR Green Supermix (2×) (Bio-Rad, USA), 10 μM of each primer, and 1 μl of DNA template (10-fold dilution series of sample DNA) or distilled water (negative control) were prepared, and samples were run in triplicates. The cycling parameters for real-time PCR were as follows: 98°C for 3 min, and then 40 cycles of denaturation at 98°C for 15 sec, annealing at 55°C for 30 sec, extension at 72°C for 30 sec, and final extension at 72°C for 1 min. The serial log-concentration of E. coli genomic DNA was used to generate 16S rRNA standard curves [28]. Bacterial loads in the sample were calculated by comparing the Ct value with the standard curve. The regression coefficient (r2) for the standard curve was 0.9959.

Data Analysis

Raw sequences obtained were trimmed and chimera-checked using the MOTHUR software. The merged reads, where barcode and primer sequences were removed, were trimmed on the basis of average quality score (<25) and the number of homopolymers (>8). The UCHIME function in the MOTHUR software was used to remove chimeric sequences [29]. Determination of operational taxonomic units (OTUs) was processed on the CLC Genomics Workbench (CLC Bio, Denmark) on the basis of 97% sequence similarity. Taxonomy assignment of non-chimeric reads was determined via the SILVA database (ver. 123) with a 80%confidence threshold. The above processes were also performed using the CLC Genomics Workbench. Validated reads (11,000 per sample) were used to compare the α-diversity (observed OTUs, Chao1, Shannon and Good’s coverage) using QIIME (Quantitative Insights Into Microbial Ecology). Principal coordinate analysis based on the UniFrac distance was implemented to illustrate seasonal β-diversity. Differences in bacterial composition were evaluated using permutational multivariate analysis of variance (PERMANOVA). The network was made using the Cytoscape software (ver. 3.4.0), and was visualized using prefuse force directed layout with correlation coefficients between core bacteria [30]. Since the relative quantity of core bacteria is nonparametric, the Spearman correlation was calculated using the SAS ver. 9.4 software; Spearman’s correlation coefficient (|r| > 0.6) with FDR-corrected significance level under 0.01 was used.

Statistical Analysis

The Student’s t-test or the Mann-Whitney test was conducted after Shapiro-Wilk normality to analyze all data aside from network co-occurrence. Spearman’s correlation analysis was performed to construct the network co-occurrence of core genera.

Comparison of Diversity Indices and Bacterial Population among Samples

A total of 2,235,280 reads (average 73,325 reads for April samples and 75,760 reads for July samples) were analyzed after quality checking the 30 lettuce samples. The read numbers were normalized to 11,000 by random subsampling (Table 2). The average number of observed OTUs in the April samples (1,450 ± 628.8) were determined to be higher (p = 0.0005) than those of the July samples (788.2 ± 187.6). The number of observed OTUs at site A in April (2,169 ± 575.9) was the highest, whereas samples at site B in April (586.7 ± 141.8) had the lowest number of OTUs. Unlike other sites, only OTUs obtained from lettuce samples collected at site B were higher (p = 0.1695) in July (785.3 ± 148.8) than those collected in April (586.7 ± 141.8). The microbial diversity of the lettuce microbiota was evaluated by Chao1 richness (Table 2). The average estimate of Chao1 richness from samples in April (1,955 ± 677.3) was higher (p = 0.0012) than that in July (1,275 ± 267.6). The highest Chao1 richness value was obtained from site A samples in April (2,730 ± 502.9), whereas the lowest was obtained from site B samples in April (945.4 ± 219.8). The diversities obtained by the Shannon index (Table 2) of April samples (4.9 ± 1.95) were significantly higher than those of July samples (3.63 ± 0.89, p = 0.0292). Site A samples in April (6.95 ± 1.67) were the most diverse group, whereas site E samples in July (2.64 ± 1.26) demonstrated the lowest diversity index.

Table 2 . Summary of diversity indices obtained from sequencing..

Sampling time	Sample	Analyzed reads	Normalized reads	Observed OTUs	Estimated OTUs (Chao1)	Shannon diversity index	Good’s coverage
April	SiteA_1	44,866	11,000	2505	3019.58	8.30	0.98
	SiteA_2	44,194	11,000	1504	2148.95	5.08	0.98
	SiteA_3	59,017	11,000	2498	3020.41	7.47	0.99
	SiteB_1	89,185	11,000	423	693.43	2.14	0.99
	SiteB_2	39,668	11,000	663	1098.12	3.91	0.99
	StieB_3	13,263	11,000	674	1044.66	5.54	0.97
	SiteC_1	102,411	11,000	1468	1933.63	4.61	0.99
	SiteC_2	77,424	11,000	1373	1837.33	4.33	0.99
	SiteC_3	148,143	11,000	1017	1816.48	1.07	0.99
	SiteD_1	70,046	11,000	2212	2703.39	7.22	0.99
	SiteD_2	96,852	11,000	1726	2390.30	4.65	0.99
	SiteD_3	88,016	11,000	1606	2123.98	4.05	0.99
	SiteE_1	81,446	11,000	1598	2076.95	5.45	0.99
	SiteE_2	44,212	11,000	1367	1791.90	6.28	0.99
	SiteE_3	101,153	11,000	1123	1625.32	3.38	0.99
July	SiteA_1	91,723	11,000	895	1472.77	4.11	0.99
	SiteA_2	59,029	11,000	629	961.65	3.87	0.99
	SiteA_3	62,771	11,000	798	1181.63	4.22	0.99
	SiteB_1	95,582	11,000	947	1526.57	4.74	0.99
	SiteB_2	76,072	11,000	654	1262.91	3.54	0.99
	StieB_3	76,487	11,000	755	1277.45	4.20	0.99
	SiteC_1	58,530	11,000	844	1261.41	3.39	0.99
	SiteC_2	75,037	11,000	702	1103.71	3.09	0.99
	SiteC_3	58,080	11,000	600	1030.50	2.58	0.99
	SiteD_1	110,471	11,000	1147	1613.84	4.56	0.99
	SiteD_2	95,677	11,000	949	1471.60	4.19	0.99
	SiteD_3	90,779	11,000	879	1646.01	4.03	0.99
	SiteE_1	65,870	11,000	589	1034.60	1.99	0.99
	SiteE_2	58,241	11,000	994	1539.77	4.09	0.99
	SiteE_3	62,056	11,000	441	741.30	1.83	0.99

The bacterial amounts were quantified by 16S rRNA gene analysis (Fig. 2). Results showed that total bacterial amounts in the July samples were significantly higher than those in April (p < 0.0001). On average, the total bacterial cell numbers were 2.26 × 10⁷ and 4.66 × 10⁸ in April and July, respectively. Samples collected from site D showed the greatest variance in total bacterial population between April (1.87 × 10⁷) and July (1.13 × 10⁹). Site E samples exhibited the least difference in total bacterial population between the two time points (1.34 × 10⁷ in April and 6.29 × 10⁸ in July). In the April group, the diversity estimators from site A samples were the highest, but total bacterial counts were the lowest. In contrast, site D samples showed the highest diversity estimators in the July group, which were accompanied by the highest bacterial load. Between April (13.8°C) and July (25.5°C), the average temperature showed an increase by more than 10°C (Fig. 1), which may have an effect on bacterial loads.

Figure 2. Comparison of total bacterial loads in the lettuce at each time and sampling site. The amounts of total bacteria were determined using quantitative real-time PCR. April and July are represented as blue and red, respectively. The average bacterial loads at each time (A) and sampling site (B, C) are indicated The asterisks (***) indicate significant difference between the April and July groups (p < 0.0001). Error bars indicate the standard error of the mean.

In previous studies, it has been reported that a low bacterial load in lettuce was correlated with high species diversity in US lettuce [31]. Therefore, such trend also seems to be resulted from crop features as well as seasonal temperatures.

These results showed that the microbiota of lettuce from South Korea vary with season and region. The locations of sites B, C, and D were adjacent to each other, whereas site A and site E were located farther away. Despite the close proximity of sites B, C, and D, the species diversity of site B differed from the other two sites. This might be due to field conditions such as UV light, temperature, humidity, and water, which all have an effect on the phyllosphere microbiota [32]. Results indicated that the quantity of a particular genus is influenced by the time and site of collection, which leads to variations in the composition of the microbiota. Since a single taxon can vary with the geographical location and time of year, regional differences also seem to play a role in the microbial diversity of lettuce [33].

Microbiota Composition Differences between April and July Samples

Weighted and unweighted PCoA plots were generated on the basis of the UniFrac distance to compare the differences in bacterial composition between lettuce collected in April and July (Fig. 3). Clearly distinguished clustering patterns (PERMANOVA, p = 0.00001) were observed on the unweighted matrix according to sampling time of these samples. On the weighted PCoA plots, similar separated patterns based on sampling time were observed, which was confirmed by PERMANOVA (p =0.00001). These results may explain the effect of season on the microbiota composition. In the clustering tree analysis based on weighted UniFrac for lettuce collected in the United States, differences in clustering due to seasonal variations were attributed to environmental features, such as air and soil temperature, solar radiation, and humidity [31].

Figure 3. Weighted/unweighted UniFrac-based principal coordinate analysis (PCoA). (A) Weighted and (B) unweighted. The April and July groups are represented as blue and red, respectively. The percentage contributions to the variance of the data from principal components 1, 2, and 3 (PCo 1, PCo 2, and PCo 3) are listed along axes representing them.

Phylum, Classs and Family Levels

The OTUs from 2,236,280 non-chimeric reads were analyzed at each taxonomy level to compare the microbiota of samples collected in April and July. Firmicutes and Proteobacteria were the dominant phyla in the April and July samples, respectively (Fig. 4A). Proteobacteria (average relative abundance of 49.20 ± 17.38%) were more abundant (p = 0.4021) than Firmicutes (39.27 ± 18.11%) in the April samples. However, Firmicutes (52 ± 19.67%) were slightly more abundant (p = 0.6293) than Proteobacteria (45.80% ± 19.40) in the July samples. Actinobacteria contributed to 10.7 ± 6.57% in the April samples and 2 ± 0.28% in the July samples. Proteobacteria, Firmicutes, and Actinobacteria have previously been reported as the predominant phyla in the plant phyllosphere [34].

Figure 4. Comparison of the taxonomic composition at the phylum, class, and family levels. (A) Comparison of bacterial composition at the phylum level. The legend only shows the phyla found at more than 0.5% in at least one sample. (B) The eight dominant classes (defined as a >1% based on total read abundance in at least one sample) are represented. Others indicate the remaining phyla and classes. (C) Comparison of dominant families (>5% of the microbiota in at least one sample) between each site in April and July.

At the class level (Fig. 4B), the most dominant classes in the collected samples were Bacilli (43.10 ± 18.21%) within Firmicutes, and Gammaproteobacteria (42.27 ± 16.94%) within Proteobacteria. Actinobacteria were expressed (p = 0.0019) in greater proportion in the April samples (10.37 ± 5.58%) than in the July samples (1.98 ± 0.66%). The Clostridia class, which belongs to Firmicutes, was distinctly observed at site B in July. Bacilli are highly utilized in agricultural biotechnology, and many of their products are currently on the market [35]. Therefore, the many Bacilli species found in the samples may have been due to their use as fertilizers. It has been shown that various genera, including pathogens, belong to Gammaproteobacteria, and can be found in plants such as peanuts, flowers, and poplar [11, 36, 37]. Lastly, since lettuce grows close to the soil, it was speculated that anaerobic Clostridia found in the samples may have been derived from soils buried in the lettuce [38].

The microbiota of the lettuce at the family level (Fig. 4C) were determined within 16 dominant families (defined as a >5% of total read abundance in at least one sample). Among them, Bacillaceae, Bacillales family XII group (incertae sedis), Enterobacteriaceae, and Pseudomonadaceae were the dominant families in lettuce collected in both April (average relative abundance 17.58%, 14.06%, 8.59%, and 6.77%, respectively) and July (32.46%, 27.51%, 6.45%, and 5%). Studies on microbial communities of lettuce roots showed that Pseudomonadaceae are core bacteria, and are correlated with Sphingomonadaceae [39]. However, the abundance of Sphingomonadaceae was found to be very low in our study samples (average relative abundance under 1%), except for samples collected at site A in April.

Micrococcaceae, Oxalobacteraceae, and Planococcaceae were found in samples collected in April (6.42%, 4.73%, and 2.29%), rather than (p < 0.0001, p = 0.0026 and 0.0107, respectively) July (0.88%, 0.5%, and 0.69%). Micrococcaceae members are a major component of bacterial air spores, and are highly expressed in the phyllosphere microbiota [40]. Oxalobacteraceae can be found in various environments, and some species are known as mild pathogens or opportunistic human pathogens [41]. Planococcaceae members have been reported to exist in extreme environments such as deep sea sediments, marine solar salterns, glaciers, permafrost, Antarctic deserts, and sea ice brine [42].

According to a report on bacterial families of the microbiota in lettuce cultivated in the USA, the relative abundance of Enterobacteriaceae in lettuce collected during the summer was higher than that in the winter owing to increased solar radiation and humidity [34]. On the contrary, Oxalobacteraceae accounted for a large proportion of microbial species in the same crop collected during the winter as opposed to summer. These findings were similar to our results of Korean lettuce collected in April and July.

Genus Level

Heatmap analysis was used for prevalent genera (over 1% detected in at least one sample) to compare the microbiota of lettuce collected in April and July (Fig. 5A). The most dominant genera (average relative abundance of >1%) observed in the April samples included Arthrobacter (5.94%), Bacillus (8.06%), Exiguobacterium (14%), Pantoea (12.95%), Pseudomonas (6.74%), and Xanthomonas (9.73%). In lettuce collected in July, Acinetobacter (5.50%), Bacillus (6.18%), Enterobacter (17.72%), Enterococcus (11.23%), Exiguobacterium (27.53%), and Pseudomonas (5%) were identified as the major genera. Bacillus, Exiguobacterium, and Pseudomonas were the dominant genera observed in lettuce sampled at both time points.

Figure 5. Heatmap analysis of the frequency of detected genera and the potential pathogens detected in lettuce in sampling sites in April and July. (A) The clustering was performed using Spearman`s rank correlation coefficient. Different colors indicate the relative abundance of each genus. Sample names indicate the sampling site and time (e.g., 4A; sample collected from site A in April). (B) The relative abundance of potential pathogens at each time is represented. The potential pathogens were determined by comparison with sequences in the PATRIC database (http:// www.patricbrc.org). (C) The presence of potential pathogens in lettuce was determined by PCR and electrophoresis. The PCR products were separated on 1.2% agarose gels except E. casseliflavus (2%). Lane M : DNA size marker (ExcelBand).

Exiguobacterium occupied a high percentage of samples collected in both April and July (14% and 27.53%, respectively). It is highly resistant to temperature, pH, and salt [43, 44]. Therefore, it can exist anywhere in nature, and has been reported to be a root microorganism that promotes plant growth [11]. Pseudomonas (April with 6.74%and July with 5%) discovered in the collected lettuce was inferred to exist as an opportunistic phytopathogen, similar to species found in the phyllosphere of grasses [45]. Bacillus (April with 8.06% and July with 6.18%) belonging to Bacilli is distributed in large areas; it exerts antifungal activity and acts as a biological control agent. Species have been used in plant health management [46, 47]; however, some reports have shown that some plants also carry the human pathogenic B. cereus [48]. Arthrobacter was found in lettuce collected in April, but not in July (p = 0.0036). Arthrobacter has been reported to be capable of producing pigments that contribute to cell membrane stability at low temperatures, and can survive in cold environments [49]. Enterobacter was found to be more abundant (p < 0.0001) in July (17.72%) than in April (0.28%). Enterobacter sp. is a biological control agent for many fungal phytopathogens [50]. It promotes plant growth, and is used for agriculture [51]. Because this bacterial species prefers humid environments [52], it was not surprising that more bacteria of this species were identified in July. Acinetobacter, known as a cause of nosocomial infection, was also found to be more abundant (p = 0.0044) in July than in April [53], because it prefers humid environments.

Observation of Potential Pathogen Species

The sequences obtained in lettuce samples at different sampling times and location were compared with those of reported pathogenic bacteria in the PATRIC (http://www.patribrc.org) database. Several pathogenic bacteria, including Bacillus spp., Enterococcus casseliflavus, Klebsiella pneumonia, and Pseudomonas aeruginosa, were found (Fig. 5B). Bacillus spp. included B. anthracis, B. cereus, B. mycoides, B. pseudomycoides, B. thuringiensis, and B. weihenstephanensis. This group including B. cereus, one of the most frequently reported foodborne pathogens [54], was observed in both the April (0.08%) and July (0.05%) samples. On the other hand, various kinds of pathogenic bacteria were detected more in July lettuce than in April lettuce. E. casseliflavus, an opportunistic premise plumbing pathogen, was more frequently (p = 0.0024) observed in lettuce collected in July (0.29%, April 0.02%) [55]. K. pneumoniae, which usually causes nosocomial infections in hospitals [56], was found only in July lettuce (0.1%). This suggests that a feces-based fertilizer was used during cultivation, because Klebsiella bacteria are normally found in the human intestines and feces. There were several reports showing that food products with feces containing enteroinvasive K. pneumoniae may be a potential source of life-threatening foodborne illness [57]. The proportion of P. aeruginosa that can cause a disease in humans as well as plants was also higher in July lettuce (0.0628%) than in April lettuce (0.0005%), although not statistically significant (p = 0.0937).

However, the identification of taxonomies on the species level using 16S rRNA is restricted. Considering the results of genus characterization, therefore, the PCR was performed using primers for detection of common bacterial pathogens (Fig. 5C). Eight primer sets for pathogen detection were used (5 sets for Bacillus spp. detection; 1 set each for E. casseliflavus, K. pneumoniae, and P. aeruginosa detection; respectively). In the case of Bacillus spp., it was determined only when all five primers set amplified positive bands. Although not found in April, Bacillus sp. was detected in one sample of the B site in July lettuce. E. casseliflavus was detected at both sampling times. It was detected at A, C, D, and E sites in April and at B and D sites in July. K. pneumoniae and P. aeruginosa were detected only in July, with K. pneumoniae at the B and D sites, and P. aeruginosa at the A, B, C, D, and E sites. These results showed similar tends with the pathogen analysis using the bacterial bioinformatics database.

These results demonstrate that the potential risks of food poisoning through lettuce consumption might be higher in July. Although the proportion of pathogenic bacteria in lettuce was relatively low, these bacteria could be added to the screening list of pathogens on lettuce. As we expected, the potential risks of these samples seemed to be assessed indirectly when looking at pathogens present in the microbiota.

Network Co-Occurrence Revealed the Correlation between Core Genera

Microbial communities are affected not only by the environmental factors, but also by the relationship between microbiomes. A previous report has shown that closely related species, also known as microbial hubs, form microbial networks that play an important role in plant health [17]. With the present study, microorganisms in the microbial community of foliar may not all be “residents,” but may be “transient” [58]. A method for analyzing interactions between “residents” and “transients” has recently been developed [39]. The core bacterial groups (present in more than 80% of samples with more than 0.1%in average relative abundance) in lettuce sampled in April and July were analyzed on the basis of correlations between genera. Unlike what we predicted by considering only season-dependent temperature changes, the bacterial network formed by different genera in lettuce collected in April was more complicated than that collected in July.

The bacterial network of the April group was composed of 33 nodes (genera) and 125 edges (Fig. 6A). Among various genera, Brevibacterium, Brevundimonas, Deinococcus, Enterococcus, and Peptoclostridium were found to correlate with a large number (over 10 edges) of core genera. The highest positive correlation was observed between Enterococcus and Lactococcus (r = 0.96416, p < 0.0001) among the networks formed by core genera in the April samples. However, only Xanthomonas, Bacillus, and Arthrobacter belonging to the predominant genera of April lettuce were part of this network. Bacillus and Xanthomonas correlated with a small number (1 and 3 edges) of other genera, but Arthrobacter correlated with a relatively large number (6 edges) of genera. This result is consistent with a previous report showing the network of lettuce root microbiota [39]. In July samples, the network was composed of 32 nodes and 53 edges (Fig. 6B); Acinetobacter, Enterococcus, Paracoccus, and Tatumela correlated with the largest number (5 edges) of core genera in the July network. Among these core bacteria, Lysinibacillus showed the strongest positive correlation (r = 0.94995, p < 0.0001) with Geobacillus, whereas Enterococcus (r = -0.87009, p < 0.0001) showed the strongest negative correlation with Acinetobacter. The six genera (Acinetobacter, Bacillus, Enterobacter, Enterococcus, Exiguobacterium, and Pseudomonas) dominant in the July samples were included the network. The dominant genera were involved in a large number of networks, and were associated with many core bacteria, similar to that previously observed in flower microbiota [12]. Genera with fewer connections in the network can be thought of as single organisms that do not require complex connections to invade and colonize [39]. Among potential pathogenic genera, the positions of Bacillus and Staphylococcus were different in the April network. Staphylococcus was involved in the main network with four edges, whereas Bacillus was separated from the main network with three genera. In the July network, only Acinetobacter demonstrated multiple connections (5 edges), whereas Escherichia-Shigella, Pantoea, and Pseudomonas had a single connection. These results may help explain the correlations between resident and transient microorganisms in lettuce cultivated in South Korea. In addition, it is likely that foodborne pathogens recognized well in lettuce belong to the transient type. In Fig. 7, Brevibacterium, the resident microorganism, showed strong negative correlation with both Acinetobacter (r = -0.8206, p = 0.0002) and Pseudomonas (r = -0.7077, p = 0.0032), which belong to potential pathogenic genera. This can be used as a basis for biomarker research to find harmful bacteria. These results suggest that the interactions between resident microorganisms and potential pathogens in the network can be used as basic information to prevent food poisoning.

Figure 6. Network of co-occurring genera of the lettuce in April and July. The network of co-occurring genera was constructed in Cytoscape using correlation data of lettuce in the (A) April and (B) July groups. The nodes (bacterial genera) are colored according to modularity class. A concentration represents a strong (Spearman’s correlation coefficient |r|>0.6) and significant (p < 0.01) correlation. The size of each node is weighted average relative abundance of genus. The edges colors represent the type of correlation (green: positive, red: negative) and the edges length and width were weighted normalize value.

Figure 7. Correlation of relative abundance between Brevibacterium and the potential pathogenic genera (A) Acinetobacter and (B) Pseudomonas.

16S rRNA gene-based analysis of the microbiota in lettuce collected in South Korea showed different microbiome compositions not only between April and July, but also between sites. A number of potential pathogens, including well-known foodborne pathogens, were detected. The network co-occurrence of core bacteria was more complex than expected in April, as compared with that in July. These correlations suggested that new pathogens may emerge on the premise of ecology changes and microbe-microbe interactions [59]. A more detailed understanding of the microbiota in lettuce should be developed for the production of safe lettuce and prevention of food poisoning due to raw consumption of lettuce. Further studies are needed to better understand the interactions between the core genera and their potential pathogenicity for safe intake of lettuce.

This work was supported by a grant from the National Research Foundation of Korea funded by the Korean Government (NRF-2014R1A1A1002580) and the Ministry of Food and Drug Safety, Republic of Korea (14162MFDS972).

The authors have no financial conflicts of interest to declare.