Next-generation sequencing identies seasonal and geographical differences in the ruminal community composition of sika deer (Cervus nippon) in Japan

To understand the nutritional conditions of culled wild sika deer (Cervus nippon) and their suitability as a pet food component, we compared the ruminal community compositions of deer living in different food habitats (Nagano winter, Nagano spring, and Hokkaido winter) using next-generation sequencing. Twenty-nine sika deer were sampled. Alpha and beta diversity metrics determined via 16S and 18S ssrRNA amplicon-seq analysis showed compositional differences. Prevotella, Entodinium, and Piromyces were the dominant genera of bacteria, fungi and protozoa, respectively. Moreover, 66 bacterial taxa, 44 eukaryotic taxa, and 46 chloroplastic taxa were shown to differ signicantly among the groups by linear discriminant analysis effect size (LEfSe). Total RNA-seq analysis showed 397 signicantly differentially expressed transcripts (q < 0.05), of which 258 were correlated with bacterial amplicon-seq results (Pearson correlation coecient > 0.7). The amplicon-seq results indicated that deciduous broadleaf trees and SAR were enriched in Nagano, whereas graminoids, Firmicutes and fungi were enriched in Hokkaido. The ruminal microbial community were corresponded with different food habits, related to the severe snow conditions in Hokkaido in winter and the richness of plants with leaves and acorns in Nagano winter and spring. These ndings are useful for understanding the nutritional conditions of wild sika deer. In the present study, using 16S and 18S ssrRNA amplicon-seq and total RNA-seq analysis, we revealed the difference in the population-level microbiome and transcriptome diversity between seasons (spring/winter) and locations (Nagano/Hokkaido) in wild sika deer in Japan. The diversity was shown by alpha diversity metrics in Euk and beta diversity metrics in Bac and Euk. The results of Lefse showed the enrichment of specic Bac and Euk taxa, and a difference in food habits was shown by the enriched Chl taxa. Given the previous reports and results in the present study, clues about the relationship between the ruminal community (bacteria, fungi, and protozoa) and food habits (limited vs. abundant feeding conditions) were provided. The high correlation between the results of amplicon-seq and total RNA-seq analysis showed the accuracy of the present data and ability to provide information about live cell function (locomotion ability and glycolysis) of the ruminal bacteria. These data are useful not only for understanding the nutritional conditions in wild sika deer in Japan but also for understanding microbial community dynamics during the fermentation process in ruminants.


