The study site is located at an altitude of 1,121 m in the eastern Eurasian steppe at the grassland ecological protection and sustainable utilization research station of the Institute of Grassland Research, Chinese Academy of Agricultural Sciences, Inner Mongolia, China (116°32′E, 44°15′N). The site has a semiarid continental climate with a mean annual temperature of −0.1°C and annual precipitation of 350–450 mm. The coldest month is January (mean temperature: −22.0°C, minimum: −41°C), the hottest month is July (mean temperature: 18.3°C, maximum: 38.5°C), and the annual accumulated temperature is ~2,100–2,400°C.

The major soil types of the study site are calcic chestnut and calcic chernozem. The vegetation is typical of this region and is dominated by two perennial grasses, Leymus chinensis and Stipa grandis. Other common species include the perennial plants Cleistogenes squarrosa and A. frigida and the annual plants Chenopodium glaucum and Salsola collina. The study site was not grazed from 2007 to 2014. Sheep grazing began in 2014, the same year as insect collection. The C. albonemus investigation and transcriptome were conducted in 2015, from early July to middle August.

A relatively flat area with homogenous soil conditions was enclosed in 2014, and a total of sixteen 1.33-ha enclosures (100 × 125 m) were constructed with 1.5-m high iron netting to prevent the movements of sheep in and out of the enclosures. Fifteen enclosures were randomly assigned to five treatments: control, light grazing, moderate grazing, heavy grazing, and overgrazing with three replicates per treatment (one enclosure held the meteorological station). According to local grazing management practices, the number of sheep in the control, light grazing, moderate grazing, heavy grazing, and overgrazing plots were 0, 4, 8, 12, 16 per enclosure, which corresponding to 0, 2, 6, 9, 12 individuals per hectare, respectively (Schönbach et al., 2011). Sheep were allowed to graze in the enclosures from June 10, 2014 to September 10, 2015. All animals used in the experiment were provided by the Chinese Academy of Agricultural Sciences.

Vegetation was sampled at maximum biomass in early August 2015. In five randomly selected quadrats (1 × 1 m) within each plot, we evaluated the following attributes of each plant species: cover (estimated visually), height (determined by using measuring tape), density (number of plants/quadrat), and biomass (g/m²). Above-ground biomass was measured by clipping standing plant material to 1 cm above ground level using shears. The litter was combed out, and the plants were separated by species, stored temporarily in paper envelopes, and then dried in the laboratory for 48 h at 80°C to obtain the dry weight. Plant samples were not taken within about 10 m of the enclosure wall to avoid the effect of heat from the galvanized iron netting.

We evaluated insect diversity from 2014 to 2015. In Xiwu Banner, Inner Mongolia, C. albonemus occurs from late June to late August. According to its life cycle dynamics, samples were collected three times from early July to mid-August (every 20 days) in 2015 corresponding to C. albonemus's early, middle, and late stages. We used the sweep method, which collects all insects, in each of the 15 plots to estimate insect species richness and abundance. In each plot, 200 nets were used to collect insects, and each collection was replicated three times. Insects were collected at least 10 m from the plot boundary to minimize edge effects. We checked the nets visually to ensure that we gathered all insects. The insects were collected only under favorable conditions (sunny days with minimal cloud cover, calm or no wind) from 09:00 to 15:00 h, and the plots were randomly sampled. The contents of the nets were preserved in Ziploc bags. All individuals identified as C. albonemus were counted.

Grasshoppers were randomly selected from samples collected in late July in each of the five grazing intensity treatments for a total of 20 samples (2 females and 2 males per treatment × 5 treatments = 20 samples). Grasshoppers were frozen in liquid nitrogen (Air Liquide, Voyageur 12) and stored at −80°C until RNA extraction.

