Potential Pathways and Genes Involved in Lac Synthesis and Secretion in Kerria chinensis (Hemiptera: Kerriidae) Based on Transcriptomic Analyses

Lac is a type of natural resin secreted by lac insects and is widely used in the military and other industries because of its excellent adhesion and insulation properties. The main ingredients of lac are lactones and lactides, which are formed from hydroxy fatty acids and sesquiterpene esters. In this study, we measured lac secretion rates by the insect Kerria chinensis at different developmental stages and identified lac secretion-minimum and lac secretion-active stages of the insect. We then analyzed transcriptomes of lac secretion-minimum and lac secretion-active stages of the insect. Based on expression profiles of genes in different stages of the insect, we identified pathways and genes that are potentially involved in lac synthesis and secretion in K. chinensis. Our study lays a foundation for future studies to reveal the molecular mechanisms and pathways of lac synthesis and secretion in this beneficial insect.


Introduction
Lac is a mixture of natural resins produced in and secreted by insects from the genius Kerria (Hemiptera: Kerriidae) [1]. The use of lac by humans can be dated back to centuries ago. Initially, lac was mainly used as wood finishing and as an adhesive for ceramics. At present, lac is used in the production of music instruments, insulation of electronics, coating of pharmaceutical pills, painting, material conservation, and various applications in the military industry [1]. Due to its strong adhesiveness, good insulation, moisture resistance, smooth film texture, non-toxic properties, and tastelessness, lac is expected to be continuously widely used in various industries in the future [2,3].
The composition of lac derived from different insect species varies. Host trees and environmental conditions can also affect the chemical composition of lac secreted even from the same insect species. In general, lac is mainly composed of lac resins, shellac waxes, laccaic acids, and other organic chemicals [2]. Lac resins are mixtures of lactones and lactides, which are formed from hydroxyl fatty acids and polyterpene esters [2].
The biosynthetic pathway for polyterpene esters is the mevalonate pathway (MVA) in the cytoplasm [4,5]. In the mevalonate pathway, acetyl-CoA is converted into mevalonic acids by the action difference between crude lac secretion on the initial collection date and crude lac secretion on the final collection date by the number of days. Three different stage female samples were measured, and three independent biological replicates were conducted for each measurement.
The average amount of lac secreted by individual insects was calculated by Formula (1), whereas the individual secretion rate was calculated by Formula (2).
In Formula (1), M represents the amount of individual secretion of the lac insects; m 1 represents the total weight of lac and insects in a lac block (containing lac and insects); m 2 represents the weight of insects in the lac block; and N represents the number of insects in the lac block.
In the Formula (2), v represents the individual secretion rate; M n represents the amount of individual secretion in the nth time; M (n+1) represents the amount of individual secretion in the (n+ 1)th time. T n represents the nth acquisition time; and T (n+1) represents the (n+1)th acquisition time.

RNA Extraction and Sequencing
Total RNA from L, A1, and A2 stages was extracted using Trizol Reagent (Sangon Biotech, Shanghai, China) according to the method provided by the manufacturer. RNA quality was monitored on 1% agarose gels. RNA quantification and integrity were checked using a NanoPhotometer ® spectrophotometer (IMPLEN, Westlake, CA, USA), and a Qubit ® RNA Assay Kit (Sangon Biotech, Shanghai, China) in a Qubit ® 2.0 Flurometer (Thermo Fisher Scientific, Waltham, MA, USA), and an RNA Nano 6000 Assay Kit (Sangon Biotech, Shanghai, China) on an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), respectively. A total amount of 3 µg RNA per sample was used for mRNA enrichment using magnetic beads with Oligo (dT). mRNA species were fragmented into short pieces with average size 150~200 bp. The fragmented mRNA pieces were used as templates for reverse-transcription with random hexamer primers and M-MuLV Reverse Transcriptase (RNaseH)to synthesize the first strand cDNA. Second strand cDNA was synthesized subsequently using DNA Polymerase I (Sangon Biotech, Shanghai, China) and RNase H (Sangon Biotech, Shanghai, China). Remaining overhangs were converted into blunt ends via exonuclease/polymerase of polymerase I (Sangon Biotech, Shanghai, China). After adenylation of 3 ends of the DNA fragments, AMPure XP beads were used for fragment size selection. cDNA libraries were obtained by PCR amplification, and the PCR products were purified by AMPure XP beads. Quality of the libraries was measured by Qubit2.0, Agilent2100, and q-PCR. After cluster generation, the libraries were sequenced on an Illumina Hiseq TM 2500 platform, and paired-end reads were generated.