Introduction
In Japan, the extinction of apex predator wolves, global warming, reduced snowfall, and the declining popularity of hunting have increased the number of sika deer (Cervus nippon) by approximately 10-fold from 1990 to 2014, with the population estimated to be 3.05 million animals [1]. Because deer feed on plants for survival, an excess of deer disturbs forest ecosystems, damages farm crops, and increases the prevalence of tick-borne diseases in humans in cities [1,2]. The Japanese government has tried to control their population by culling and using deer resources for some purposes, such as "game meat" for human consumption [1]. However, the high costs of disposing of carcasses are an economic burden on local authorities [3]. To solve this problem, the use of deer resources as pet food has recently been attempted to cover facilities costs [3]. and in association with diet in livestock and wild deer [13]. A few amplicon-seq studies have examined the effect on rumen bacterial, fungal, and protozoal composition in wild deer. In one study. Li et al. [17] reported that tannin-rich plants altered the rumen microbiome and fermentation patterns in captive sika deer in China. Wilson et al. [13] reported an association between diet and rumen microbiota in wild roe deer, but protozoa were observed in only 1 % of the population and were unevenly distributed. These studies examined the factors that drive changes in rumen microbial diversity. Thus, understanding ruminal community compositions is important for nutritional conditions and growth.
To promote the use of culled sika deer for human and companion animal consumption in Japan, we were interested in examining their ruminal community diversity to understand host nutritional conditions. To understand the quality of wild deer meat, it is necessary to collect more information about wild deer biology and ecology. In general, the amount and type of feed eaten by wild deer depend on the location and season. Feed intake increases in spring and reaches a peak in summer [21]. Stanisz et al. reported that deer meat obtained in winter showed a greater brightness and a reduced redness compared to meat obtained in summer [22]. Thus, seasonal variation in food habits is likely an important determinant of deer meat quality. Nonetheless, little information about the relationship between food habits and meat quality is available to date.
Many wild sika deer inhabit Mount Asama, Nagano, the central northern area of Honshu Island (36'N, 139'E), and Shiretoko, the northernmost island of Hokkaido (44'N, 145'E), Japan. In Shiretoko, the winter is long and cold (below zero) with heavy snowfall, so sika deer experience limited access to food [10,23].
In Nagano, as reported in Ashio, Tochigi (36'N, 139'E) [24], and Akagi, Gunma (36'N, 139'E) [25], oak forests and conifer stands are dominant, and sika deer can obtain much more food in winter in Nagano than in Shiretoko. Similar to the Poland study [22], it is possible that geographical and seasonal variations in food habits affect meat quality in Japanese sika deer.
Total RNA sequencing (RNA-seq) aims to describe microbial ecology and the transcriptome simultaneously. It has been widely used in studies of marine ecosystems [26,27], soil microbes [28,29] and rumen microbiomes [30,31] to obtain information from all domains of microbial habitats, including eukaryotes, archaea, and bacteria. In particular, unlike the genomic DNA sequencing approach, this RNAbased transcriptome approach can describe differences in functional gene expression patterns among live samples.
In this study, we examined the effects of season (spring/winter) and location (Nagano/Hokkaido) on the population-level microbiome and transcriptome diversity in wild sika deer using 16S and 18S ssrRNA amplicon-seq and total RNA-seq analysis.
The mean read count per sample, OTUs per sample, alpha diversity metric and related statistics for season and location groups (NW, NS, and HW) of the ruminal community (Bac, Euk, and Chl) are summarized in Table 1. For all Bac, Euk, and Chl amplicon-seq, rarefaction curves were constructed to ensure su cient sequencing depth (89,000, 24,000, and 9,000 reads, respectively) to evaluate the dominant members of the ruminal community ( Fig. 1). In Bac, signi cant differences were not observed among the groups; however, in Euk, signi cant differences were observed among the groups in terms of Shannon index (p < 0.01) and Pielou's evenness (p < 0.01). In particular, NW showed the lowest eukaryotic richness (q < 0.05) and evenness (q < 0.02). In Chl, signi cant differences among the groups were observed only in the observed OTUs (p < 0.05) and Shannon index (p < 0.05). (Table 1).
Taxon-level differences in ruminal community composition using amplicon-seq Linear discriminant analysis (LDA) effect size (LEfSe) was used to compare the taxon abundances among sampling seasons and locations within the species level (level 7). LEfSe identi ed 66 Bac taxa, 44 Euk taxa, and 46 Chl taxa with LDA scores above 3 and signi cant differences (p < 0.05) among NS, NW, and HW (Tables S4 to S7). Relevant features on these taxonomic or phylogenetic trees were visualized by cladograms ( Fig. 4 A, B, and C).
Of the 66 Bac taxa, 36 were included in HW, whereas 15 were included in both NS and NW. p_Firmicutes was enriched in HW, occupying in almost half of the circle, including the 30 enriched taxa in the cladogram (Fig. 4A, Table S4). In p_Firmicutes, c_Clostridia and o_Clostridiales were enriched in HW. In clade o_Clostridiales, of the 13 taxa, 5 genera of f_Lachnospiraceae (including g_Lachnoclostridium 10, g_Lachnospiraceae XPB1014 group, and g_Lachnospiraceae AC2044 group) and 5 genera of f_Ruminococcaceae (including g_Ruminococcaceae UCG011, g_Ruminococcus 1, and R. avefaciens) were enriched in HW. In the same clade of o_Clostridiales, 9 taxa were enriched in NW, including identi able f_Ruminococcaceae, g_Ruminococcus 2, g_Ruminococcaceae UCG_014, and s_Ruminococcus 1 bacterium V9D2012. NS enriched only the g_Lachnospiraceae NK3A20 group in o_Clostridiales. In p_Firmicutes, g_Kandleria and K. vitulina were enriched in NS. The enriched species of g_Succiniclasticum were observed in HW and NS.
The 3 single rows of the clades (total 14 taxa) within c_Saccharimonadia, c_Gammaproteobacteria, and c_Oligo exales were enriched in HW, NS, and NW, respectively. Of the remaining 10 taxa, 4 descendent species (UWB11, UWCM, F. succinogenes, and uncultured) in p_Fibrobacteres were enriched in HW. Three taxa within o_Gastranaerophilales and 2 species within g_Treponema 2 were enriched in NS. NW was enriched in c_Mollicutes RF39NW, identi ed alone in the clade.
Differential expression transcriptome analysis of ruminal community composition using total RNA-seq analysis A total of 2,998 predicted ORF sequences (predicted transcripts) were obtained from a total of 4,356,106 read counts (mean 161,337 ± 41,251 SD read count per sample). The normalization and DE analysis for the predicted transcripts were performed using TCC-baySeq. Of the 2,998 predicted transcripts, 397 statistically signi cant DE transcripts (q < 0.05) were obtained by clearing 56 nonDE transcripts. Of the 397 DE transcripts, the numbers in each class (HW > others, NS > others, NW > others, others > HW, others > NS, and others > NW) were 193, 70, 18, 53, 38, and 25, respectively. The results of the homology search using BLAST+ against SwissProt showed the annotation of 119 DE transcripts among the 397 DE transcripts. The summary of predicted taxa for the DE transcripts annotated by Swiss-Prot in each class order is described in Table 3, and the raw numerical data for the DE analysis are provided in Table S7.
Correlation between the results of 16S ssrRNA amplicon-seq and total RNA-seq analysis To examine the association between the results of amplicon-seq and total RNA-seq, the differential taxa (DE taxa) were identi ed by TCC baySeq at the species level (level 7) of the 16S ssrRNA amplicon-seq result. Then, the normalized counts of the annotated DE transcripts were compared with those of the DE taxa for analysis of covariance. A visual summary of the correlation is provided in Fig. 5, and the raw numerical data used for the analysis of covariance are provided in Table S8

