DNA metabarcoding Uncover the Diet of Subterranean Rodents in China

Objective In the past, the zokor, which lived in northern China, caused great harm to agriculture and forestry production due to its large and sophisticated diet. Since the rat lives underground for most of its life, researchers know little about its dietary habits. Further understanding of its diet in the field is of important meaning for developing green and sustainable control strategies for the rat. Methods Longde County in Liupan Mountain area of Ningxia Hui Autonomous Region was selected as the interest area to capture zokor and investigate the species of habitat plants.We selected chloroplast trnL gene and eukaryotic internal transcription spacer 1 (ITS 1) primers to amplify DNA from the gastric contents of zokor,and then sequenced on Illumina Miseq PE300 platform. Results The gastric contents of Eospalax smithii (n=16) and E.cansus(n=9) were analyzed by operational taxonomic units (OTU) clustering and amplicon sequence variants(ASVs).The OTU clustering method obtained 2,995 OTUs, and the ASV method obtained 4,657 ASVs. The ASV method was better than the OTU clustering method, and the ASV method was adopted in the subsequent analysis. The food list of 32 families, 80 genera and 154 species was obtained by ASV method after the error was removed. The food composition of zokor was evaluated by relative abundance(%RA) method and frequency of occurrence(%ROO) method. At the Family level, it was found that zokor mainly fed on Asteraceae, Poaceae, Rosaceae, Pinaceae, Brassicaceae, Apiaceae, etc. At the Genus level, they are mainlyEchinops, Littledalea,Artemisia,Picea, Cirsium, Elymus and so on. The diet alpha diversity of E.cansus was slightly higher than that of E.smithii (P > 0.05). The correlation coefficient between Sobs index of alpha diversity and body weight of zokor was −0.382 (P = 0.059). The diet beta diversity proved that most zokors (22/25) clustered together, with low heterogeneity. They fed positively on Calamagrostis, Cirsium, Echinops, Medicago, Sanguisorba and Taraxacum. We found that zokor mainly fed on the roots of perennial herbs(PH), which were rich in water, carbohydrate, fat and protein, which provided an important energy source for its survival. Conclusion High-throughput sequencing(HTS) based DNA metabarcoding technology has well revealed the diet of zokor, which are generalist.

