The CTLs investigated in this study were Dexue CTLs cultivated in Shifang, Sichuan Province, China. Shifang is known as the cigar production base in China, with a history of nearly 400 years of cigar tobacco planting and cigar production. The agricultural processed CTLs, including freshly harvested, air-cured, and agricultural fermented CTLs were collected, 500–1,000 g each sample, then CTLs were transferred into sterile bags, and stored at −20°C until further study. All CTLs were performed in triplicate. Freshly harvested, air-cured, and agriculturally fermented Dexue Nos.1, 3, 4, and 7 CTLs were marked accordingly (D1X, D3X, D4X, D7X, D1C, D3C, D4C, D7C, D1F, D3F, D4F, and D7F).

The CTLs (5.0 g) was suspended in 100 mL of sterile PBS and shaken for 2 h at 200 rpm, after which the supernatant was centrifuged at 10,000× g for 30 min. Total genomic DNA of each sample was extracted using an EZNA® Soil DNA Kit (Omega, USA) according to the manufacturer’s instructions. The genomic DNA was amplified and sequenced using the V4–V5 hypervariable region of 16S rRNA genes (forward primer 515F: 5’-GTGCCAGCMGCCGCGGTAA-3′; and reverse primer 907R: 5’-CCGTCAATTCMTTTRAGT TT-3′) (Parada et al., 2016) and the ITS1 hypervariable region of internal transcribed spacer (ITS) genes (forward primer ITS1F: 5’-CTT GGT CAT TTA GAG GAA GTA A-3′; and reverse primer ITS2R: 5’-GCT GCG TTC TTC ATC GAT GC-3′) (Usyk et al., 2017).

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).

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.

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.

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.

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% 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).

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 on the left side of the plot. The VIP values were Septoria, Filobasidium, Alternaria, Sampaiozyma, Penicillium, Mycosphaerella, and Symmetrospora in the first and second principal components (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 derivative biosynthesis were the top five functions significantly positively correlated with Filobasidium.