Discussion
The dominant ruminal species in Japanese sika deer Our study examined the ruminal community composition of wild sika deer using next-generation sequencing and subsequent analysis. The dominant genus of Bac taxa ( Fig. 3A and Table S1) revealed by amplicon-seq in the present study were consistent with the ruminal composition reported previously in sika deer [16,17,19] and related species (wild roe deer [13], elk [14,18], reindeer [20], ruminants [6]). Of the previous reports, Henderson G et al. [6] reported the dominant core rumen bacteria in ruminants from 32 species and a total of 742 individual animals. The core bacteria were similar across all ruminants; the top mean relative abundances of bacterial groups in 60 deer were reported here to be almost 70 % and were composed of Prevotella, unclassi ed Clostridiales, Bacteroidales, Ruminococcaceae, Lachnospiraceae, Veillonellaceae, and Ruminococcus. Our results indicated that Japanese sika deer have similar components of dominant bacterial species as other ruminants, especially wild ruminants.
In the analysis of the dominant Euk taxa ( Fig. 3B and Table S2), our results showed that g_Piromyces and g_Entodinium were assigned as dominant genera. These genera have been reported as the dominant rumen fungi and protozoa in ruminants [6], sika deer in Japan [10], China [17], and red deer in New Zealand [32].
Thus, we considered the results of our study to be informative to understand the diversity of ruminal community compositions associated with different food habits developed over a long history of adaptation to natural habitats in Japanese sika deer. In addition, it could be useful for understanding the rumen composition of livestock ruminants.
Genera enriched in the rumen among the groups The compositions were different among the groups, and the order of most to least diversity was NW, NS, and HW, as shown by alpha diversity metrics in Euk and Chl and beta diversity metrics in Bac, Euk, and Chl (Figs. 1 and 2, and Tables 1 and 2).
Of the Bac taxa, LefSe revealed that g_Prevotella 1 and Succinivibrionaceae were enriched in NS ( Fig. 4A and Table S4). Prevotella ruminicola, P. albensis, P. brevis, and P. bryantii have been reported to have expanded ruminal niches under starch-, pectin-, xylan-, and sugar-rich conditions [5]. It has been reported that Prevotella and Succinivibrionaceae are more abundant in animals fed diets containing concentrate than in those fed forage-rich diets [6]. These taxa are possible major producers of propionate and the propionate precursor succinate and so are responsible for the greater levels of propionate formed under concentrate-rich diets [6]. It has been reported that the proportion of Prevotella and the ruminal concentration of ammonia decreased on day 28 after the start of cellulose-rich corn silage feeding in sika deer. A positive association between Prevotella and the ruminal concentration of butyrate was also shown previously [16]. On the other hand, in the present study, cellulolytic bacteria such as Lachnospiraceae and Ruminococcaceae [33] were enriched in HW (Fig. 4A and Table S4). In particular, R.
avefaciens and F. succinogenes, enriched in HW, are well-studied cultivable bacteria 5, 34,35]. Consistent with our results, it has been reported that the sika deer rumen ora contain a higher proportion of R. avefaciens in winter than in summer at Shiretoko [11]. Under winter conditions in the same habitat, a shortage of food corresponds to a decrease in ruminal concentrations of total VFAs, propionate, butyrate, ammonia and mineral contents in sika deer [10,36]. Given the accordance between the results of previous reports and the present study, it was supposed that the amount of nutrients in the rumen, such as cellulose, ammonia and total VFA, could be associated with the proportion of the Bac taxa enriched in Nagano and Hokkaido, such as Prevotella and cellulolytic bacteria, respectively ( Fig. 4A and Table S4).
In Euk taxa, the relative abundance of Entodinium in HW was remarkably low in comparison to NW and NS ( Fig. 4B and Table S5). Our results are consistent with the fact that Entodinium density in winter decreases to 1/10 of that in summer at Shiretoko [10]. In winter, due to the limited availability of deciduous broadleaf tree leaves, sika deer eat foods such as ground plants, bark, twigs, and lichens, which are less abundant and have poorer crude protein than tree leaves [10]. The sika deer in HW therefore could be considered adapted to utilizing poor-quality diets, because rumen protozoa have been reported to have minor contributions to ber digestion [10].
The whole transcriptomes of Entodinium caudatum have been reported using RNA-Seq analysis [37]. E. caudatum is considered to prefer starch since E. caudatum has more transcripts involved in the use of starch than cellulose and hemicellulose. In addition, E. caudatum ferments sugars to VFA. In the present study, C. glandium was enriched in NW ( Fig. 4B and Table S5). In winter in Nagano, sika deer might prefer to eat acorns because weevil larvae often grow in acorns over the winter [38]. It is known that acorns themselves are rich in starch and lipids [39]. The ruminal concentration of total VFAs was reported to be higher in summer and autumn than in winter [10]. It is known that rumen protozoa affect the meat quality of ruminant mammals in terms of the fatty acid composition of the meat. It has been reported that rumen protozoan membranes contain high concentrations of conjugated linoleic acids (CLAs) and vaccenic acid [40], and the total ciliate number in the rumen has a positive correlation with cis-9, trans-11 CLAs and vaccenic acid depositions in lamb meat [41]. These ndings suggested that the rumen Entodinium could play an important role in storing nutrients in the muscle for severe winter seasons after consuming fallen leaves and acorns in autumn.
In Euk taxa, the anaerobic fungi of f_Neocallimastigaceae and g_Piromyces, well-established species with utility for lignocellulose bioprocessing [42], were more enriched in HW than in NW and NS (Fig. 4B and Table S5). It has been reported that the transcripts of biomass-degrading enzymes are repressively controlled by bacteria via glucose stimulation [42]. The comparison of transcriptomes among anaerobic and aerobic fungi and rumen and nonrumen bacteria showed that P. rhizin ata had almost twice as many putative pectinases as the other anaerobic lignocellulolytic fungi [43]. These putative pectate lyases were conserved among the anaerobic fungi, whereas these enzymes were much less commonly found in the other groups of organisms, particularly bacteria [43]. Under winter conditions at Hokkaido, nutritional shortages could be required for Piromyces to degrade pectin within the bark. In our results, putative glycolysis and glycogenesis estimated DE transcripts in rumen bacteria observed in HW (Table 3 and Fig. 5). It was considered that the number and function of rumen lignocellulolytic fungi might be regulated by cellulolytic bacteria to overcome severe winter conditions.

