Expansion of Cyclophyllidea Biodiversity in Rodents of Qinghai-Tibet Plateau and the “Out of Qinghai-Tibet Plateau” Hypothesis of Cyclophyllideans

The Cyclophyllidea comprises the most species-rich order of tapeworms (Platyhelminthes, Cestoda) and includes species with some of the most severe health impact on wildlife, livestock, and humans. We collected seven Cyclophyllidea specimens from rodents in Qinghai-Tibet Plateau (QTP) and its surrounding mountain systems, of which four specimens in QTP were unsequenced, representing “putative new species.” Their complete mitochondrial (mt) genomes were sequenced and annotated. Phylogenetic reconstruction of partial 28S rDNA, cox1 and nad1 datasets provided high bootstrap frequency support for the categorization of three “putative new species,” assigning each, respectively, to the genera Mesocestoides, Paranoplocephala, and Mosgovoyia, and revealing that some species and families in these three datasets, which contain 291 species from nine families, may require taxonomic revision. The partial 18S rDNA phylogeny of 29 species from Taeniidae provided high bootstrap frequency support for the categorization of the “putative new species” in the genus Hydatigera. Combined with the current investigation, the other three known Taeniidae species found in this study were Taenia caixuepengi, T. crassiceps, and Versteria mustelae and may be widely distributed in western China. Estimates of divergence time based on cox1 + nad1 fragment and mt protein-coding genes (PCGs) showed that the differentiation rate of Cyclophyllidea species was strongly associated with the rate of change in the biogeographic scenarios, likely caused by the uplift of the QTP; i.e., species differentiation of Cyclophyllidea might be driven by host-parasite co-evolution caused by the uplift of QTP. We propose an “out of QTP” hypothesis for the radiation of these cyclophyllidean tapeworms.


INTRODUCTION
Cestoda is a class of parasitic worms in the flatworm phylum Platyhelminthes that parasitize the intestines of all major groups of vertebrates, including fish (Kuchta et al., 2020), reptiles (Yudhana et al., 2019), birds (Okulewicz, 2014), and mammals (Carlson et al., 2020), and cause severe, mild or no symptoms of infection. Their larvae usually achieve development in one or two intermediate hosts (Kuchta et al., 2020) and can cause severe symptoms and even death in animals and humans; species of Taenia and Echinococcus cause the highest health and economic impact.
To date, there are approximately 4,800 described species in the class Cestoda, belonging to 833 genera and 19 orders, of which the Cyclophyllidea is the most species-rich order in the class, with more than 3,100 valid species, distributed among 16 families (Sharma et al., 2016;Caira and Jensen, 2017;Kuchta et al., 2020). Each free-living metazoan species is considered to harbor at least one protozoan or metazoan parasite species (Poulin and Morand, 2004) with which they interact and usually co-evolve (Ebert and Fields, 2020). Parasites are probably one of the largest groups of undescribed organisms, as most are cryptic whether in their parasitic or free-living forms (Dobson et al., 2008;Larsen et al., 2017;Okamura et al., 2018;Carlson et al., 2020). It is estimated that there are a total of about 100,000-350,000 helminth species in vertebrates around the world, of which 85-95% are undiscovered or recorded to science (Carlson et al., 2020).
Many of the world's biodiversity hotspots are located in large mountain systems, and their role in the evolutionary diversification of organisms is manifold (Hoorn et al., 2018;Rahbek et al., 2019a,b). The Qinghai-Tibet Plateau (QTP) and its surrounding mountain systems of the Eurasian continent, have yielded arguably the biggest and probably the most biologically diverse area of montane species (Päckert et al., 2020). In the past decade, several new Taeniidae species have been described from wild rodents on the QTP (Xiao et al., 2005;Wu et al., 2021). In terms of species diversity and sheer population sizes, rodents are perhaps the most important intermediate and definitive hosts of tapeworms, and are the most widely distributed and diverse group of mammals; about 43% of all species (Singla et al., 2008;Wu et al., 2018). Changes in climate and vegetation caused by the uplift of the QTP may have promoted local adaptations such as the evolution of cold-and hypoxic-tolerant rodent species, which in turn may have led to the co-evolution and radiation of their parasites (Wang et al., 2018;Wu et al., 2021). We hypothesize that there may be many undiscovered tapeworm species parasitizing the rodents in the QTP. Although the number of tapeworm species has been underestimated in general, some practices may lead to the erroneous proposal of new species (Carlson et al., 2020;Kuchta et al., 2020).
Morphological distinctions have been used for the description of many tapeworms, however, the homoplasy in morphology poses quite a challenge to infer their evolutionary lineage (Scholz et al., 2021), where morphological features may be significantly affected by different intermediate host sources (Lymbery, 1998). Many undescribed species may be genetically different but morphologically indistinguishable. The inability to distinguish such cryptic species affects accurate assessment of host range and estimates of total diversity (Carlson et al., 2020). The 28S and 18S nuclear ribosomal RNA genes (rDNA) are relatively conserved within species and are often used to differentiate different species (Wickström et al., 2005;Nakao et al., 2013;Scholz et al., 2021). The mitochondrial (mt) genome is largely haploid and uniparentally inherited, so their effective population size is four times smaller than that of the nuclear genome (Toews and Brelsford, 2012). Since the process of lineage sorting is inversely proportional to the effective population size, this means mitochondrial (mt) genomes will complete this process faster than nuclear genomes (Toews and Brelsford, 2012). Thus, the genetic nature of a mt genome makes it likely more sensitive than any single nuclear marker to distinguish closely related species and study their phylogenetic relationships (Lee et al., 2007).
Therefore, in this study, we sampled parasites in rodents from the QTP and its surrounding areas. In total, seven Cyclophyllidea specimens were collected and characterized using molecular tools. By comparing new (mitochondrial, nuclear 28SrDNA, 18S rDNA) with published homologs (GenBank) we found four isolates to be markedly different in DNA sequence thus likely representing "putative new species." The complete mt genomes of these unknown taxa were sequenced and annotated, and their taxonomic status was analyzed and verified through phylogenetic reconstruction of five datasets containing a total of 320 species from 10 families. By combining evolutionary divergence time analyses of mt genes of classified cyclophyllideans in NCBI Taxonomy Database 1 and the paleogeography of QTP, we thus speculate an "out of QTP" theory for cyclophyllideans.