Data Analysis
After sequencing, raw reads were subjected to credibility analysis. Sequences with a linker, sequences with N greater than 10%, and low-quality sequences (the number of bases with a mass value of Qphrd ≤ 5 more than 50% of the entire reads) were removed and discarded, and the remaining reads were referred as Cleans Reads. Q20, Q30, GC-content, and sequence duplication levels of the clean reads were calculated. All the downstream analyses were based on clean reads data. Transcriptome assembly was accomplished based on the de novo splicing of valid sequences [18]. The transcript sequence information was stored in FASTA format as a reference sequence for subsequent analysis using Trinity.

De Novo Assembly and Gene Annotation
Assembled unigenes were annotated by blasting against the NCBI non-redundant protein (Nr) database, NCBI nucleotide database (Nt), Swiss-Prot protein database (Swiss-Prot), Protein families database (PFAM), Kyoto Encyclopedia of Genes and Genomes (KEGG) database, Clusters of Orthologous Groups (COGs) database, and Gene Ontology (GO) database [19,20]. The e-value thresholds for Nr, Nt, and Swiss-Prot databases were 1e-5 and for KOG/COG database 1e-3. GO functional annotation was done using blast2go (b2g4pipe_v2.5). Annotation against the KEGG database was carried out using the KEGG Automatic Annotation Server (KAAS) with default parameters [21].

Quantitative PCR (qPCR) Analysis
To validate the results of RNA-seq analysis, several genes were selected from K. chinensis transcriptome data to perform qPCR experiments. The list of gene-specific primers is shown in Table S1. Each reaction was carried out with 1 µL of a 1/20 (v/v) dilution of the first cDNA strand, 0.5 mM of each primer in a total volume of 25 µL. Amplification was performed with iQ SYBR Green Supermix (Sangon Biotech, Shanghai, China) on an iCycler real time detection system (Bio-Rad, Hercules, CA, USA) as the following thermal circles: 95 • C for 30 s followed by 40 cycles of 95 • C for 5 s and 60 • C for 30 s. The unigene expression levels were calculated with the 2 −∆∆CT method [22]. Actin was selected as a reference for normalization of template concentration. Three independent biological replicates were carried out for each treatment.

Analysis of Differential Expression Gene (DEG)
Read counts were used to estimate the levels of gene expression. In order to reduce false positives, trimmed mean of M values (TMM) was used to standardize data processing. Differentially expressed genes were identified by DEGseq with q value < 0.005 and log2 (fold change) > 1 [23,24].