Seasonal differences in the food habits of wild sika deer in Japan
The alpha and beta diversity metrics of Chl taxa were different among NW, NS, and HW, suggesting seasonal changes in the food habits of sika deer in Japan (Figs. 1 and 2, and Tables 1 and 2). The present study indicates that the NW diet is enriched in broadleaf trees (Fagales, Morus, and Ilex), which bear many acorns and winter buds (Fig. 4C, and Table S6), whereas the NS diet is enriched in spring-to summer-owering deciduous broadleaf trees (Populus, R. pseudoacacia, and A. quinate), evergreen trees (C. camphora, Fokienia) and forbs (P. tenella) (Fig. 4C, and Table S6). The food habits of sika deer living close to Nagano has been described previously [24,25]. In winter to spring in the wildlife reserve area of Ashio, Tochigi, the predominant rumen contents of sika deer were culm sheaths, dead leaves, barks and twigs, conifers, and graminoids (from most to least prevalent) [24]. In autumn to winter at mount Akagi, Gunma, the diet was composed of graminoids, dead leaves of broadleaf trees, nuts, and berries [25]. These results indicated that the sika deer in Nagano use mostly broadleaf tree leaves and nuts during winter to spring.
In the present study, the rumen composition of HW was enriched in graminoids (Poales) and forbs (Panax) (Fig. 4C, and Table S6). In winter, at Shiretoko, Hokkaido, bark and twigs, graminoids, and leaves were the main foods in the rumen contents of sika deer [10]. In February, at Ashoro and Onbetsu, Hokkaido, the rumen contents were composed of the broadleaf tree leaves, bark and twigs, graminoids and forbs [23]. Among the Poales, Sasa bamboo (Sasa nipponica) has been reported to be a predominant food species, especially in winter, for wild sika deer in Japan [10, 23,24,25]. S. nipponica has evergreen leaves, and its nutritional content is moderate compared with that of other plants [23]. The Chl taxa enriched in HW indicated that Poales is one of the important nutrients for sika deer in Hokkaido during severe winter periods with heavy snow.
The associations among nutritional source, fermentation, and rumen community composition Given the discussion above, it is hypothesized that rumen bacteria (e.g., Lachnospiraceae, Ruminococcaceae) and anaerobic fungi (P. nnis) enriched in HW prefer to grow in cellulose-rich and starch-and lipid-poor conditions, while abundant starch and lipids are fermented by Prevotella, Succinivibrionaceae and ciliates (e.g., Entodinium) to generate large energy stores. It is suggested that the rumen bacteria and anaerobic fungi in HW extract glucose and xylose from the fermentation of cellulose, hemicellulose and lignin while saving and using small energy sources. It has been reported that the number of rumen bacteria increases and that of Entodinium decreases in winter based on morphological examinations of deer [10,44]. It is known that signi cant numbers of ruminal bacteria can be consumed by protozoa, resulting in an inverse relationship between protozoal and bacterial densities [5]. Total RNA-seq analysis showed that many glycolysis-and glycogenesis-estimated DE transcripts in the rumen were observed in the HW > other group (Table 3). Furthermore, the majority of the DE taxa, highly associated with the predicted DE transcript were included in the taxa enriched in HW (Fig. 5). Our results suggested that the bacterial activity was enhanced under rumen conditions in HW. The decrease in protozoa indicated that the interaction of rumen bacteria with anaerobic fungi compensates for the presence of protozoa.