Introduction 1 vision degeneration, smell and hearing developed. Underground rats need a lot of energy 4 to dig and maintain a complex system of tunnels, especially in hard, dry soil [3]. This 5 high energy consumption forces the underground rat to expand its food forage and 6 become a food generalist [2,4]. In addition, Heth analyzed the feeding behavior of moles 7 (Spalacidae) and concluded that animals that forage underground, with large search 8 costs and diverse food resources, must collect all edible food species found when digging 9 burrows.This, combined with the non-directed search pattern, will produce a generalist 10 foraging behavior [5]. 11 Myospalacinae is a subterranean rodent endemic to East Asia [6], and mainly dis-12 tributed in a wide area north of the Yangtze River Basin in China. According to 13 our research and literature records, there are currently 9 species of zokor in China, 14 including 6 species of Eospalax and 3 species of M yospalax [6]. The average daily 15 food intake of E.cansus is 105.7g for herbaceous plants and 60.5g for tree roots, and 16 there is a highly prominent positive correlation between daily food intake and body 17 weight [7]. The distribution of the zokor in China also conforms to Bergmann's rule, 18 that is the body weight of the zokor increases with the latitude. The body weight 19 of E.rothschildi living in Qinba Mountains(QBM) was markedly lower than that of 20 E.f ontanierii, E.cansus, and E.smithii living in Loess Plateau(LP). It is also signifi-21 cantly smaller than M.aspalax in the northern Mongolian Plateau(MP). At the same 22 time, the body weight of M.aspalax was remarkably higher than that of the three zokors 23 in the LP.(W QB < W LP < W MP ,P < 0.01; See Fig S1). In QBM, zokor lives mainly 24 in farmland and alpine meadow. In LP, the species and population of zokor were the 25 highest owing to less woodland and more cultivated land. In MP, the woodland is scarce, 26 the range of movement of zokor is large, and the harm to the grassland by feeding on 27 grass roots is larger. Studies have shown that the zokor has caused serious losses to 28 the industry of new woodland [8][9][10], crops and herbs [11,12]. The damage of zokor 29 is particularly serious on LP. Chinese scientists' research on the zokor began with the 30 prevention and control of the disaster of zokor. 31 The gnaw of zokor causes serious losses in new woodland and young forest. In the 32 light of our investigation and literature records, in terracing fields without prevention 33 and control measures, the damage rate of seedlings in new woodland can be as high as 34 82.5% [8]. With the deepening of the concept of ecological development, the control 35 of zokor has shifted from pesticide poisoning and physical killing to the concept of 36 protecting the purpose tree species [13]. The purpose of this study was to explore the 37 diet of zokor, to understand the interaction mechanism between zokor and plant, to 38 explore the status and function of zokor in the ecosystem, and to objectively analyze 39 the degree of harm of zokor, so as to provide scientific basis for making effective rodent 40 control measures [14]. 41 At present, traditional methods are used to study the diet of zokors. For instance, 42 captive feeding [16], observation of food accumulation in caves [21], observation of 43 foraging behavior and microhistology of gastric contents [22,23], etc. The above methods 44 are limited in terms of season and food diversity, and can only reflect the dietary trait 45 of zokor in part of the time period, it is even more difficult to reveal the diet in wild 46 zokor. Microhistology is the highest resolution technique in the current study of diet 47 of zokor. This method requires accurate analysis of partially digested plant fragments 48 but may overemphasize the indigestible groups [24]. However, it is still limited by the 49 researchers' ability of plant classification and heavy workload, and these factors make 50 the study on eating habit of zokor not comprehensive and objective. 51 In recent years, with the improvement of molecular biology database and the emer-52 gence of high-throughput sequencing (HTS). DNA metabarcoding, which integrates the 53 idea of DNA barcode and the HTS technology, is born at the historic moment. DNA 54 metabarcoding is not restricted by prey species, and can identify prey at the taxonomic 55 Species level, especially suitable for diet research of difficult-to observe or special habits, 56 and allows parallel processing and sequencing of large samples [25]. 57 In the selection of diet samples, researchers generally choose different sampling 58 methods in line with the study species. The most common type is feces, because it is 59 the easiest way to obtain and has the minimal damage to predators, especially some rare 60 and mysterious species [20,24,26,27].However, due to the digestion and decomposition of 61 animals, the food DNA in feces degrades too much. Therefore, regardless of the species, 62 the best sample is the stomach contents, where the degradation degree of food in the 63 stomach is the lowest and relatively more food DNA sequences can be obtained [28]. DNA 64 metabarcoding technology mainly uses universal primers to amplify food DNA fragments 65 of animals, so the selection of primers is particularly crucial. For herbivores, chloroplast 66 gene fragments such as trnL(UAA) intron and its above P6 Loop [4,18,24,29,30], 67 rbcL [31][32][33] and trnH-psbA [18,46] are usually selected. Amplified the fecal DNA of 68 herbivorous geese by trnH-psbA, rbcL, matK and trnL respectively, resaechers found 69 that the success rate of trnL amplification was the highest [18]. Another study suggests 70 that the combination of three primers can improve the food variety by 30% [26]. The use 71 of trnL (UAA) P6 loop combined with ITS1-F/ITS1P oa-R and ITS1-F/ITS1Ast-R can 72 improve the species resolution of Poaceae and Asteraceae [24]. Therefore, it is extremely 73 vital to combine different types of primers for the study of the diet of zokor, which is a 74 broad-feeding animal.