Data Availability
We have deposited our data in the Sequence Read Archive (SRA) (http://www.ncbi.nlm.nih.gov/ sra/); the accession for our submission is PRJNA489372.

The Characteristics and Dynamic Regularity of K. chinensis
Male and female larvae can only secrete a small amount of lac. The female is able to secret vast amounts of lac in the adult stage, whereas the male cannot secret lac any more both in the pupal and adult stages [1]. The lac secretion of K. chinensis in the whole summer generation was about 24.9 mg/female ( Figure S2). Through the determination of individual gum secretion rate in each stage, the results showed that the individual gum secretion rate was slower in the larval stage and increased gradually in the adult stage. In the early adult stage (A1), the individual fastest secretion rate was 4.82 × 10 −1 mg/d ( Figure 1B), and the secretion decreased gradually in the mid-late adult stage (A2).

Overview of Transcriptomic Analyses
Based on the characteristics and dynamic regularity of K. chinensis observations, a strategy of transcriptomic analyses was designed to identify candidate genes and pathways that are potentially involved in lac synthesis and secretion ( Figure S3). The strategy was conducted to identify genes that are differentially expressed between lac secretion-minimum stage (L) versus lac secretion-active stages (A1 and A2). Those genes that are commonly up-or down-regulated in A1 and A2 in comparison with L are potentially involved in lac synthesis and secretion. Those genes that are differentially regulated between A1 and A2 in comparison with L are likely involved in other developmental regulations.
Illumina sequencing of cDNA libraries generated a total of 24,897,365 raw reads covering a total of 126,059,825 nucleotides. After removing adapter sequences, ambiguous nucleotides and low-quality sequences, 24,455,197 high quality reads were retained. The mean length of reads was 486 nt. The high quality reads were assembled into 183,899 unigenes, which had average size of 532 nt and N50 670 nt. Detailed statistics of spliced transcripts and assembled unigenes are shown in Table S2 and Figure S4. These statistics were comparable with those of high standard transcriptomes from other well known insect species.
To annotate the unigenes, a blastx search was carried out against the NR database with a cutoff E-value smaller than 10 −5 (Table 1). A total of 54,037 unigenes (29.38% of all unigenes) returned above cutoff alignments when searched against the NR nucleotide database. Specifically, 71.6% of matched sequence alignments had E-values ranging between 1.0E −5 and 1.0E −60 , and the remaining alignments were with E-values less than 1.0E −60 ( Figure S5). The similarity of 7.2% sequence alignments ranged from 96% to 100%, 24.0% ranged from 80% to 95%, 46.9% ranged from 60% to 80%, 21.8% ranged from 40% to 60%, and only 0.1% ranged from 18% to 40% ( Figure S5). The species distribution of the first hits was 10% with Acyrthosiphon pisum, 9.3% with Zootermopsis nevadensis, 7.2% with Nasonia vitripennis, 4.2% with Diaphorina citri, and 3.7% with Bombyx mori ( Figure S5). In addition, the unigenes were similarly analyzed against other databases, and the results are summarized in Table 1 as well. Of the 183,899 unigenes, 50,406 were assigned to specific GO terms and classified into 19 subcategories of 'biological processes', 15 subcategories of 'cellular components', and 12 subcategories of 'molecular functions' ( Figure S6). The subcategories of 'cellular process', 'cell', and 'binding' were dominant in each of the three major categories, respectively. Analysis of the unigenes against the KEGG database, a network diagram of cell metabolic pathways, assigned 30,464 unigenes specific pathways ( Figure S7). Among them, 15,460 unigenes were annotated to 'metabolism' (40.07% of the total), 3839 to 'environmental information processing' (9.95%), 7929 to 'genetic information processing' (20.55%), 7162 to 'organismal systems' (18.56%), and 4011 to 'cellular processes' (10.47%).

Genes Expressed Differentially between Insects at Different Stages
There were 2774 genes expressed differentially between A1 and A2, with 990 genes up-regulated and 1284 down-regulated. There were 2331 genes expressed differentially between A2 and L, with 706 genes up-regulated and 1625 down-regulated. There were 1330 genes expressed differentially between A1 vs. L, with 937 genes up-regulated and 393 down-regulated (Figure 2A). Combinational analyses revealed that 202 genes were commonly up-regulated and 1187 genes commonly down-regulated in the two lac secretion-active adult stages in comparison with the lac secretion-minimum larval stage ( Figure 2B). The remaining genes were differentially regulated among the three stages. Since the commonly up-and down-regulated genes are most likely involved in lac synthesis and secretion, the composition of these genes were analyzed in detail in the following sections.

Reduced Protein Synthesis in Lac Secretion-Active Stages
As shown in Figure 2C,D, the largest proportions of genes that were both commonly up-and down-regulated are un-annotated. These unique genes represent opportunities to identify novel genes with new functions. The roles of these un-annotated genes in lac synthesis and secretion remain to be studied.
The second largest category of genes that were commonly down-regulated was 'protein metabolism', which included proteins involved in protein synthesis, folding, modification, and degradation ( Table 2, Table S3). There were 21% of commonly down-regulated genes belong to the 'protein metabolism' category. On the other hand, there were only 8% of commonly up-regulated genes belong to this category. Most of the commonly down-regulated genes encode various ribosomal proteins, followed by chaperones, and proteases or protein degradation-related, indicating that protein synthesis was at much lower levels in lac secretion-active stages of the insect. In comparison, most of the commonly up-regulated 'protein metabolism' genes encode proteins involved in folding and degradation (Table S4). The combination of down-regulation of genes in protein synthesis and up-regulation in protein degradation indicated an overall reduction in protein synthesis in lac secretion-active stages. This could be a way for the insect to save energy and nutrients for the production and secretion of lac.

Reduced and Unbalanced Central Metabolic Pathways
Overall, metabolism via central metabolic pathways was reduced in insects at lac secretion-active stages based on the proportion of genes that are commonly up-or down-regulated (Table 2, Figure S10). There were 70 down-regulated genes involved in central metabolic pathways and respiratory chain with average of 5.59 log2 fold down-regulation. In comparison, there were only 22 genes up-regulated with only 1.25 log2 fold increase. However, not every step in the three central metabolic pathways was affected evenly. All (37) identified genes encoding components in the respiratory chains were down-regulated with similar fold changes, indicating that overall ATP production for energy use was lowered. However, there were genes for some enzymatic steps, which were actually up-regulated (Table 2, Figure S10). This uneven regulation of central metabolic pathways suggested that these pathways were regulated for accumulation of certain intermediates for biosynthesis. Presumable accumulated intermediates include glucose-6-phosphate, ribulose-5-phophate, glyceraldehide-3-phosphate, sedoheptulose-7-phosphate, glucose-6-phosphate, fructose-1,6-bisphosphate, citrate, succinyl-CoA, and malate. These are activated sugars and other intermediates that could be used for lac synthesis directly or indirectly.

Shifted Metabolism of Fatty Acids, Lipids, and Terpenes
In contrast to the overall reduction in protein metabolism and unevenly affected central metabolic pathways, the metabolic pathways of fatty acids, lipids, and terpenes were shifted in the lac secretion-active stages in comparison with the lac secretion-minimum stage ( Table 2). There were 24 genes commonly up-regulated and about an equal number (23) of genes down-regulated, even though the magnitude of down-regulation was much bigger than that of up-regulation. The up-and down-regulation of genes within the same functional categories indicated a shift in formation of different products of the same categories at different stages. We hypothesize that the genes up-regulated in the lac secretion-active stages are involved in lac production and secretion.

Candidate Genes Potentially Involved in Lac Synthesis and Secretion
Based on the correlation between gene expression levels and lac accumulation rates as well as the putative functions of genes, candidate genes potentially involved in lac synthesis and secretion were identified, including genes involved in fatty acid synthesis: As shown in Figure 2E, Table S5, the expression of these genes was well correlated with the secretion of lac. The genes putatively involved in fatty acid synthesis may produce hydroxyl fatty acids contained in lac. The genes putatively involved in terpenoid synthesis may participate in production of sesquiterpene acids, which are major components of lac [1]. The involvement of UDP-dependent glycosylation in lac synthesis remains very speculative. The reason we think it might be involved directly or indirectly is the apparent accumulation of activated sugars based on unbalanced central metabolic pathways. The same types of genes are specifically up-regulated as well in wax secretion insects [15]. Among these differentially expressed genes, 14

Discussion
A total of 205,581 transcripts were obtained in this study, among which 183,899 unigenes were successfully annotated by seven databases including Nt, Nr, Swiss-Prot, COG, GO, KOG, and KEGG. The most sequences were annotated in A. pisum by homologous sequence alignment in the Nr database; they may have a relatively close genetic relationship and the A. pisum has abundant genomic information [25]. There are 50,406, 36,724, and 30,464 unigenes annotated to GO, COG, and KEGG databases, respectively, which provide rich data resources for deeper studies on the mechanism of lac secretion of lac insects in the future.
The secretion rate of lac was different between different individuals. The secretion rate of Kerria lacca gradually decreased in the adult stage; the highest increase of secretion rate of each female was 2.54 × 10 −4 mg/d in the early adult, which gradually decreased at middle and late stages [26]. At larval and late adult stages of Kerria yunnanensis, the individual secretion rate was relatively slow. In the middle adult stage, the maximum secretion rate was 1.46 × 10 −1 mg/d, which gradually decreased at late stages [27]. Through the study on the secretion rate of K. chinensis, it was found that the mid-late stage had a maximum secretion rate of 4.82 × 10 −1 mg/d, and the secretion rate gradually decreased at later stages. Therefore, the study proved that the expression of genes associated with secretion of lac in the adult stage was significantly different (especially higher than other stages), which will provide a foundation for screening the follow-up genes associated with lac biosynthesis.
Based on our results, we propose a hypothesis for lac synthesis in K. chinensis ( Figure S11). The model consists of three major physiological changes that provide the basis for lac synthesis, including savings of nutrients and energy, production of intermediates and NADPH necessary for lac biosynthesis, and a shift in metabolism of fatty acids, lipids, and terpenes that results in reduction of fatty acids, lipids, and terpenes for other physiological uses but increase of these compounds for lac synthesis. To save nutrients and energy, the overall metabolism is lowered in lac secretion-active stages, particularly in protein synthesis and structural components. The central metabolic pathways, including pentose phosphate pathway, glycolysis, and TCA cycle, were unequally affected, leading to accumulation of activated simple sugars and other intermediates for lac synthesis. The shift in metabolism of fatty acids, lipids, and terpenes results in reduced production of these compounds in other part of tissues that are used for normal growth and development, but elevated production in the lac secretion glands for the production of lac. Our model hypothesis provides directions for further experimental tests in the future regarding lac synthesis in the insect K. chinensis.

Conclusions
In this paper, the amount of individual secretion and rate of Kerria chinensis as the object of study of lac secretion in different periods were measured to reveal the lac secretion characteristics and dynamics. By analyzing the transcriptomes of different stages of K. chinensis, we speculated the correlation of gene functions through the link of gene expression level and lac secretion rate; the results showed that 28 candidate genes that may be involved in lac synthesis were preliminarily screened; meanwhile, the results provided a scientific basis for the verification of molecular mechanisms related to lac secretion.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2075-4450/10/12/430/s1. Figure S1: The pathways of terpenoid and fatty acid biosynthesis. (A) The procession of terpenoid biosynthesis. ATOT (acetoacetyl-CoA thiolase); HMGS (hydroxymethylglutaryl-CoA synthase); HMGR (hydrox-methylglutaryl-CoA reductase); MK (mevalonate kinase); PMK (phosphomevalonate kinase); MPD (diphos phomevalonatedecarboxylase); IDI (isopentenyl diphosphate isomerase); FPPS (farnesyl diphosphate synthase); TPS (terpene synthase). (B) The procession of Fatty acid biosynthesis, ACC (acetyl-coenzyme A carboxylase); FAS (fatty acid synthase); ELO (elongase of very long chain fatty acid); FAD (fatty acid desaturase), Figure S2: Lac secretion of females at different stages, Figure S3: The strategy to identify shellac-secretion related genes, Figure S4: Length distribution of unigenes of K. chinensis, Figure S5 Figure S6: GO classification analysis of unigenes in unigene. GO functions is showed in Y-axis. The X-axis shows the number of genes that have the GO function, Figure S7: KEGG classification analysis of unigenes in all-unigene. KEGG functions is showed in Y-axis. The X-axis shows the percentage of genes that have the KEGG function, Figure S8: Metabolism of terpenoids, polyketides, and lipid metabolism pathway. (A) Metabolism of terpenoids and polyketides pathway. (B) Lipid metabolism pathway. Metabolism pathway is shown in Y-axis. The X-axis shows the number of genes, Figure S9: COG function classification of unigenes. The horizontal coordinates are function classes of COG and the vertical coordinates are numbers of unigenes in one class. The notation on the right is the full name of the functions in X-axis, Figure S10: Central metabolic pathways. (A) Pentose phosphate pathways. (B) Glycolysis. (C) Citruate acid cycle. Blue numbers: the number of genes related to the pathways of pentose phosphate pathways, Glycolysis and citruate acid cycle enzymes in transcriptome data, Figure S11: The lac synthesis in K. chinensis; Table S1: The quantitative PCR primers of 14 candidate genes in K. chinensis, Table S2: Assembly quality of K. chinensis transcriptome, Table S3: The commonly down-regulated genes metabolic pathways, Table S4: The commonly up-regulated genes metabolic pathways, Table S5: FPKM of 28 putative genes related to fatty acid synthesis, terpenoid biosynthesis, and UDP-driving glycosylation.
Author Contributions: H.C., X.C., and M.-S.C. designed and led the project; W.W. collected samples and afforded photos of fresh collection; W.W., J.Z., and X.L. designed the figures; W.W., P.L., and Q.L. played major roles in the writing of the manuscript and its publication. All Authors reviewed the manuscript, which is not being considered for publication, in whole or in part, in another journal.