Limitations
We could not examine fermentation parameters in the rumen. The Nagano and Hokkaido sika deer were different subspecies (C. n. centralis and yessoensis, respectively), and an effect of host genetic background on the rumen microbiome and metabolites in sika deer and elk hybrids has been reported [15]. Further research is warranted to examine the physiological status and quality of game meat from sika deer in Japan.

Correlation between total RNA-seq and amplicon-seq analysis
High correlations between the annotated DE transcripts and DE bacterial taxa were observed in the present study (Fig. 5, and Table S8). The correlated DE taxa were consistent with the genera indicated as enriched by LefSe, such as f_Lachnospiraceae (E. ruminantium), g_Lachnoclostridium 10, g_Lachnospiraceae XPB1014 group, g_Fibrobacter, and g_Ruminococcus 1 in HW and g_Prevotella 1 in NS (Table S4 and Table S8). The correlated agellin genes could indicate that rumen bacteria have motility and can adhere to plant tissues [45]. The correlation with the EMP pathway, LDH, PTS, and sucrose catabolism could indicate that rumen bacteria utilize glucose, sucrose, and pyruvic acid for glycolysis resulting in the production of VFAs [46,47]. Biological information converted from annotated transcripts has been used to study the seasonal prevalence of functional marine microorganisms [27]. The high quantitative expression gene data from the total RNA-seq analysis could indicate the accuracy of genomic data from amplicon-seq from different perspectives. Although the Ribo-tag approach for ssrRNA data is now being developed [28], total RNA-seq may be superior to amplicon-seq when information is required from a live subject. Meanwhile, it has been indicated that both approaches show varying results [30,48], so the combination of different approaches is desirable to permit an evaluation of the accuracy of ssrRNA data.