For each of the five grazing treatments, equal amounts of body tissue (head, thorax, abdomen, legs, and ovaries) from the four individuals were combined and homogenized. Total RNA was extracted using TRIzol reagent (Invitrogen, CA, USA) following the manufacturer's instructions. RNA quality (degradation and contamination) was determined by agarose gel electrophoresis, and purity was determined by using a NanoDrop™ 2000 spectrophotometer (Thermo Fisher Scientific). RNA concentration was determined by using a Qubit H 2.0 Fluorometer (Life Technologies), and RNA integrity was determined by using an Agilent 2100 Bioanalyzer (Agilent Technologies). RNA was extracted using RNAprep pure Tissue Kit (TIANGEN Biotech Co., Ltd., China). The RNA samples were enriched for mRNA using magnetic beads conjugated to oligo (dT) and fragmented into 400- to 600-bp fragments, which were used as a template for first-strand and second-strand cDNA synthesis. The double-stranded cDNA was then purified using AMPure XP beads. After end repair of the double-stranded cDNA, a poly A tail was added, followed by the ligation of sequencing adapters. The cDNA fragments were then selected based on size (150–200 bp) using AMPure XP beads and amplified by PCR. The PCR products were purified using AMPure XP beads to generate the cDNA libraries, which were sequenced using the Illumina HiSeq 2000 platform and the NGS Fast DNA Library Prep Set for Illumina. The paired-end method was used, and the sequencing read length was 200 bp.

The G+C content was measured for each sequencing cycle to determine whether the A+T and G+C levels differed. Using the Illumina HiSeq™ 2000 platform, the relationship between the sequencing error rate (e) and base quality (Q_phred) can be described as follows:

The relationship between base call accuracy and Phred score was calculated using Illumina Casava version 1.8.

To produce clean reads for subsequent assembly and analyses, adapter sequences and low-quality data were removed from the raw data as follows: (1) remove/trim the adapters, (2) discard data for which the percentage of Ns (bases that could not be identified) exceeds 10%, and (3) discard low-quality data (for which the percentage of Q_phred < 5 bases exceeds 50%).

Transcript assembly was carried out using Trinity software version v2012-10-05 (Iyer et al., 2011) at min kmercov = 2; default settings were used for the remaining parameters. The assembly process was previously described (Grabherr et al., 2011). The sequences assembled by Trinity were used as reference sequences for the subsequent analysis. Trinity combined reads with a certain length of overlap to form longer fragments without N (contigs), which were subjected to further sequence clustering to form longer sequences without N. For each gene, the longest assembled sequence [more than one assembled sequence (transcript) for each gene] was regarded as a unigene.

BLAST searches against the NCBI non-redundant-redundant (NR) and nucleotide sequence (NT) databases, SWISS-PROT, PFAM, KEGG, and KOG were performed with a cut-off of 1e^-5. GO terms were assigned using Blast2GO version 2.5 (Götz et al., 2008) by searching the NR database.

Data from each treatment were combined for the reference assembly. Supplementary Table 8 lists the software and parameters used for non-reference transcriptome assembly and analysis.

Clean reads for each sample were mapped onto the reference transcriptome using RSEM software (Li and Dewey, 2011). The read count for each gene was converted to FPKM using the estimation method (Mortazavi et al., 2008). To verify the expression profile of each sample, an FPKM density plot was generated. To analyze read count data and identify differentially expressed genes under different grazing intensity treatments, FPKM in the different grazing intensity plots were compared using DEGSeq and a cutoff value of p adj < 0.005 (Mortazavi et al., 2008).

GO enrichment analysis of the differentially expressed genes was carried by using the GOseq procedure, which is based on the Wallenius non-central hypergeometric distribution (Young et al., 2010), to adjust for gene length bias.

To identify the main biochemical and signal transduction pathways in which differentially expressed genes were involved, pathway enrichment analysis was performed using the KEGG database. KEGG items were mapped using hypergeometric test (Young et al., 2010). FDR corrections were performed for correcting q-value using Benjamini and Hochberg (Shringarpure, 2012). Downstream products of differentially expressed genes were evaluated to identify the substrate associated with the response to grazing intensity. The differentially expressed genes were filtered with q < 0.005 & |log2 (foldchange)| > 1 (Supplementary Table 8).

Principal component analysis was used to evaluate relationships between grasshopper number, vegetation variation, and grazing intensity using the program CANOCO 4.5 (Ter Braak and Smilauer, 2002). Experimental plots were sampled for grass-biomass per species data. The Monte Carlo permutation test was used (number of permutations 999, full model) to indicate the main factors and correlation. The principal component analysis plot (Figure 1) was constructed using CanoDraw 4.5. Other figures were constructed using Microsoft Excel and Origin 8.0. Transcriptomic analyses parameters were listed in Supplementary Table 8.