75
For a long time, the food sequence analysis of metabarcoding was mainly based on 76 the OTU clustering method. The idea was to generate operational taxon Unit (OTU) by 77 clustering the sequencing results after amplification and sequencing of molecular markers 78 (barcodes). Each OTU includes sequences within the sample that differ very little 79 (usually within 3%). There is no specific classification information of OTU and it cannot 80 be matched one-to-one with species. Simulation studies have found that the number 81 of OTU obtained by OTU clustering analysis is often far more than the actual number 82 of species, which indicates that OTU clustering analysis has a high number of false 83 positives and is too dependent on reference database to be compared between different 84 samples [34].Recently, researchers have developed new methods to resolve amplicons 85 sequence variants(ASVs) from Illumina-scale amplicons data without imposing arbitrary 86 distinct thresholds that define molecular OTUs [34]. Many microbiology-related studies 87 have compared OTU clustering method with ASV method. The ASV method maintains 88 the extensive ecological pattern observed by OTUs in the analysis of fungal diversity, 89 including diversity sequencing, and at the same time improves the reproduceability and 90 data sharing of fungal studies, which can reproduce the biogeographic information hidden 91 by OTU clustering method [35].Under the condition that the sequencing depth is sufficient 92 to capture the complexity of the community, the ASV method is superior to the OTU 93 clustering method in terms of community diversity, especially in terms of fungal-related 94 sequences [36]. ASV method has the strongest correlation to infer pollutants and sample 95 biomass, and the ASV method has good sensitivity and precision, which excel traditional 96 methods [37]. The data set of QIIME2 analyzed by ASV could identify the sample 97 collection point more accurately after removing the genera with lower relative abundance 98 (< 1%) [38]. In terms of intestinal microbiota of aquaculture shrimp: Compared with 99 traditional OTU clustering methods, ASV method has practicality in finding a wider 100 range of vibrio species diversity [39]. A number of recent environmental DNA studies 101 have also chosen to use ASV methods in data processing. Such as visualization of seasonal 102 variations in airborne pollen [40], food choices of endangered small mammals [41], and 103 food composition of native and alien invasive fish [42]. It is a necessary attempt to use 104 ASV method to process and analyze the data of diet sequencing.
woodland and young forest is extraordinarily serious. The research on local zokor is 108 mainly focused on prevention and control, and the basic research on the diet of zokor 109 is still in the traditional direct observation method [8][9][10]. It is helpful to provide 110 theoretical basis for the control of zokor in this region and even in China to study the 111 diet of zokor by using the cutting-edge DNA metabarcoding technology.

121
Study area flora 122 We investigated the flora in the study area on September 3, 2020, October 26, and April 123 18, 2021. The plant species were identified according to the morphological character-124 istics of F lora of China and F lora of N ingxia. References to the PPG I system for 125 ferns [42], Christenhusz [43] for gymnosperms, and APG IV for angiosperms [44]. Due 126 to the limited survey time, the common plants in the same area were also referred to in 127 F lora of Liupan M ountain and A Guide to Common P lants on the Loess P lateau. 128 Seven plots with an area of 1ha were set within the main activity area of zokor. The 129 species, quantity and height of arborous layer, shrub layer and herbal layer were investi-130 gated in the plots with an area of 1ha.

131
Acquisition of stomach contents of zokor 132 As zokor is a local pest, our samples came from a "rat hunt" by the local forestry 133 department. We obtained the wild zokor with the approval and assist of Forest Pest 134 Control station of Natural Resources Bureau of Longde County, Ningxia Hui Autonomous 135 Region, and the research area is not a nature reserve. The rat is trapped in a ground 136 arrow trap, which causes the rat to become unconscious within a short time of being hit, 137 greatly reducing its pain. The harvesting of zokor and the collection of stomach contents 138 was approved by the Forest Pest Control Station in Longde County. All experiments were 139 approved and supervised by the Science Ethics Committee of Northwest A & F University. 140 Dissection was performed after recording relevant body indexes of zokor. The digestive 141 tract was stripped out, and the stomach contents of zokor were clamped with sterilized 142 scissors and forceps into a 2mL cryopreservation tube. After collection, the contents were 143 immediately stored in liquid nitrogen gas phase for subsequent experiments. Among the 144 captured zokors, 25 were selected for diet analysis. Among them, E.smithii(n= 16; 8 145 male ,8 female), E.cansus(n = 9; 6 male,3 female).( Table 1.) The average body weight 146 of male was 315.0g (±40.7 SD, n =14) and the female zokor was 204.2g (± 21.2SD, n = 147 11). The body weight of the two species of zokor in the study area collected by us was 148 analyzed (n=44). There was no significant difference in body weight between different 149 species of zokor, and the body weight of male zokor was significantly higher than that of 150 female (P <0.01) (data not shown).

162
Chloroplast trnL-c-F trnL-h-R and eukaryotes ITS 1 primers were selected for PCR 163 amplification. Transgen AP221-02:T ransstart FastPFU DNA Polymerase,20µL reaction 164 system: 5 × T ransStart FastPfu buffer(4µL),2.5 mM dNTPs(2µL), forward primer(5µM) 165 0.8µL, reverse primer(5µM) 0.8µL, T ransStart FastPfu polymerase 0.4µL,BSA 0.2µL, 166 template DNA 10 ng, up ddH 2 O to 20µL, 3 replicates for each sample.   [45]was used for quality control of the original sequence, and FLASH 179 (v1.2.7) [46] was used for splicing. Filter the reads with a tail mass value of 20 or less 180 bases and set a window of 50bp. If the average mass value in the threshold is lower than 181 20, remove the rear bases from the threshold. Filter the reads with a tail mass value of 182 50bp and remove the reads containing N bases. According to the overlap relationship 183 between PE reads, pairs of reads were merged into a sequence, and the minimum overlap 184 length was 10bp. The maximum error ratio allowed in overlap area of spliced sequences 185 was 0.2, which screened the non-conforming sequences. The samples were differentiated 186 according to the barcode and primers at both ends of the sequence, and the sequence 187 direction was adjusted. The number of mismatches allowed by barcode was 0, and the 188 maximum number of primer mismatches was 2.

