Analysis of the structure and metabolic function of microbial community in cigar tobacco leaves in agricultural processing stage

The agricultural fermentation processing of cigar tobacco leaves (CTLs), including air-curing and agricultural fermentation, carried out by tobacco farmers has rarely been studied. In this study, we have investigated the microbial community in the CTLs during air-curing and agricultural fermentation by 16S rRNA and ITS gene high-throughput sequencing. The results showed that the richness of microbial communities gradually increased with the development of agricultural fermentation, which means that not all microorganisms in CTLs come from the fields where tobacco grows, but gradually accumulate into CTLs during the fermentation process. Enterobacteriaceae, Chloroplast, and Alternaria were the dominant genera in the air-cured CTLs. Aquabacterium, unclassified Burkholderiaceae, Caulobacter, Brevundimonas, and Aspergillus were the dominant genera in the agriculturally fermented CTLs. Acinetobacter, Methylobacterium, Sampaiozyma, and Plectosphaerella first significantly increased, and then significantly decreased during agricultural processing. The changes in microbial communities are mainly related to their different functions during fermentation. This means that when the fermentation effect of the original microbial community in cigar tobacco leaves is not ideal, we can optimize or design the microbial community based on the fermentation function that the microbial community needs to achieve. These results may help adjust and optimize the agricultural fermentation process of CTLs, and help develop the quality of CTLs and increase the income of tobacco farmers.


Introduction
Tobacco (Nicotiana tabacum L.) is the most widely cultivated non-food crop in the world (Banožić et al., 2020). Cigar is a tobacco product that has been dried and fermented. Fresh cigar tobacco leaves (CTLs) after harvest usually require air-curing and agricultural fermentation performed by farmers before they can be sold to industrial companies. Air-curing is done to remove moisture from CTLs and change the color of CTLs from green to yellow. The most important thing is that the macromolecular substances in CTLs (starch, cellulose, proteins, etc.) begin to degrade into reducing sugars and amino acids. Agricultural fermentation can 10.3389/fmicb.2023.1230547 Frontiers in Microbiology 02 frontiersin.org be considered as a natural continuation of air-curing, where macromolecular substances continue to degrade and reducing sugars and amino acids begin to be converted into various aromatic substances. Air-curing and agricultural fermentation are important agricultural processes that determine the quality of cigar leaves (Jin, 1982;Toharisman et al., 2008). However, relying solely on experience in agricultural fermentation can easily cause uneven quality of tobacco leaves from different batches. Microbes have been found to play an important role in the fermentation of CTLs (Reid et al., 1944). Based on cultureindependent molecular biology techniques such as polymerase chain reaction denaturing gradient gel electrophoresis (PCR-DGGE) (Giacomo et al., 2007) and Illumina MiSeq sequencing (Zhang et al., 2018a,b) have been used to characterize microbial communities associated with CTLs, Bacillus, Staphylococcus, Penicillium, Debaryomyces, Jeotgalicoccus, Lactobacillus, Weissella, Yania, Pseudomonas, and Acinetobacter have been identified in many tobacco samples (Li et al., 2009;Du et al., 2016;Zhang et al., 2018a,b;Liu et al., 2021). However, few studies have tracked the succession of structure and function of microbial communities in CTLs during agricultural processing stage.
In this study, we investigated the structure and function of microbial communities in CTLs during the agricultural processing by 16S rRNA and ITS gene high-throughput sequencing. Exploring the dynamic change of microbial community during agricultural processing stage would be helpful to master the change pattern of microbial communities in CTLs, and investigating the function of microbial communities would be helpful to understand the mechanism of microbial community on CTLs. These results may have important contributions to improving the quality of Chinese CTLs.
Deoxyribonucleic acid libraries were validated by TruSeq Nano DNA LT (Illumina, USA) and quantified using a PicoGreen dsDNA Assay Kit (Invitrogen, USA). Amplicons were pooled in equal amounts, and sequencing was performed using a 2 × 300 paired-end (PE) configuration using an Illumina MiSeq sequencing system according to the manufacturer's instructions (Illumina, USA). Operational taxonomic units (OTUs) of qualified sequences were identified using the clustering program VSEARCH version 1.9.6 and the SILVA version 132 database (Rognes et al., 2016) with 97% similarity. The alpha diversity, including the Chao1 and Shannon values, was analyzed using QIIME version 1.9.1 (Caporaso et al., 2010).