The dominant plant species, grass cover (%), plant height (cm), and plant biomass (g) differed significantly (P < 0.05) between the following five grazing intensity levels: control (no grazing), light grazing (3 sheep/hectare), moderate grazing (6 sheep/hectare), heavy grazing (9 sheep/hectare), and overgrazing (12 sheep/hectare). In addition, grazing intensity influenced plant succession (Figure 1). Plant species that were dominant in the control plots (L. chinensis, S. grandis, and C. squarrosa) decreased significantly in biomass with increased grazing intensity (Figure 1). The shift in dominant plants occurred under moderate grazing intensity, with L. chinensis, S. grandis, and C. squarrosa replaced by Carex korshinskyi and Artemisia sieversiana as the dominant plant species in the heavy grazing and overgrazing plots. Compared with the control (no grazing) plots, species richness increased in the light grazing and moderate grazing plots and then decreased significantly in the heavy grazing and overgrazing plots (Figure 1).

The total insect species diversity significantly declined along the grazing intensity gradient and over time (from 2014 to 2015), showing a significant linear relationship with sheep density (Figure 2A). C. albonemus abundance also varied according to grazing intensity (Figures 2B,C), demonstrating that grazing had a significant impact on the temporal dynamics of C. albonemus. In the five grazing intensities, C. albonemus abundance was generally lowest in the early sample (early July) and highest in the middle sample (late July; Figure 2). A quadratic correlation was observed between abundance and grazing intensity in the early and middle samples (early sample: df = 14, F = 8.47, P = 0.0051; middle sample: df = 13, F = 3.79, P = 0.0561). The abundance of C. albonemus was greatest in the moderate and heavy grazing plots in the early sample, and was greatest in the moderate grazing plots in the middle sample. In contrast, a negative linear relationship was observed between C. albonemus abundance and grazing intensity in the late sample (middle August; df = 14, F = 7.78, P = 0.0153). These results demonstrate that heavy grazing and overgrazing decreased the C. albonemus population size.

Sequencing of the C. albonemus adult transcriptome yielded more than 42,359,852 clean reads from the 43,962,497 raw reads, and a total of 154,164,934 nucleotides (transcripts; Supplementary Tables 1, 2). A set of 190,722 transcripts and 132,710 unique sequences were generated, with N50 values of 1,548 and 1,286, respectively (Supplementary Tables 2, 3; Supplementary Figure 1). As expected, half of the sequences annotated in the NCBI non-redundant protein sequence database matched those of insect species, including Tribolium castaneum (15.4%), Acyrthosiphon pisum (6.9%), and Pediculus humanus (6.5%; Supplementary Figure 2). Of the 39,090 (29.45%) unique sequences that were annotated through BLAST searches against the seven indicated databases (NR, NT, SWISS-PROT, PFAM, KEGG, KOG, and GO), 28,999 (21.85%) were annotated using the Gene Ontology (GO) database, and 4,530 were annotated using the nucleotide database (NT; Supplementary Table 4). In the EuKaryotic Orthologous Groups (KOG) database, 13,712 annotated genes were assigned to 26 groups; most of the annotated genes were categorized into the following groups: R, general functional prediction only (4,602 genes); T, signal transduction mechanisms (1,606 genes); and O, posttranslational modification, protein turnover, and chaperones (1,081 genes; Supplementary Table 5). The Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway database was used to annotate 8,600 genes, most of which were categorized as being involved in signal transduction (939 genes), translation (650 genes), or carbohydrate metabolism (550 genes; Figure 3A). The transcriptome data of C. albonemus was submitted to the NCBI Sequence Read Archive (SRA) database (ID: SRP058368).