Conclusion
In the present study, using 16S and 18S ssrRNA amplicon-seq and total RNA-seq analysis, we revealed the difference in the population-level microbiome and transcriptome diversity between seasons (spring/winter) and locations (Nagano/Hokkaido) in wild sika deer in Japan. The diversity was shown by alpha diversity metrics in Euk and beta diversity metrics in Bac and Euk. The results of Lefse showed the enrichment of speci c Bac and Euk taxa, and a difference in food habits was shown by the enriched Chl taxa. Given the previous reports and results in the present study, clues about the relationship between the ruminal community (bacteria, fungi, and protozoa) and food habits (limited vs. abundant feeding conditions) were provided. The high correlation between the results of amplicon-seq and total RNA-seq analysis showed the accuracy of the present data and ability to provide information about live cell function (locomotion ability and glycolysis) of the ruminal bacteria. These data are useful not only for understanding the nutritional conditions in wild sika deer in Japan but also for understanding microbial community dynamics during the fermentation process in ruminants.

Animals and Sampling
Twenty-nine, 16 male and 13 female rumen samples were obtained from adult sika deer, 24 living in mount Asama, Nagano and 5 living in Shiretoko, Hokkaido. Twelve samples from Nagano were collected in the winter season between December 2019 and January 2020 (referred to as NW), and the other 12 samples were collected in the spring season during April 2020 (referred to as NS). All 5 Hokkaido samples were collected in the winter season during February 2019 (referred to as HW). The mean temperature (maximum/minimum) at Nagano was 7.0/-4.0°C in December-January and 16.0/-2.0°C in April, and that at Hokkaido was − 2.2/-9.6°C in February-March. The monthly snow depths were 16 cm in April 2019, 27 cm in December 2019, 26 cm in January 2020 at Karuizawa, Nagano (36'N, 138'E) and 123 cm in February 2019 at Utoro, Hokkaido (44'N, 145'E), which were close to the sampling area. The hunted sika deer were immediately dressed in the eld at Hokkaido, and the culled deer were moved to a slaughter facility operating in Komoro city, Nagano. One milliliter of rumen juice was carefully put in and suspended in a sterile 15 ml conical tube lled with 9 ml RNAlater solution (Thermo Fisher Scienti c Japan, Tokyo, Japan) by a clean plastic spoon. Then, the tubes were sent to the laboratory at Azabu University at 4°C within 3 days. The tubes were kept in a -80°C freezer until DNA and RNA preparation.