189
We used OTU clustering and ASV methods to compare the processed data and 190 compare the differences in species annotation between the two methods. OTU clustering: 191 bioinformation statistical analysis was performed on OTU at 97% similar level, and 192 UPARSE7.0.1090 was used for OTU clustering. NT V20200604 database was selected, 193 and RDP Classifier2.11 was used for sequence classification annotation. Usearch7 was 194 used for OTU statistics, Mothur1.30.2 was used for Alpha diversity analysis, Qiime1.9.1 195 was used to generate abundance tables for each taxonomic level, and the distance of 196 beta diversity was calculated. ASV methods: Based on the default parameters, the 197 DADA2 [48]plug-in in QIIME2 process [47] was used to de-noise the optimized sequence 198 after quality control splicing. The representative sequences of ASV were classified 199 and identified using the "NR/NT Collection" (https://blast) of GenBank database 200 NCBI. https://blast.ncbi.nlm.nih.gov/Blast.cgi. Taxonomic annotation of ASV in the 201 "Nucleotide Collection (NR/NT)" of GenBank database NCBI was performed using the 202 multi blast method in QIIME2 (V2020.2). Sequence consistency: 0.8 Sequence coverage: 203 0.8. Through the diversity cloud analysis platform (QIIME2 process) of Majorbio Bio-204 Pharm Technology Co. Ltd. for subsequent data analysis. Two indexes, %RA (Relative 205 abundance) and %FOO (Frequency of Occurrence) [20,27], were used to evaluate the 206 food composition of zokors. % RA refers to the percentage of the occurrence frequency of 207 the sequence number (ASVs) of a certain food species in the total occurrence frequency 208 of all food ASVs. The calculation formula is as follows: N i is the sum of the ASV occurrence times of all the food of the animal. %FOO 210 refers to the percentage of the samples containing ASV of a certain food in the total 211 number of samples, and the calculation formula is: represents the number of samples with ASV of class i food, and N is the number of 213 effective samples [27]. Plant species in the study area 216 We surveyed the flora in the study area in September 2020, October 2020 and April 217 2021, respectively. In the study area, after the policy of returning farmland to forest 218 and grass was implemented, the main tree species planted were P icea crassif olia, 219 P inus sylvestris, Larix gmelini, Betula platyphylla, Armeniaca sibirica and P opulus 220 davidiana. HTS showed that effective food DNA was amplified from all gastric content samples. 225 We used the traditional OTU clustering method to obtain 2,694,529 sequences. After 226 denoising, 1,986,045 sequences were obtained, and 2,995 OTUs were obtained by cluster-227 ing. 99.80% OTU annotated as Eukaryota. In Eukaryota's sequence, Fungi accounted 228 for 67.01%, Viridiplantae 31.32%, Metazoa 1.62%, and unclassified 0.05%. We also used 229 the ASV method to obtain 3,514,769 sequences, and 1,852,644 valid sequences were 230 obtained by DADA2 denoising and 4,657 ASVs. 99.81% of the effective ASV sequences 231 were annotated to Eukaryota, and Fungi accounted for 69.55%, Viridiplantae 29.52%, 232 Metazoa 0.82% and unclassified Eukaryota's sequences 0.11%.

233
Studies have shown that zokor is a pure herbivorous animal, with a well-developed 234 cecum for decomposing cellulose in plants [50]. Meanwhile, we classified the fungi groups 235 obtained, and found that they were mainly pathogenic fungi of plants, endophytic fungi, 236 fungi in soil, and fungi intrinsic to the stomach of zokor, but not large edible fungi. 237 Metazoan analysis found that the main species were ticks and mites on plants, as well as 238 Canis. The above species may be absorbed by the zokor when feeding on the roots and 239 stems of plants, and whether there is active feeding needs further study. In the analysis, 240 we deleted all the above groups and did not do the analysis, which mainly focused on 241 Viridiplantae.

242
OTU vs ASV taxonomic analysis 243 We used the OTU clustering method and the ASV method to make taxonomic annotation 244 on the amplified feeding sequences, and compared the types annotated by the two methods 245 at the Order, Family, Genus and Species levels. At the level of Order, Family and Genus, 246 there is a high degree of coincidence between the annotated information obtained by 247 the two methods, but at the level of Species, the number of overlapped species of the 248 two methods decreases, and the annotated information of the two methods at the level 249 of Species is greatly different. We found that ASV method can improve the diversity 250 of food species than OTU clustering method, that is, more species can be annotated. 251 Therefore, ASV method was selected for subsequent dietary analysis. ( From the Pan species curve, we can see that the curves of species at Order and Family 254 level have been flat, but at Genus and Species level, the curves are still in an increasing 255 stage, indicating that our sample size cannot fully reflect the food composition of 256 the study area, and further expansion of the sample size is needed in the follow-up 257 studies. (Fig.1a) According to the Core species curve, with the increase of sample size, 258 at 4-5 samples, there are Core species at Order and Family level, while the number of 259 Core species at Genus and Species level tends to zero. (Fig.1b In the 7 study areas, there were 4 shared food families (Asteraceae, Poaceae, Rosaceae 262 and Convolvulaceae) (Fig.2a,3a) , and 4 shared genera (Artemisia,Echinops, Cirsium 263 andConvolvulus) (Fig.2b,3b). The main plant genera wereArtemisia, Echinops and 264 Cirsium(3b). In addition, there were significant differences in food species of zokor 265 in different areas. Combining with plant species investigation in the study area, the 266 differences in food species of zokor in different study areas were consistent with the 267 differences in plant species among plots, which indicated that zokor was a typical 268 opportunist who first fed on plants near the burrow [21]. Among the food composition of the two species of zokor, at the Family level, 17 plant 270 families were feeding for both species of zokor (Fig.3c), and at the Genus level, 31 plant 271 genera were feeding for both species ( (Fig.3d). The Venn diagram indicated that the 272 dietary overlap of the two species was high, especially at the Family level.

Diet of zokor after the combination of two primers 274
After combining the results annotated by the two primers, the total food composition of 275 25 zokors was obtained. The food composition of zokors was evaluated by %RA and 276 %FOO, respectively. We excluded species with %RA less than 0.01% from the analysis , 277 on the one hand, because their ASV numbers were too low and their %FOO was too low 278 to introduce bias if they were included in the analysis [51], on the other, the importance 279 of rare food categories may be overestimated if the analysis is conducted only by ordering 280 the occurrence.

281
By combining %RA and %FOO methods, we evaluated the diet of the zokor at the 282 Family and Genus levels. We found that (1) at the Family level, %RA of Asteraceae was 283 38.16%, while %FOO was 100%, indicating that every zokor fed on Asteraceae. %RA of 284 Poaceae was 22.52%, %FOO was 96%, and only one zokor did not feed on grassy weed. 285 This indicates that the Asteraceae and Poaceae play an extremely important role in 286 the zokor's diet, next for Rosaceae, Pinaceae, Brassicaceae and so on. The tendency of 287 %RA and %FOO at the Family level showed relatively good consistency. (Fig.4);(2) at 288 the genus level, the %RA of Echinops is 14.45%, which is the highest among all plant 289 genera, but the %FOO is only 56%. The %RA and %FOO of Littledalea are 11.86% and 290 40% respectively. %RA of Artemisia was 9.85%,surprisingly, %FOO was the highest 291 among all plant genera, 84%. In addition, some plant genera with relatively important 292 values of %RA and %FOO include Cirsium, Elymus, Potentilla, Bupleurum, etc. (Fig.5). 293 The detailed %RA and %FOO tables are shown in the attached tables.(Appendix 1) 294   Fig.4 and Fig.5 respectively showed that at the Family level, the zokor mainly fed 295 on plants of Poaceae, Asteraceae and Rosaceae, at the Genus level, the zokor mainly fed 296 on plants of the Echinops, Artemisia, Picea and Cirsium. Fig.6 and Fig.7 respectively 297 showed the diet of zokor at the Family and Genus level by %RA method. Circos diagram of food composition of zokor 299 We produced Circos diagram to show the detailed correspondence between the diet of 300 zokor in different groups. Fig.8 is the correspondence analysis of diet of two species of 301 zokor at the Genus level of food, and Fig.9 is the correspondence analysis of zokor's 302 diet at the Genus level in 7 plots.  As Sobs index, Chao 1 index and Ace index have almost no difference in the calculated 306 diversity value, so only the Sobs index was selected for analysis. Coverage index was 307 always 1 because ASV which only appeared once was removed in the sequence annotation, 308 so the analysis of the Coverage index was not carried out.

309
At the Family level, the zokor's diet Sobs index ranged from 3.67-7.25 (mean 5.75), 310 and Shannon ranged from 0.32 to 0.82 (mean 0.54), indicating low diversity. At the 311 Genus level, the zokor's diet Sobs index ranged from 9.33 to 14.40 with an average value 312 of 11.35, and the Shannon index ranged from 0.36 to 1.10 with an average value of 0.86, 313 indicating that the diversity of diet was low.
The diversity of zokor's diet differs greatly among different study areas, which may 315 be related to the plant species distributed in different study areas. Shannoneven index 316 reflects the evenness of eating habits of zokor, which has the same trend as the Shannon 317 index. In the study area with a high Shannon index, the evenness of food species is also 318 high. (Table 3.) 319 In comparison of food species between E. smithii and E. cansus, the food diversity 320 of E. cansus was slightly higher than that of E. smithii in terms of Sobs index, Shannon 321 index and Shannoneven , but the sample size was smaller than in E. smithii. (9<16). 322 Zokor was the dominant species in this study area. Our results suggest that the E. 323 cansus avoids interspecific competition with its relative E. smithii by adopting a more 324 diverse diet. Due to its large population, the E. smithii consumes more available food 325 resources in this area. Kruskal-Wallis rank sum test and FDR multiple test correction were used to test the 328 difference between index groups, and there was no significant difference between Family 329 and Genus and among different species of zokor (P >0.05) (data not shown).

330
Correlation analysis of diet alpha diversity indexes and body indexes 331 We measured the body indexes of zokor, including body weight and body length, and 332 used Pearson double-tailed significance test to find the correlation between diet diversity 333 index (Sobs index, Shannon index, Shannoneven index) and body indexes. The results 334 showed that the body weight and body length of zokor were slightly negatively correlated 335 with the alpha diversity indexes, and there was no significant difference. The correlation 336 coefficient between Sobs index and body weight was -0.382 (P =0.059). Our results 337 seem to indicate that diet diversity may decrease with greater body size. Larger zokor 338 are better at obtaining food than smaller ones, so the variety of food available may be 339 relatively simple.  Fig S2) These results may indicate that (1) the plant species in our study area 345 may have low heterogeneity, and the zokor's diet in different study areas was mostly 346 clustered together; (2) The interspecific competition between the two zokor species may 347 not be strong, and the food between the two species is similar.

348
Analysis of food selectivity 349 Through investigating the research area of the plant resources, form a list of plants in the 350 study area (Appendix 2), compare the zokor feeding selectivity index of sample area 351 plant selection in the sample area coverage of more than 0.1% of the plant, comparative 352 DNA metabarcoding on the plants of the Genus of zokor feeding ratio, selectivity index 353 calculation. Ivlev's [52] selectivity index was used to calculate the plant selectivity of 354 zokor. E i = (r i -p i )/(r i + p i ); r i : refers to the proportion of i plants in the feeding 355 composition. p i : refers to the proportion of i plant biomass to the total biomass in the 356 plant survey quadrate. E i is between -1 and +1, if E i >0 indicates that the animal has 357 a positive choice for the plant. E i <0 means negative selection, and E i = -1 means no 358 selection.

359
It can be seen that the zokor feeds positively on the genera Echinops, Bupleurum, 360 Crisium, Brassica, Calamagrostis, Medicago, Sanguisorba and other Poaceae, which 361 account for 21.87% of the plant resources in this area, but 53.22% in the food of the 362 zokor. For the genera Elymus, Leymus, Artemisia, Poa, Potentilla, and Convolvulus they 363 were widely distributed in this area, accounting for 74.15% of the plant resources, but 364 only 16.43% of the zokor's diet. (Fig.10  Plants are classified in terms of their feeding parts 375 We briefly discussed and analyzed various types of food according to the plant parts that 376 zokor feeds on. The rat mainly feeds on taproots(54.60%) and rhizomes(26.41%), (Table 377 4.) which are rich in water and energy, which can meet the high energy consumption 378 caused by zokor's digging habit [50]. In this study, the zokor's diet endemic to China was determined by DNA metabarcoding 381 based on HTS, which was the first report on the food composition of E. smithii. Our 382 results suggest that DNA metabarcoding is feasible for studying the feeding habits of 383 rodents that live entirely underground and feed mainly on plant roots. In the study area, 384 E. smithii is sympatric with E. cansus, and both belong to Eosplax. Studies have shown 385 that the diet of seven species of Ctenomys in Brazil has a high overlap [4].Our results 386 also confirmed that there was no significant difference in food composition between 387 E. smithii and E. cansus, with schoener niche overlap index of 0.69. In this area, 388 two species of zokor mainly fed on the roots of PH, which mainly belonged to genus 389 Echinops,Artemisia,Cirsium, Bupleurum and Elymus, and also fed on the roots of Picea 390 and Pyrus. Since the samples we selected were the gastric contents of zokor and ITS 391 primers were used to amplify the DNA of the gastric contents, many fungi were also 392 amplified because the digestion degree of food in the stomach of mammals was the 393 lowest, but the main groups were common plant pathological fungi, endophytic fungi, 394 intestinal fungi, etc. (data not shown). It's either eaten by the rat as it feeds on the 395 roots, or it's already present in the stomach, which may not be the result of the zokor's 396 choice of food. This was also the case in the metabarcode feeding studies of lemmings' 397 stomach contents [28].

398
Interestingly, ITS also amplified Canis. Since there were no stray dogs distributed in 399 the study area, it was preliminarily speculated that it might be Canis lupus, which may 400 be that zokor ingested wolf feces in the process of digging and eating. The common sense 401 now is that birds of prey (Falconiformes,Strigiformes) and small predators (Canidae, 402 Mustelidae,Felidae) are natural enemies of the zokor, and it's easy to understand that 403 they are, and the zokor is also a rat. But this view is one-sided, zokor rarely appears 404 on the ground. Its temporary herbivorous burrow is 10-30cm from the ground, and its 405 permanent cave is 50-210cm from the ground [13]. Whereas raptors, who hover high in 406 the air and use their keen vision to spot their prey, apparently couldn't spot the zokor, 407 there is no study or report on prey of zokor by raptors. Small, slender carnivores, such 408 as the Mustela sibirica, that can burrow into the zokor's hole to prey on the zokor, are 409 more plausible, but research on this is lacking. Further research is needed to determine 410 whether the wolf we have found indicates that it is a natural enemy of the zokor.

411
And we found something interesting in the species of plants that we annotated, 412 Hydrodictyaceae(EC-2,ES-11,EC-6), Polypodiaceae (ES-2), these plant species are unique 413 in that they will only live where there is watery or humid. Streams in the study area 414 confirmed the possibility of the Hydrodictyaceae. Both EC-2 and EC-6 belong to LD-F 415 plot, while ES-11 and ES-2 belong to LD-D plot, and their distance from streams is 180m 416 and 600m, respectively. This suggests that the zokor is adapted to moist conditions, and 417 that it travels considerable distances in search of food.
found that zokor is a typical opportunist, that is, it feeds on most of the plants in the study 420 area, which is consistent with many studies on the feeding habits of herbivores [4,21,24,26]. 421 Since most of the Picea, Pinus mongolicus and Pinus tabulaeformis in the study area 422 have carried out wire nets in their roots, which can effectively prevent the biting of 423 zokor, they then turned to feed on the rhizome of PH such as Poaceae and Asteraceae 424 widely distributed in the area. The zokor's food storage showed that almost entirely the 425 asexual part of the plant.

426
It was autumn when we caught the rat, which is the season for storing food. We also 427 found two storehouses, one of which had five small pantries, and the other only three 428 small pantries, which were neatly stuffed with plant roots. It can be mainly divided 429 into gramineous weeds and forbs. Gramineous weeds account for 28.18%-34.23% of 430 raw weight, and forbs account for 64.84%-69.74%. Different from the E. baileyi [21], 431 the zokor warehouse we found contains only the hypogeal parts of plants, but not the 432 aerial parts. It has been found that the zokor spends part of its time on the ground, 433 which may allow it to feed on the aboveground parts of plants [59]. In autumn, the 434 density of Ctenomys was positively correlated with the availability of plants with reserve 435 organs [60]. Similarly, the zokor becomes significantly more active in autumn, actively 436 searching for plant reserve organs (rhizomes, tubers, taproots, etc.). In food selection, 437 taproot (54.60%) and rhizome (26.41%) were the main food types, and these types of 438 roots provided water, carbohydrate, fat and protein for zokor. For example, dandelion 439 roots contain 22.59% more carbohydrates than aerial parts [57]. The root tuber of 440 Stachys geobombycis contains 23g of sugar, 4.1g of protein, 0.3g of fat and nearly 72% of 441 water per 100g [58]. 442 Studies have shown that grassland management such as grazing and mowing to 443 reduce litter accumulation can alleviate the negative impact of nitrogen deposition on 444 plant diversity by reducing the asexual reproduction of dominant species [56]. Our study 445 area is dominated by PH, and the zokor's feeding on the asexual reproductive organs of 446 these herbaceous plants may reduce the asexual reproduction of dominant species and 447 indirectly slow down the decrease of plant diversity. In our study area, and in a larger 448 area of woodland, the zokor is considered a pest mainly because it gnaws on the roots 449 of afforestation plants. The contribution of woodlands to soil and water conservation 450 is well known, however, the role of the zokor in the ecosystem is little-known. Only 451 when the density of zokor was too high, it was forced by intraspecific and interspecific 452 competition to cause harm to forest. Most of the tree species in the study area were 453 covered with wire nets, and our results and field survey showed that zokor did little 454 harm to these trees. According to the literature, zokor can optimize the composition of 455 meadow community with appropriate density [22]. At the same time, as the basic species 456 in the food web, zokor plays a positive role in maintaining the stability of the ecosystem 457 and promoting energy flow. We should not destroy zokor blindly and unilaterally. Our 458 results showed that the zokor did not cause great harm to the tree species in the study 459 area. It is not that the zokor's habits have changed,outside the study area, there are 460 still many pines that were eaten to death. In the subsequent afforestation, we should 461 set up the idea of protecting the target species and tolerate the existence of zokor to a 462 certain extent, which is a win-win situation for the development of forestry economy 463 and ecological restoration.

464
At present, most studies on diet based on HTS still choose OTU clustering method 465 for sequence annotation [4,18,20,24,[26][27][28]. We found that the ASV method was similar 466 to the OTU clustering method in the identification ability at the Order, Family and 467 Genus levels, but at the Species level, ASV method had 16 more species than OTU. As 468 with fungal and microbial studies [35,36], the ASV method also improves the diversity 469 of herbivore analysis.

5).
In this study, 65 species of plants belonging to 56 genera and 24 families were 472 investigated in the study area, and 154 species of plants belonging to 80 genera and 473 32 families were detected by DNA metabarcoding technology.(Appendix 3) It can 474 be directly found that DNA metabarcoding can greatly improve the identification of 475 food species and improve the understanding of zokor's diet. There may be two reasons 476 why zokor fed more species than plants in the study area: (1) As a typical burrowing 477 animal, zokor mainly feeds on the roots of plants, possibly taking in some soil in the 478 process. In the study on the feeding habits of Ctenomys, it was found that 23 plant 479 families and 58 MOTUS could be recovered by extracting DNA from soil alone, and the 480 diversity of plant families and MOTUS in soil samples was higher than that in the feces 481 of Ctenomys [4]. Therefore, there may be false positives in our experimental results, that 482 is, the actual variety of food resources consumed may be smaller than that recovered by 483 DNA metabarcoding. (2) The time when we investigated the plant species in the study 484 area was autumn, when vegetation and trees were dying, and some species were not 485 investigated. Indeed, in the Pan species graph, our specimen has not flattened the observed species 487 curve, which requires more samples to be added for food analysis of zokor. In our 488 study area of 300 ha with similar habitat types, we need to collect samples over a large 489 geographical range and a long time period to avoid bias [63]. Fortunately, our samples 490 were not collected at the same time or within a few days. Our sampling time span 491 was nearly two months. In general, larger sample sizes and technical repetitions of a 492 single sample (multiple extraction, multiple sequencing) are certainly beneficial [64,65]. 493 However, due to the increasing workload and cost, researchers cannot continuously 494 increase the sample size and increase the repetition. As a result, they need to find the 495 best balance between sample size and biological and technical replication, taking into 496 account the workload and cost [63]. The results of our study are credible at the level of 497 Genus classification and can explain the dietary characteristics of zokor under natural 498 conditions. Many herbivorous animal feeding studies are also conducted at the level of 499 Family [30]and Genus [26,66]. In the study of feeding habits of DNA sequencing, there 500 may be deviations/errors in every link, such as sample preservation, sample contamination, 501 DNA extraction, PCR amplification, primer selection, library preparation, selection of 502 sequencing platform, error removal, sequence taxonomic assignment, etc., which will 503 affect the final result. Even so, the classification accuracy of diet studies based on HTS 504 is much higher than that of traditional microscopical methods, and the workload is much 505 less than that of microscopical methods.