Transcript levels of C. albonemus differed according to grazing intensity, as assessed by fragments per kilobase per million (FPKM) mapped reads. The FPKM distributions of the five grazing intensity plots were similar, with the histogram showing a peak on either side of zero. The highest transcript level was observed in the heavy grazing plots (Supplementary Figure 3). The highest number of differentially expressed genes was observed in the moderate grazing plots (Figure 3B), and the lowest number in the heavy grazing and overgrazing plots. Results of gene cluster analyses showed similar gene expression patterns in grasshoppers collected in the control plots (no grazing) and moderate grazing plots, with greater numbers of genes identified (Figure 4A). Gene expression patterns were similar between the heavy grazing and overgrazing plots, and more downregulated genes were observed in C. albonemus collected in these plots. Analysis of genes involved in stress response associated with grazing pressure and grasshopper population dynamics showed the upregulation of certain genes with increased grazing intensity (Figures 4B,C), such as polyprotein-like protein, ankyrin repeat domain-containing protein, eupolytin, and late trypsin (Supplementary Data 1). Genes that were downregulated with increased grazing intensity (Figures 4D,E) included yellow-g protein, cadherin-related tumor suppressor-like protein, ribosomal protein S5, and cuticle protein (Supplementary Data 1). To identify the biological functions of these gene products, the differentially expressed genes were mapped to the GO and KEGG databases. More than 20 pathways were enriched by grazing treatment (P < 0.05) including galactose metabolism, lysosome, starch and sucrose metabolism, and other glycan degradation (Table 1).

A total of 28,999 unigenes of C. albonemus were subcategorized into 50 GO classes (Supplementary Table 6). Most of these transcripts were assigned to biological processes (47.65%), cellular components (29.89%), or molecular function (22.46%). In the biological processes category, many of the transcripts appear to be involved in cellular processes (16,580 genes, 22.04%). In the cellular components category, many of the transcripts appear to be involved in cell biology (8,967 genes, 19.00%). In the molecular function category, many of the transcripts appear to be involved in binding (16,593 genes, 46.81%; Supplementary Table 6). Compared with gene expression in the control plots (no grazing), the highest number of upregulated genes were observed in the light grazing plots. They included genes involved in serine-type, peptidase activity, serine hydrolase activity, peptidase activity acting on L-amino acid peptides, and hydrolase activity (Supplementary Table 7). Genes upregulated in the moderate grazing plots were involved in chitin binding, chitin metabolic processes, glucosamine-containing compound metabolic processes, carbohydrate derivative binding, and amino sugar metabolic processes. Genes upregulated in the heavy grazing plots included those involved in serine-type endopeptidase activity, serine-type peptidase activity, serine hydrolase activity, multicellular organism reproduction, and lipid transporter activity. Genes upregulated in the overgrazing plots included those involved in structural constituents of eye lens, single-organism developmental processes, anatomical structure development, sensory organ development, and developmental process metabolism. Compared with gene expression in the control plots, downregulated genes in the moderate grazing plots belonged to only two categories: structural constituents of cuticle and structural molecule activity (Supplementary Table 7). Downregulated genes in other grazing treatments were involved in structural constituents of ribosome, ribosome biogenesis, ribonucleoprotein complex biogenesis, translation, ribosome, chitin binding, and chitin metabolic processes (Supplementary Table 7).

KEGG pathway analysis provided insight into the transcriptional responses to grazing intensity. Upregulated genes in the light grazing plot were primarily involved in galactose metabolism, lysosome, starch, and sucrose metabolism, and other glycan degradation (Table 1). Only three genes involved in galactose metabolism (E1.1.99.1, betA, CHDH) were significantly upregulated in the moderate grazing intensity plot. Genes encoding proteins that are processed in the endoplasmic reticulum (e.g., MAN2B1, HSP) were significantly enriched in the high grazing and overgrazing plots (Table 1). Most of the downregulated genes in plots other than the moderate grazing plots were associated with ribosomes; they included RP-S26e, RPS26, RP-L15e, RPL15, RP-L23e, RPL23, RP-S23e, and RPS23, RP-L12e (q < 0.05; Table 1). Downregulated genes in the moderate grazing plots included LCT and malZ, which are associated with galactose metabolism (q < 0.05). Other downregulated genes in heavy grazing plots included those involved in lysosomes, starch and sucrose metabolism, and amino sugar and nucleotide sugar metabolism.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The handling Editor declared a shared affiliation, though no other collaboration, with several of the authors XQ, XH, JM, ZZ.