Animal ethics
Collecting rumen samples from sika deer was approved by the ethics committee of animal experiments at Azabu University (approval 150917-1) and conducted in accordance with the ethical standards of The Mammal Society of Japan (http://www.mammalogy.jp/en/guideline.pdf).

DNA and RNA preparation
One hundred microliters of the RNAlater-xed rumen content suspension was collected by a pipet chip. The suspension was centrifuged (21130 x g, 1 min), and then the supernatant was discarded. The remaining pellet was used for DNA and RNA extraction steps with a DNeasy PowerSoil kit (QIAGEN, Venlo, Netherlands) and RNeasy PowerBio lm kit (QIAGEN), respectively. The extracted DNA and RNA were eluted and dissolved in 100 µl of water. DNA and RNA concentrations were measured by a Qubit DNA and RNA HS kit (Thermo Fisher, Waltham, MA, USA).

Construction and quanti cation of DNA and RNA libraries
The prepared DNA was used for the construction of 16S and 18S ssrRNA gene sequence libraries using the "16S Metagenomic Sequencing Library Preparation" protocol (Illumina, San Diego, CA, USA) with some modi cations. Brie y, PCR enzyme "KOD plus" (TOYOBO, Osaka, Japan) was used with a standard Read Archive under the accession numbers DRA011262 for 16S ssrRNA gene sequence library, DRA011263 for 18S ssrRNA gene sequence library, and DRA011264 for Total RNA-seq library.

Sequence analysis
The obtained 16S and 18S ssrRNA gene sequence reads were analyzed using the Quantitative Insights Into Microbial Ecology (QIIME2) pipeline [51] with DADA2 [52] for quality control. The generated operational taxonomic units (OTUs) were assigned using a naïve Bayes classi er for annotation [53]. Each target region was speci ed by primer sequences to strains in the naïve classi er with silva132_99.fna of the Silva database, release 132 [54]. Annotation was performed on taxonomy_7_levels.txt in the same database. When the identi ed plant species were not distributed in Japan, genus-level or family-level taxa were represented in the manuscript (raw data are shown in Tables   S6 and S9).
Mapping was performed by bowtie2 v.2.3.5.1-linux-x86_64 [60] with option -U and -local. The resulting SAM les were transformed with SAMtools into bam les and sorted. Sorted BAM les were used to obtain counting information by eXpress v.1.5.1-linux_x86_64 [61]. Count data were truncated with a custom script to remove reads with fewer than 10 counts. The expressed genes were annotated for the expected ORF sequence by TransDecoder software. The created expected ORF sequences were annotated by BLAST + against the Swiss-Prot database.

Statistical analysis
In amplicon-seq analysis, based on the taxonomy generated, we excluded reads originating from archaea, chloroplasts, and host mitochondria. We removed singleton features and elements with fewer than 10 reads. The ltered feature tables and computed relative abundance per taxonomic level were used for microbial community diversity analysis. Rarefaction curve analysis of the obtained data was performed to estimate the completeness of microbial community sampling. Subsequently, alpha diversity (Shannon's diversity index, observed OTUs, Faith's phylogenetic diversity, and Pielou's evenness) was computed for differences in the number of species and species richness per sample), and beta diversity metrics [Jaccard distance, Bray-Curtis distance, unweighted UniFrac distance, weighted UniFrac distance and generated principal coordinates analysis (PCoA)] were plotted using Emperor to estimate microbial community similarities between groups. Group signi cance between the alpha and beta diversity indices was calculated with the QIIME2 pipeline using the Kruskal-Wallis (KW) test and permutational multivariate analysis of variance (PERMANOVA), respectively [51]. The Benjamini-Hochberg false discovery rate (FDR) correction was then used to identify each index differently represented between groups. P-and q-values < 0.05 were considered statistically signi cant. Linear discriminant analysis (LDA) effect size (LEfSe) [62] was performed to identify the microbial taxa that were signi cantly associated with the groups. The LEfSe algorithm consisted of the KW rank sum test for the difference in classes and for LDA in the relevant features. The parameters were set at p = 0.05 and LDA score = 3.0 for the computation.
For the total RNA-seq analysis, all analyses were performed using R functions [63]. Tag Count Comparison (TCC) baySeq [64] was used for normalization and differential expression (DE) analysis on the multigroup RNA-seq count data. The TCC package was generated from original TbT methods (TMM-baySeq-TMM pipeline), consisting of a combination of the trimmed mean of M values (TMM) normalization [65] in edgeR [66] and DE gene detection in baySeq [67]. In this strategy, normalization of count data and DEG estimation are iterated to avoid false positives; we repeated this cycle three times in this study [68]. DE transcripts between groups were de ned as transcripts at FDR (q-value) lower than 0.05.
The Pearson correlation coe cient between amplicon-seq data and transcriptome data was calculated using the "cor" function of R [63]. The obtained matrix was used to draw a heatmap chart with the "gplots" package of R. In this heatmap chart, we visualized correlation factors between 0.7 and 1.0.   Cladograms visualized by LefSe analysis. Red, green, and blue indicate signi cantly different groups (Hokkaido winter, Nagano spring, and Nagano winter, respectively), with the diagram and species classi cation at the phylum, class, order, family, and genus levels shown from the inside to the outside.
Yellow circles represent species with no signi cant difference. (A) Bacterial taxa in 16S ssrRNA, (B) eukaryotic taxa without chloroplasts and host mitochondria in 18S ssrRNA, (C) chloroplast taxa in 18S ssrRNA. Alphabets in the cladograms indicate the signi cantly different taxa of the clade, including more than 2 taxa.

Figure 5
Pearson correlation coe cients of rumen microbial communities between amplicon-seq and total RNAseq. The taxonomy_7_levels.txt in the SILVA database was used for annotation of bacterial OTUs in the amplicon-seq data. BLAST+ against SwissProt was used for the annotation of expected ORF sequences in the total RNA-seq data. The correlation factor is set between 0.7 to 1.0. The color bar indicates a statistically signi cant order among the groups (q < 0.05) on the differential taxa (DE taxon) analyzed by TCC baySeq from amplicon-seq.in

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. KawaraiTable20210805.xlsx KawaraiSupTable20210805.xlsx