Sample Collection
Rodents were live-trapped in meadows in Tibet, Qinghai, Sichuan, Gansu and Xinjiang province or autonomous region of China in 2013 and from 2018 to 2020. Rodents were euthanized and dissected according to the Ethics Statement mentioned below, cysticerci and host livers were collected from the enterocoelia and thorax, and adults were extracted from the intestines. Detailed sample collection data and host identities are described in Supplementary Table 1. After detaching the lesions, parasite specimens and host livers were kept in 75% (v/v) ethanol for molecular identification.

DNA Isolation, Amplification, and Sequencing
DNA samples of hosts and parasites were extracted using Blood and Tissue Kit (Cat. No. 69504,Qiagen,Germany) as instructed by the manufacturer, and were amplified and sequenced for identification by conserved primers of cytb gene of small mammals (Fan et al., 2011) and cox1 gene of tapeworm (Bowles et al., 1992), respectively. By means of highly similar BLAST search in the nucleotide collection (nr/nt) database, 2 four of the Cyclophyllidea specimens with less than 95% identity of the most similar cox1 sequence were identified as putatively unknown species (Xiao et al., 2005). The mt genomes of four of the putatively unknown specimens were sequenced and assembled according to the following procedure: firstly, the DNA of the four species was amplified and sequenced using primers published in Wu et al. (2021); missing fragments were amplified with newly designed primers using the program Oligo 6.0 with the method described in Wu et al. (2021), until entire circular mtDNAs were amplified and sequenced. A list of primers for each species can be found in Supplementary Table 2. The 18S and 28S rDNA fragments of these four species were also amplified and sequenced with conserved primers Yan et al., 2013) for further species identification and phylogenetic analyses. Primers were synthesized by Tsingke Biotechnology (Xi'an, China). Standard 25 µl PCR protocol was used to amplify the DNA fragments. The PCR products were purified and sequenced according to protocols in Wu et al. (2021). SeqMan software was used to assemble the cytb gene sequences, 18S and 28S rDNA partial sequences and mt genomes (see Supplementary Table 3 for GenBank accession nos.).

Mitochondrial Genome Annotation
The four new mtDNAs were annotated preliminarily by Geseq 3 with the reference of the most closely related species (Mesocestoides corti, Anoplocephala magna, Hydatigera krepkogorski, and Moniezia expansa) identified by the phylogenetic analyses in Figures 1, 2, whose mt genome annotations are available in GenBank. Putative tRNA genes were verified using ARWEN 4 using default parameters (Laslett and Canbäck, 2008). The positions of their open reading frames (ORF) and rRNA genes were further checked and modified using SnapGene (v3.2.1) based on alignments with the reference of the most closely related species mentioned above. SnapGene (v3.2.1) was used to translate the protein-coding genes (PCGs) into their amino acid sequence with echinoderm/flatworm mitochondrial code (NCBI translation table 9) and to illustrate the annotated mt genomes.

Phylogenetic Analyses and Sequence Identity
As ingroups for phylogenetic reconstruction, we combined mt genes and 18S or 28S rDNA sequenced in this study with those of other classified cyclophyllideans available in GenBank. Three species (Caryophyllaeus brachycollis, Anindobothrium anacolum, Spirometra erinaceieuropaei) of Eucestoda belonging to different orders were chosen as outgroups.
The 28S rDNA fragments of classified cyclophyllideans in NCBI Taxonomy Database are available for most families in GenBank except Taeniidae, while that of the 18S rDNA were only available for Taeniidae, so here the evolutionary trees of 28S and 18S rDNA fragments were constructed separately. Based on the results of the 28S rDNA phylogeny, a simplified 28S evolutionary tree was reconstructed without affecting the topological structure of the tree by removing the species in the same genus and clade. To complement and validate the inferred evolutionary relationship between Cyclophyllidea species, except the species of family Taeniidae, the phylogenies with cox1 and nad1 fragments were reconstructed.
In summary, 5 datasets containing 320 species in 10 families were assembled for phylogenetic analyses (see Supplementary  Tables 4-8 for GenBank accession nos. and species, genera and families used in each dataset). Since the limited availability of data on different genes of the same species, we performed separate phylogenetic reconstruction for each dataset to cover more species. Datasets were aligned using MAFFT v7.487 with auto option (Katoh and Standley, 2013). The alignments were trimmed by using Trimal v1.2 under the automated1 option (Capella-Gutiérrez et al., 2009) to preserve the same sequence regions and exclude the ambiguously aligned sites. All phylogenetic trees were constructed with maximum likelihood (ML) inference using IQ-TREE v2.1.4 (Nguyen et al., 2015) with ultrafast bootstrap 1,000 replicates in Ubuntu 20.04.2 LTS operating system, where the best-fit models were automatically selected by ModelFinder (Supplementary Table 9; Kalyaanamoorthy et al., 2017) and the best number of threads were also selected under AUTO option (using the command line: iqtree -s alignment_file -m -MFP -bb 1000 -nt AUTO). All other parameters were set to their default values. The ML trees were then visualized on the IToL web server (Letunic and Bork, 2016).
The related species of the putatively unknown species were identified based on phylogenetic reconstruction, and the percentage identity of sequences in the 5 datasets between the putatively unknown species and their related species was calculated by BLAST (Supplementary Table 10).

Analysis of Divergence Times
Species divergence times were calculated using BEAST v2.6.2 (Bouckaert et al., 2014) based on two datasets without shared species: all available mt PCGs dataset of 54 species (Supplementary Table 11) and cox1 + nad1 fragments dataset of other classified 54 species (Supplementary Table 12). These two datasets were aligned and trimmed with MAFFT v7.487 and Trimal v1.2 as described above, and were partitioned according to different genes. Model selection of partitions was identified by Partition Finder v2.1.1 (Lanfear et al., 2012) with the set of "linked" branch lengths, "beast" models, "aicc" model selection, and "greedy" search. Partition schemes and substitution models can be found in Supplementary Table 9. The Strict Clock model was chosen to ignore the rate differences between the branches in the mode and the gamma category count was set to 4. Other settings, such as substitution rate and shape, in the site model were evaluated in the analysis. The Calibrated Yule model (Heled and Drummond, 2015) was used as the tree prior, as it is a simple model of speciation that is generally appropriate when considering sequences from different species. Time calibration was calibrated with the divergence date between T. saginata and T. asiatica (∼1.14 Mya) and divergence date between Schistosoma japonicum and S. mansoni (∼56.10 Mya), which agreed with reported fossil evidence from shark coprolites and previously estimated dates based on mt genes (Wang et al., 2016). Posterior probability estimates were sampled every 1,000 iterations over a total of 10,000,000 iterations per MCMC run. Other options were run on their default values. Tracer (v1.7.1) was used to summarize posterior probabilities. Trees were annotated via TreeAnnotator (v2.1.2) using a maximum clade credibility tree and median heights settings with 10% burn-in. The number of divergence nodes in every 2 Mya was summarized based on the evolutionary divergence time trees of two datasets, respectively.

Species Identification and Phylogenetic Relationships
The hosts of parasites were confirmed by BLAST searches (Supplementary Table 1; see Supplementary Table 3 for GenBank accession nos. of host cytb gene). Two of the parasite species were cysticerci and adult worms, each retrieved from the abdominal cavity and intestinal tract of two vole hosts (Neodon leucurus) collected from the same pasture location near the Shigatse City of Tibet Autonomous Region (Supplementary Table 1). Phylogenetic analyses shown in Supplementary Figure 1 and Figures 1A-C confirmed that these two species likely belong to the genera Mesocestoides and Paranoplocephala, and were labeled as Mesocestoides sp. RKZ08 and Paranoplocephala sp. RKZ13, respectively. The cysticerci collected from the liver of plateau pika (Ochotona curzoniae) from Xietongmen county of Tibet and zokors (Eospalax fontanierii) from Xiahe county of Gansu province were confirmed to be the same species by alignment of cox1 and 18S rDNA segments, which was in the monophyletic group of the genus Hydatigera in the phylogenetic trees (Figure 2), and was marked as Hydatigera sp. XHPW10. The final species was an adult tapeworm collected from the intestine of plateau pikas from Shiqu county, Sichuan province, and occurred within the genus Mosgovoyia (Supplementary Figure 1 and Figures 1A-C), and was thus named Mosgovoyia sp. SQ20. These four species were identified as newly sequenced and "putative new species" as their cox1 and nad1 sequences having less than 95% identity with available related taxa (Supplementary Table 10; Xiao et al., 2005). The degree of differentiation of their mt genes was higher than that of the nuclear genes 18S and 28S rDNA (Supplementary Table 10), reflecting their differentiation likely results from the so-called deep mitochondrial divergence (DMD, Zhang et al., 2019). In addition, three known cysticerci of Taeniidae, T. caixuepengi, T. crassiceps, and Versteria mustelae, were also identified in the present study (see Supplementary  Table 13 for their cox1 fragments).

General Features of the Mitochondrial Genome of "Putative New Species"
The complete mt genomes of the four "putative new species" were 13,361 bp (GenBank ID: MW808979), 13,730 bp (GenBank ID: MW808980), 14,148 bp (GenBank ID: MW808981), and 13,776 bp (GenBank ID: MW808982) in length. Each of them contains two rRNA genes (rrnS and rrnL) and 12 proteinencoding genes (atp6, cytb, nad4L, cox1-3, and nad1-6), which are arrayed in the typical order of mt genomes of cestodes. They each contain the 22 typical tRNA genes set of cestodes, and share a common set of anticodons. The order of tRNA genes is roughly the same, except between nad6 and nad5 genes. Species Paranoplocephala sp. RKZ13 had a repeat sequence of trnL and trnR in the highly variable region between nad6 and nad5, suggesting that it had an insertion sequence in this region that made its mt genome longer than the other three species (Figure 3 and Table 1).
Flatworms use a unique set of mt code for protein translation (Nakao et al., 2000;Telford et al., 2000). In addition to ATG, GTG was also used as an initiation codon in a small fraction of coding genes in their mt genomes. For the termination codon, all species used only TAA and TAG; TGA was not identified as a termination codon ( Table 1).

Divergence Times Analysis
The divergence time based on the two datasets used for time calibration is consistent with previous genome-based analysis results (Figure 4; Wang et al., 2016). Three of the "putative new species, " Mesocestoides sp. RKZ08, Mosgovoyia sp. SQ20, and Hydatigera sp. XHPW10, might have originated from a similar phase in the late Miocene, while Paranoplocephala sp. RKZ13 diverged during the Pleistocene (Figure 4). By summarizing the number of divergent time nodes in the divergence time trees over time, it was found that the trees of mt PCGs and cox1 + nad1 produced similar differentiation rate trends: there was a slight acceleration of the evolutionary rate in the period 14-24 Mya, and a marked acceleration during the period 4-10 Mya (Figure 5). However, compared with cox1 + nad1, the differentiation rate of mt PCGs was faster in the 14-24 Mya but slower in 4-10 Mya, which may have been caused by the species bias used in both datasets.

DISCUSSION
Given the rich diversity and the large rodent population in western China and the few published reports of tapeworms in rodents, except for Taeniidae (Xiao et al., 2005;Zhao et al., 2014;Wang et al., 2018;Zhang et al., 2018;Wu et al., 2021), the present  Table 1. knowledge of tapeworm biodiversity in rodents in western China suggests far greater biodiversity yet to be uncovered.
Here, the mt genes and 18S or 28S rDNA fragments of four (two larvae and two adults) unidentified Cyclophyllidea species differed from their related species (Supplementary Table 10). However, due to the specimen distortion and insufficient specimen encountered in this study, it is not clear whether they have been described morphologically. All four species showed apparent discordance percentage identity with the related species in their mtDNA and 18S or 28S rDNA, which is the common DMD pattern across the animal kingdom, and demonstrated in Echinococcus granulosus sensu stricto (Kinkar et al., 2017), possibly due to the parthenogenetic inheritance of mitochondria, gene flow and recombination in the nuclear genome (Toews and Brelsford, 2012).
to the recently described cysticerci species in plateau pika from Qinghai province (Wu et al., 2021). Here the T. caixuepengi larva was also found to be parasitic in plateau pikas from Xietongmen, Saga and Sa'gya county of Tibet and Qilian county of Qinghai (Supplementary Table 1), and so far, not found in other sympatric rodent species. Due to the absence of comparative studies, it is currently not clear if this species has preference only for plateau pikas as its intermediate host.
The larvae of T. crassiceps and V. mustelae were, respectively, found in Jeminay and Xinyuan county of Xinjiang autonomous region (Supplementary Table 1), which were also distributed on the northeast QTP Zhao et al., 2014). Wide geographic distributions of identical species suggest that their endemic geographic range should be far beyond the available survey data. Except for Mesocestoides sp. RKZ08 identified in this study, the Me. litteratus is another species of the genus Mesocestoides reported in Qinghai and Heilongjiang provinces of China (Wang et al., 2006;Li et al., 2013). However, there is no   Figure 4A. The blue curve represents the change of divergence nodes number based on Figure 4B.
previous record of Paranoplocephala spp. and Mosgovoyia spp. available in China except for Paranoplocephala sp. RKZ13 and Mosgovoyia sp. SQ20 found in this study. Complete mt genomes of the four "putative new species" were sequenced and annotated, and the sequences were clearly different from all available mt genomes sequences; however, they were similar in length, gene order and composition as other cyclophyllideans with respect to rRNA, tRNA, and proteinencoding genes (Le et al., 2000;Jeon et al., 2005Jeon et al., , 2007Wu et al., 2021). Different arrangement of genes occurred in tRNA genes between nad5 and nad6 genes, but was consistent with their respective most relative species, matching the expectation that the arrangement of mt genes partly determines the genetic relationship of parasites (Gazi et al., 2016). Moreover, repeat copies of tRNA gene between the nad5 and nad6 genes in the mt genome of E. granulosus sequenced by third-generation sequencing have been reported (Kinkar et al., 2019). We also observed a repeat sequence of trnL and trnR genes between nad5 and nad6 genes of the Paranoplocephala sp. RKZ13 mt genome (Figure 3 and Table 1), which indicates that there may be hidden tRNA gene repeats in the mt genome of tapeworms that are hard to identify due to errors in PCR amplification and Sanger sequencing, techniques which often fail to recover repeat regions (Kinkar et al., 2019).
Genetic drift and adaptive differentiation between allopatric populations is responsible for most speciation amongst plants and animals (Turelli et al., 2001). For parasites, however, host association is a key driver in their evolution. Host switching among sympatric populations may lead to ecological isolation, so sympatric speciation of parasites is common (de Meeûs et al., 1998;Paul, 2002;Huyse et al., 2005;Wu et al., 2021). These two models of evolution are not in conflict: the adaptation of parasite to specific host is like the adaptation of animal to specific living environment; so, the evolution of the host, especially its immune system, may be viewed as equivalent to the change of the living environment for the parasite. The environmental and climatic changes caused by the uplift of the QTP are a major driving force for the evolution of associated biotas (Favre et al., 2015). Paleobotanical data suggest that the southeastern margin of the QTP was dominated by a warm and humid subtropical or tropical climate during the Miocene due to the influence of South Asian and East Asian monsoons (Sun and Wang, 2005;Jacques et al., 2011). Since the mid-Miocene, the significant rise of the Himalayas and the Tianshan Mountains, together with worldwide cooling, incurred dramatic changes in air circulation, leading to gradual aridification of the QTP and Central Asia (Miao et al., 2012;Favre et al., 2015). Finally, in the Late Miocene and Early Pliocene, QTP uplift resulted in the accumulation of global ice and the eventual disappearance of the Tethys Sea, which also contributed to the drying of Central Asia (Lu and Guo, 2013).
These timelines of climatic and environmental changes caused by the uplift of the QTP are highly consistent with the timelines of the differentiation rate of Cyclophyllidea species analyzed in this study (Figures 4, 5) and that of plateau pika analyzed in Wang et al. (2020). We speculated in this study that the differentiation of cyclophyllideans may have been driven by host evolution caused by the uplift of the QTP. During the tropical period of QTP, the optimum living environment created the biological diversity (Cai et al., 2020;Päckert et al., 2020), and the species of Cyclophyllidea gradually differentiated. With the rapid uplift of the QTP, the environment changed into a dry and cold climate (Lu and Guo, 2013), and species differentiation of Cyclophyllidea was accelerated by the rapid adaptive evolution of their hosts and geographical isolation caused by the radiation of hosts to the Palaearctic (Favre et al., 2015;Xing and Ree, 2017;Päckert et al., 2020). Finally, in the last 2 million years, Cyclophyllidea differentiation demonstrated an accelerated diversification based on cox1 + nad1 divergence tree (Figure 5), which may be related to the evolution and broad distribution of mammals in Eurasia and the associated population expansion and migration of hominids from Africa to Asia (Dennell, 2004;Rohland et al., 2005;Brugal and Croitor, 2007;Turner et al., 2008;Klein, 2009;Terefe et al., 2014;Wang et al., 2016). The Host-parasite Database of Natural History Museum (HPDNHM) found that most tapeworm were prevalent mainly in the Palaearctic, 5 which is also consistent with the viewpoint that the order Cyclophyllidea originated from the QTP. The Nearctic is another major endemic area where Cyclophyllidea species may have spread over land Bridges across the Bering Strait. However, some species parasitizing birds and other economic and companion animals tend to show a global epidemic, which may be due to the long-distance migration of birds and the spread of human trade and activities.
Phylogenetic reconstruction reveals that some classifications of Cyclophyllidea species may need to be redefined. Closely related species of tapeworm parasites often have similar host specificities and life history (Scholz et al., 2021), a pattern common in the HPDNHM and in our evolutionary analyses, and may provide a basis for revising the classification of Cyclophyllidea species; for example, the family Hymenolepididae and Anoplocephalidae can be divided into multiple families (Figures 1A-C), and Thysanotaenia congolensis should be reclassified into family Davaineidae (Figures 1A,C). However, considering that there is still a large number of undiscovered species, which may provide better support for classification, the current classification status is likely to remain for some time.

CONCLUSION
In conclusion, this study expands the biodiversity of Cyclophyllidea in rodents in QTP and its surrounding mountain systems, and suggests an "out of QTP" hypothesis for the Cyclophyllidea, wherein species differentiation was driven by the uplift of the QTP. Although beyond the scope of this study to consider the evolutionary relationships and history of the whole cyclophyllideans, the species analyzed represent 10 of all 16 families, making this the most extensive study of the evolution of Cyclophyllidea order to date. Verifying the taxonomic revision and the "out of QTP" hypothesis requires more sampling and investigation, including data on a wider geographic and host range, and molecular studies uncovering patterns of host-parasite co-evolution.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Ethics Procedures and Guidelines of the People's Republic of China. Animal Ethics Committee of Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences (No. LVRIAEC2012-007).