Data analysis
Significant differences among the freshly harvested, air-cured, and agricultural fermented CTLs groups were determined using SPSS version 19 (IBM, USA), by one-way analysis of variance (ANOVA) and Duncan's multiple comparison test (p < 0.05). In order to further investigate the changes in microbial abundance, circular visualization diagrams were used to display the overlapping and differentiating microbial taxa among the three courses of agricultural processing of CTLs at the phylum and genus levels using Circos software online (http://circos.ca; Krzywinski et al., 2009). Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) and Fungi Functional Guild (FUNGuild) were used to predict the metagenomic functions based on the normalized OTU tables (Douglas et al., 2019), which were compared using Statistical Analysis of Metagenomic Profiles (STAMP) version 2.1.3 (Parks et al., 2014). Partial least-squares regression (PLS) was conducted using SIMCA-P (version 13.0; UMETRICS, Sweden) to construct statistical models of bacteria and fungi for correlation between microbial genera (X) and functional genes (Y). The data are presented as means, which were scaled before the analysis. In the PLS model, the variable importance for the projection (VIP) reflects the importance of the variables both to explain X and to correlate to Y. Variable importance for the projection values exceeding 1 indicated the 'important' X-variables. The PLS regression coefficients (CoeffCS) are coefficients used for interpreting how strongly Y is correlated with the systematic part of each of the X-variables. Networks were explored and visualized using the interactive platform Gephi (Bastian et al., 2009) based on the CoeffCS of the VIPs, with each node representing one genus and one functional gene and edges representing a strong and significant correlation.

Microbial community diversity
In total, 1,741,054 and 2,095,371 reads of the 16S rRNA and ITS gene sequences were obtained from all CTLs. Table 1 shows the diversity indices of the bacterial and fungal communities in the CTLs. The coverage was more than 99%, indicating that Illumina MiSeq sequencing was deep enough to represent all the detected microbial communities. The Chao1 values of the bacterial communities of the air-cured and agriculturally fermented CTLs were higher than those of the freshly harvested CTLs (p < 0.05). The Shannon values of the microbes in the air-cured CTLs were higher than those in the freshly harvested and agriculturally fermented CTLs (p < 0.05). These results showed that the diversity of microbial communities gradually decreased with the progress of agricultural fermentation, which means that not all microorganisms in CTLs come from the fields where tobacco grows, but gradually accumulate into CTLs during the fermentation process.

Microbial community composition and structure
Taxonomic analysis of the reads revealed that Proteobacteria, Cyanobacteria, Firmicutes, Bacteroidetes, Actinobacteria, Ascomycota and Basidiomycota were dominant at the phylum level ( Figures 1A,B). A total of 610 bacterial genera and 473 fungal genera were detected, of which 45 and 17 had a relative abundance higher than 0.5% in at least one sample ( Figures 1C,D). The dominant bacterial genera, with relative abundance higher than 5.0% in at least one sample, were unclassified Enterobacteriaceae, Pseudomonas, Chloroplast, Acinetobacter, Pantoea, Sphingomonas, Staphylococcus, Aquabacterium, unclassified Burkholderiaceae, Methylobacterium, Caulobacter, and Brevundimonas, whose abundance accounts for 61.81-95.81% in all CTLs. The dominant fungal genera were Aspergillus, Alternaria, Sampaiozyma, and Plectosphaerella.
Differences in the dominant phyla and genera in the different courses of agricultural processing of all samples were shown in the Circos diagrams (Figures 2A,B, respectively), which clearly showed that the dominant phyla and genera were more diverse in the air-cured and agriculturally fermented CTLs than in those in the freshly harvested CTLs. At phylum level, Cyanobacteria were mainly present in the freshly harvested CTLs (D1X, D3X, D4X, and D7X), Bacteroidetes and Basidiomycota were mainly present in the air-cured CTLs (D1C, D3C, D4C, and D7C), and Ascomycota was mainly present in the agriculturally fermented CTLs (D1F, D3F, D4F, and D7F) (p < 0.05). At genus level, unclassified Enterobacteriaceae, Chloroplast, Pantoea, and Alternaria were the dominant genera in the air-cured CTLs. The Acinetobacter, Sphingomonas, Methylobacterium, Sampaiozyma, and Plectosphaerella increased first and then decreased during agricultural processing. Pseudomonas, Staphylococcus, Aquabacterium, unclassified Burkholderiaceae, Caulobacter, Brevundimonas, and Aspergillus were the dominant genera in the agriculturally fermented CTLs.

Metabolic function prediction of microbial communities
The metabolic function of microbial communities was predicted using PICRUSt2 and FUNGuild ( Figure 3). Functions of the bacteria and fungi in the CTLs involved fatty acid and lipid biosynthesis, amino acid biosynthesis, fermentation, and carbohydrate biosynthesis. The abundance of 56 bacterial functions and 20 fungal functions were higher than 1,000. Bacteria are mainly involved in amino acid biosynthesis, fatty acid and lipid biosynthesis, nucleoside and nucleotide biosynthesis, amino acid degradation, and carbohydrate degradation at Kyoto Encyclopedia of Genes and Genomes (KEGG) level 2, including fatty acid biosynthesis, proteinogenic amino acid biosynthesis, purine nucleotide biosynthesis, proteinogenic amino acid degradation, and sugar degradation at level 3. Fungi are mainly involved in fatty acid and lipid biosynthesis, fatty acid and lipid degradation, electron transfer, and respiration at KEGG level 2, including purine nucleotide biosynthesis, pyrimidine nucleotide biosynthesis, fatty acid degradation, aerobic respiration II, aerobic respiration I, and aerobic respiration at level 3.
Significant differences in the functional prediction profiles of the CTLs at different agricultural processing among the three courses were observed in KEGG level 3 (Figure 4). For the bacterial communities, a total of 38 out of 56 (41.01% increased, 28.57% Frontiers in Microbiology 04 frontiersin.org decreased) were significantly different between the freshly harvested and air-cured CTLs, and eight out of 56 (8.93% increased, 5.36% decreased) were significantly different between the air-cured and fermented CTLs ( Figures 4A,B). For the fungal communities, a total of 14 out of 20 (35% increased, 35% decreased) were significantly different between the freshly harvested and air-cured CTLs, and nine out of 20 (25% increased, 20% decreased) were significantly different between the air-cured and fermented CTLs ( Figures 4C,D).

Correlation between microbial genera and functional genes
A model based on Partial least-squares regression was used tialto analyze the correlation between microbial genera and functional genes. And two significant principal components of the total variance in the data matrix were extracted in both the bacterial and fungal models. For the bacterial model, the R2X, R2Y, and Q2 were 0.886, 0.871, and 0.809, respectively, which meant that 88.6% variation was due to these two components, with a total of 87.1% dummy Y variable per class, and 80.9% overall cross-validated R2 for these two components. The data indicated that the PLS model was suitable for this study. Samples of the freshly harvested, air-cured, and fermented cigar leaves were clearly separated on the score scatter plot ( Figure 5A, Table 2), with the freshly harvested group located on the right side of the plot, the air-cured group located on the upper left side of the plot, and the fermented group located on the left-middle side of the plot. The VIP values were unclassified Enterobacteriaceae, Pseudomonas, Chloroplast, Acinetobacter, Pantoea, and Sphingomonas in the first and second significant principal components (Table 3). The coefficients refer to the PLS model being rewritten as a regression model. The correlation network between the VIP genera and functional genes was described by significant CoeffCS, as shown in Figure 5B. There were 52 and 51 edges (functions) significantly positively correlated with unclassified Enterobacteriaceae and Chloroplast, and the top five functions of both were proteinogenic amino acid degradation, proteinogenic amino acid biosynthesis, sugar degradation, purine nucleotide biosynthesis, and vitamin biosynthesis. Eight functions were significantly positively correlated with Pantoea, and the top five functions were lipopolysaccharide biosynthesis, quinol and quinone biosynthesis, sugar acid degradation, sugar derivative degradation, and fermentation of pyruvate.
For the fungal model, the R2X, R2Y, and Q2 were 0.680, 0.826, and 0.760, respectively, which meant that 68.0% of the variation was due to these two components, with a total of 82.6% dummy Y variable per class and 76.0% overall cross-validated R2 for these two components. The data indicate that the PLS model was suitable for this study. Samples of the freshly harvested, air-cured, and fermented cigar leaves were clearly separated on the score scatter plot ( Figure 5C, Table 2), and the freshly harvested group was located on the right side of the plot, the air-cured group was located on the middle side of the plot, and the fermented group was located Plot of the phylum and genus levels' relative abundances for the bacterial and fungal communities in the cigar tobacco leaves (CTLs). Panels (A,C) represent the bacterial communities at the phylum and genus levels. Panels (B,D) represent the fungal communities at the phylum and genus levels. D1X to D7X denote the freshly harvested CTLs, D1C to D7C denote the air-cured CTLs, and D1F to D7F denote the agriculturally fermented CTLs.  Table 3). The correlation network is shown in Figure 5D. The tRNA charging and sugar derivative biosynthesis were significantly and positively correlated with Alternaria. Twenty functions were significantly positively correlated with Sampaiozyma and Symmetrospora, and the top five functions of both were fatty acid biosynthesis, phospholipid biosynthesis, sugar derivative biosynthesis, sugar degradation, fermentation of pyruvate, and fermentation to alcohols. There were 16 and 21 functions that were significantly positively correlated with Penicillium and Mycosphaerella, and the top five functions of both were proteinogenic amino acid biosynthesis, 2′-deoxyribonucleotide biosynthesis, pyrimidine nucleotide biosynthesis, purine nucleotide biosynthesis, and sugar derivative biosynthesis. Phospholipid biosynthesis, 2′-deoxyribonucleotide biosynthesis, pyrimidine nucleotide biosynthesis, purine nucleotide biosynthesis, and sugar Differences in the microbial community composition among the three groups (harvested: freshly harvested cigar tobacco leaves [CTLs]; air-cured: air-cured CTLs; and fermented: agricultural fermented CTLs) at the phylum and genus level. Panels (A,B) represent the relative abundances of the phyla and genus shown in circular visualization (Circos). The thickness of each ribbon represents the abundance of each taxon. The absolute tick above the inner segment and relative tick above the outer segment indicate the read and relative abundances of each taxon, respectively.

Frontiers in
Frontiers in Microbiology 06 frontiersin.org Function prediction of the bacteria (A) and fungi (B) in the cigar tobacco leaves (CTLs). D1X to D7X denote the freshly harvested CTLs, D1C to D7C denote the air-cured CTLs, and D1F to D7F denote the agriculturally fermented CTLs.
Frontiers in Microbiology 07 frontiersin.org derivative biosynthesis were the top five functions significantly positively correlated with Filobasidium.

Discussion
Air-curing and agricultural fermentation are important production processes that determine the quality of CTLs and finished products (Zhang et al., 2020a,b). In this study, the structure and diversity of microbial communities in CTLs at different stages of agricultural processing in Shifang, Sichuan were studied by using the Illumina miseq sequencing method based on 16S rRNA and ITS genes. The results indicated that there are significant differences in the microbial communities of CTLs at different stages of agricultural processing, and the differences in the microbial communities during different stages of agricultural processing are even greater than those between different varieties. Our research helps to grasp the trend of microbial community changes in the agricultural fermentation process of CTLs, reveal the essence of CTL fermentation, and is of great significance for producing stable quality CTLs. In addition, understanding the succession and function of microbial community helps to reasonably regulate the microbial community to produce higher quality tobacco leaves.
At the phylum level, Proteobacteria, Cyanobacteria, Firmicutes, Bacteroidetes, Actinobacteria, Ascomycota, and Basidiomycota were abundant in the CTLs, which is similar to the results of previous studies based on flue-cured tobacco leaves (Zhang et al., 2020a,b). Except for Cyanobacteria, the microbial phyla in the cigar products were associated with the results obtained in this study (Ye et al., 2021). The abundance of Cyanobacteria in the freshly harvested CTLs was Differences in the microbial function composition among the three groups (DX: freshly harvested cigar tobacco leaves [CTLs]; DC: air-cured CTLs; DF: agriculturally fermented CTLs). Panels (A-D) represent the Statistical Analysis of Metagenomic Profiles (STAMP) analysis of the relative abundance of the bacterial and fungal functions. p < 0.05 indicates significant differences between two different site types. Error bars represent Welch's t-interval.
Frontiers in Microbiology 08 frontiersin.org higher than those in the air-cured and agricultural fermented. Cyanobacteria have plant-type photosynthetic organs, which are mainly used to fix carbon dioxide through the reduction of pentose phosphate (Stal, 2015). The most dominant bacterial phylum was Proteobacteria, and the most dominant fungal phylum was Ascomycota in the CTLs from the Shifang, Sichuan province, which was consistent with those of the H382 cigar leaves from the Hainan province (Liu et al., 2021). At the genus level, there were 16 dominant genera, of which unclassified Enterobacteriaceae, Pseudomonas, Acinetobacter, Pantoea, Sphingomonas, Staphylococcus, and Methylobacterium were also the dominant bacterial genera in Mexican CTLs (Zhang et al., 2018a,b) and Hainan CTLs . The dominant genera were unclassified Enterobacteriaceae, Chloroplast, Pantoea, and Alternaria in the freshly harvested CTLs. The unclassified Enterobacteriaceae in the CTLs may be from the soil (Akita et al., 2019). Chloroplast are oxygen-producing phototrophic microorganisms that can form symbiotic relationships with plants and usually provide nitrogen fixation to plants (Rai et al., 2000). Alternaria can decompose cellulose in CTLs (Macris, 1984). After air-curing, the dominant microbes were Acinetobacter, Sphingomonas, Methylobacterium, Sampaiozyma, and Plectosphaerella. These microbes were reported to participate in the degradation of nicotine (Wang et al., 2011), lignin (Masai et al., 1999), and one-carbon compounds such as methanol and methylamine (Rossetto et al., 2011). After fermentation, the dominant microbes were Pseudomonas, Staphylococcus, Aquabacterium, unclassified Burkholderiaceae, Caulobacter, Brevundimonas, and Aspergillus. These microbes are mainly involved in the degradation of nicotine (Tang et al., 2009;Li et al., 2010), pectin, protein, starch (Schuster et al., 2002), xyloside, and oligosaccharides (Coenye, 2014), as well as the production of amino acids and fatty acids (Aro et al., 2010;Tabanelli et al., 2012) to form a series of flavoring substances and carbohydrates, thus improving the quality of CTLs. Although Aspergillus could improve the quality of cigar leaves by degrading nicotine or producing organic acids such as citric acid, it could cause mildew in cigar leaves (Yan et al., 2018). These results indicate that the microbial community in CTLs is not entirely derived from the soil where the tobacco plants are grown, but gradually recruited from the external environment during the fermentation process. The shifts in community structure Score scatter three-dimensional plots of the partial least-squares regression (PLS) (A,C) and network of correlation between the microbial genera and functional genes (B,D) based on the PLS. The PLS of the bacterial (A), fungal (C), and various CTLs is represented as a two-dimensional representation of the scores (t[1] and t[2]) on the first and second PLS components. Relationship between the functional communities and bacteria (B) and fungi (D), with significant changes in abundance. For each panel, the size of each node is proportional to the number of connections, nodes of the same color are affiliated with the same genus or function, and the thickness of each connection between two nodes is proportional to the value of the PLS regression coefficients (CoeffCS) with statistical significance (p < 0.05). D1X to D7X denote the freshly harvested CTLs, D1C to D7C denote the aircured CTLs, and D1F to D7F denote the agriculturally fermented CTLs.
Frontiers in Microbiology 09 frontiersin.org and composition mainly follow the changes in CTLs, originating from changes in their metabolic function (Zhou et al., 2020;Chen et al., 2021). Microorganisms are mainly involved in carbohydrate degradation, and biosynthesis of fatty acids, amino acids, and aroma components. Functional genes of fatty acid and lipid biosynthesis, aromatic compound degradation, and amino acid biosynthesis had higher relative abundances in the freshly harvested CTLs than in the air-cured and fermented CTLs, which meant microbes could synthesize fatty acids and amino acids and metabolize aromatic compounds in order to improve the flavor and aroma of CTLs (Wang, 2003). PLS is particularly useful when we need to predict a set of dependent variables from a large set of independent variables (Abdi, 2010). In this study, we used PLS to analyze the correlation between microbial genera and functional genes. PLS indicated that samples could be separated into three groups according to the processing process, which is consistent with the sampling results. Bacteria, mainly including unclassified Enterobacteriaceae, Chloroplast, and Pantoea, biosynthesized quinol and quinone in the freshly harvested CTLs the highest, which might promote the color change in the CTLs. Polyphenolic compounds in the CTLs can form brown pigments such as quinones through an enzymatic browning reaction in the air-curing process to make CTLs turn dark, thus affecting the color and aroma of the tobacco leaves (Guo et al., 2009). These bacteria also contributed to the degradation of total and reducing sugars through sugar degradation, sugar acid degradation, polysaccharide degradation, and sugar derivative degradation, and Alternaria contributed to sugar derivative biosynthesis, resulting in lower contents of total sugar and reducing sugar of the air-cured CTLs compared with those of freshly harvested CTLs (Wang, 2013).
In summary, our study systematically investigated the microbial community in CTLs during agricultural processing stage. The microbial communities vary widely due to fermentation time differences. The richness of microbial communities gradually increased with the development of agricultural fermentation. The changes in microbial communities are mainly related to their different functions during fermentation. This means that when the fermentation effect of the original microbial community in cigar tobacco leaves is not ideal, we can optimize or design the microbial community based on the fermentation function that the microbial community needs to achieve. These results may help adjust and optimize the agricultural fermentation process of CTLs, and help develop the quality of CTLs and increase the income of tobacco farmers.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.

Author contributions
DL and JZ conceived and designed the manuscript. QZ, TZ, ZY, SY, WC, PL, and YH handled samples and conducted experiments. QZ and TZ wrote and revised the manuscript. All authors read and approved the manuscript.

Conflict of interest
QZ, DL, ZY, SY, WC, PL, and YH were employed by China Tobacco Sichuan Industrial Co., Ltd.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.