Characterization of Adelphocoris suturalis (Hemiptera: Miridae) Transcriptome from Different Developmental Stages

Adelphocoris suturalis is one of the most serious pest insects of Bt cotton in China, however its molecular genetics, biochemistry and physiology are poorly understood. We used high throughput sequencing platform to perform de novo transcriptome assembly and gene expression analyses across different developmental stages (eggs, 2nd and 5th instar nymphs, female and male adults). We obtained 20 GB of clean data and revealed 88,614 unigenes, including 23,830 clusters and 64,784 singletons. These unigene sequences were annotated and classified by Gene Ontology, Clusters of Orthologous Groups, and Kyoto Encyclopedia of Genes and Genomes databases. A large number of differentially expressed genes were discovered through pairwise comparisons between these developmental stages. Gene expression profiles were dramatically different between life stage transitions, with some of these most differentially expressed genes being associated with sex difference, metabolism and development. Quantitative real-time PCR results confirm deep-sequencing findings based on relative expression levels of nine randomly selected genes. Furthermore, over 791,390 single nucleotide polymorphisms and 2,682 potential simple sequence repeats were identified. Our study provided comprehensive transcriptional gene expression information for A. suturalis that will form the basis to better understanding of development pathways, hormone biosynthesis, sex differences and wing formation in mirid bugs.


Results
Illumina sequencing and reads assembly. To acquire an in-depth understanding of molecular mechanisms controlling A. suturalis biology at various developmental life stages, cDNA was generated from eggs (clean reads accession number: SRX533456), 2 nd instar nymphs (clean reads accession number: SRX534296), 5 th instar nymphs (clean reads accession number: SRX534293), A. suturalis female adults (AFA) (clean reads accession number: SRX534295), and A. suturalis male adults (AMA) (clean reads accession number: SRX534294), followed by sequencing using the Illumina HiSeq TM 2000 sequencing platform. A total of 264,212,834 clean reads were generated and assembled into 677,305 contigs from the five cDNA libraries (Table 1 and Table 2). The egg cDNA library produced the most data (55,038,084  clean reads), while the AFA cDNA library produced the least (51,295,410 clean reads). Q20 score (the average quality value) was > 97% (Table 1), and SOAPdenovo was used to map raw reads 17 . Quality checks and contig assembly were carried out using the Trinity software < http://trinityrnaseq.sourceforge. net/> 18 after removal of adaptor sequences. The assembled reads representing these different developmental stages resulted in 88,614 unigenes and a total length of 70,479,040 base pair (bp), with a mean contig length of 795 bp (N50 = 1,398 bp), including 23,830 distinct clusters and 64,784 distinct singletons ( Table 2). The unigene sequence length size distributions showed that 22,230 unigenes were over 1,000 bp. Across the entire developmental stages of A. suturalis, the unigene sequence length size distribution indicated that 11,993 unigenes of eggs, 11,344 unigenes of 2 nd instar nymphs, 10,692 unigenes of 5 th instar nymphs, 10,352 unigenes of AFA, and 13,223 unigenes of AMA were more than 1,000 bp (Table S1).

GO, COG and KEGG
To better understand the biological pathways involved in A. suturalis, potential pathways were investigated using all 19,204 sequences identified in the KEGG Orthology (KO) annotation of the mirid bug, with sequences mapped to the reference authoritative pathways in the KEGG database. As a result, a total of 258 KEGG pathways were identified, with 2,971 unigenes (15.47%) being involved in metabolic pathways. These pathways played a dominant role in pathways associated with RNA transport, olfactory transduction, insect hormone biosynthesis, and ubiquitin mediated proteolysis, amongst others (Table  S2).
suturalis unigene dataset to determine the nature and frequency of SSRs within our reduced representation of its genome. Among the total of 3,293 unigenes mapped, more than 3,200 pairs of potential SSR primers were designed, although primer efficacies were not tested with respect to A. suturalis population and evolutionary genetics, nor their relationships to microsatellite DNA family 20 , transposable elements 21 , and implications to genome organization and rates of evolution 22 . By searching for di-, tri-, tetra-, penta-and hexa-nucleotide repeats, we obtained a total of 2,682 SSRs, of which the majority of SSRs have less than 24 repeat units, with 5-8 repeat units being the most common categories. Among all SSRs identified, dinucleotide repeats (47.84%) represented the most abundant microsatellite repeat units, followed by trinucleotide (45.38%), tetranucleotide (3.65%), pentanucleotide (2.65%), and hexanucleotide (0.48%) repeats, respectively (Fig. 5).
SNPs identification. 791,390 potential SNPs were identified using the program PolyBayes (Table 4)  samples recorded the least (39,225 A-G; 38,365 C-T). 40.53% of SNPs were transversions, with 14.78% being A-T transversions, followed by G-T (9.83%), A-C (9.77%), and C-G (6.14%) transversions. From the entire dataset, the average SNP frequency is 0.011 SNP/nt and the average transition and transversion frequency is 0.0067 transition/nt and 0.0046 transversion/nt, respectively. The complete list of SNPs identified across the selected developmental stages (i.e., eggs; 2 nd and 5 th nymph stages; AMA; and AFA) is provided in Table S3.
Comparison of gene expression profile among A. suturalis developmental stages. To confirm the relative expression levels of differentially expressed genes across the major developmental stages, we calculated the number of clean tags for every gene, and gene expressed differentially were identified using the IDEG6 tool 23 between every two A. suturalis life-stages that represented selected transitional developmental stages (eggs and 2 nd instar nymphs, 2 nd and 5 th instar nymphs, 5 th instar nymphs and adult stages (i.e., AFA, AMA), the reproductive stage (i.e., AFA and eggs), and sex differences (i.e., AMA and AFA)) ( Fig. 6). Furthermore, we also analyzed the change of gene expression between different developmental stages by GO and KEGG pathway enrichment analysis.
In the comparison between eggs and 2 nd instar nymphs, the expression profiles of 34,930 genes had changed. There were 19,147 unigenes that were up-regulated in the 2 nd instar nymph library and 15,783 unigenes down-regulated in the egg library (Fig. 6). Among the top ten up-regulated genes, significant matches included the gene homologous to one that encodes the odorant-binding protein 2 in the related species A. lineolatus; a cysteine proteinase gene in Dictyostelium fasciculatum, three cuticular    protein genes (Bombyx mori putative cuticle protein; C. floridanus cuticle protein 21, and C. floridanus cuticle protein 19.8), and three predicted functional genes (T. castaneum hypothetical protein TcasGA2_ TC007306; Columba livia putative serine/threonine-protein kinase kinX, and Macaca mulatta type I inositol-1,4,5-trisphosphate 5-phosphatase-like). The top ten down-regulated genes in eggs libraries included predicted functional genes; an uncharacterized protein gene LOC101240627 (Hydra magnipapillata); Proclotting enzyme gene (C. floridanus), GK14571 (Drosophila willistoni), Toutatis gene (T. castaneum), one Longitudinals lacking protein isoforms A/B/D/L (Acromyrmex echinatior), the hypothetical protein TcasGA2_TC006306 gene (T. castaneum)); the GF17401 gene (D. ananassae), and one transport protein Sec24C (H. saltator) (Table S4). Differentially expressed genes between eggs and 2 nd instar nymphs were characterized into three groups from the GO classification: (i) Biological process, (ii) Cellular component, and (iii) Molecular function. The results showed 4,118 biological processes-related genes were predominantly concentrated in the 'Cellular process' subcategory, and 3,329 genes involved in the 'single-organism process' subcategory. A total of 3,067 genes were associated with the 'cell/cell part' subcategory within the 'Cellular component' category, while 3,038 genes were involved in 'binding' within the 'Molecular function' category (Table S5a). In the pathway analysis (Table S5b), among all 19,204 unigenes involved in 257 pathways, the most differentially expressed genes from the eggs and 2 nd instar nymphs comparison were involved in 'metabolic pathways' (1,302), 'regulation of actin cytoskeleton' (444), 'RNA transport' (394) and 'focal adhesion' (376).
Between 2 nd and 5 th instar nymph libraries, 18,360 genes showed significant expression changes from the comparative analysis, with 8,391 genes up-regulated in the 5 th instar nymph library and 9,969 down-regulated genes in 2 nd instar nymph library (Fig. 6). Among the top ten differentially up-regulated and down-regulated expressed genes, 4 showed defined functions, i.e., a cytochrome c oxidase subunit 6B1 gene (A. echinatior) in up-regulated genes group; a 48 kDa subunit-like protein (M. rotundata), a Translocon-associated protein subunit delta precursor gene, and a sugar transporter (Culex quinquefasciatus) in the down-regulated genes group. Two of ten up-regulated genes have predicted functions, i.e., TcasGA2_TC005836 protein gene (T. castaneum), and CBG23181 protein gene (Caenorhabditis briggsae). Among ten most down-regulated genes, five genes have predicted functions, i.e., putative transporter SVOPL-like gene (Bombus impatiens), dolichyl-diphosphooligosaccharide--protein glycosyltransferase 48 kDa subunit-like (M. rotundata), hypothetical protein SINV_03497 (Solenopsis invicta), RNA-binding protein (Triatoma infestans), and regulator of G-protein signaling 12-like (A. pisum) (Table S4). Gene sets in this library were predicted by GO analysis (Table S6a) Table S6b).
The comparative analysis between 5 th instar and AFA libraries revealed 23,528 genes with significant expression profile changes (Fig. 6). Two of the ten most up-regulated genes in the 5 th instar library have homologies to D. ananassae GF17401 and GH24370 genes. Three genes have predicted functions, i.e., cathepsin B-like of Oryzias latipes and cathepsin L-like of Strongylocentrotus purpuratus, and opacity protein and related surface antigens of Ixodes scapularis. Three genes have defined functions: RP45 gene of Rhodnius prolixus, CBN-CPR-4 protein gene of C. brenneri, reproduction-associated glycoprotein 7 of Odocoileus virginianus, and two unannotated genes. Among the ten most down-regulated genes, one matched the reverse transcriptase gene of Aedes aegypti and one matched D. melanogaster's sallimus, isoform Z. Three genes have predicted functions, i.e., general transcription factor II-I repeat domain-containing protein 2A-like (Xenopus (Silurana) tropicalis), RNA-directed DNA polymerase from mobile element jockey-like (Metaseiulus occidentalis), and high affinity cAMP-specific and IBMX-insensitive phosphodiesterase (A. pisum) (Table S4). GO analysis (Table S7a) revealed that gene sets in this library showed high accordance with genes relating to 'Cellular process' (2,173), 'single-organism process' (1,772), 'metabolic process' (1,664), and 'biological regulation' (1,010) within the 'Biological process' category. Among the 'Cellular component' category, 'cell/cell part' (1,623), 'organelle' (1,215) and 'organelle part' (719) were the three subcategories with the highest number of gene sets. Among the 'Molecular function' category, 1,748 genes acted as catalytic activity, and 1,607 genes were involved in binding. Among 256 pathways identified (Table S7b), gene sets were correlated to 'metabolic pathway' (885), 'regulation of actin cytoskeleton' (302), 'RNA transport' (243), 'cell cycle' (131), 'RNA polymerase' (90), and 'RNA degradation' (58).

Quantitative Real-time PCR (qRT-PCR) analysis.
To better validate the sequencing data, various genes related to wing formation (flightin, titin, myosin heavy chain (MHC), laminin), hormone biosynthesis (neurofibromin), sex difference (vitellogenin, thioredoxin, KLHL10 (Kelch-like protein 10), and detoxification (Cytochrome P450) were chosen randomly and quantified by the qRT-PCR method. Results from the qRT-PCR experiment supported those obtained from transcriptome analysis and demonstrated similar tendency in up-or down-regulation of A. suturalis genes (Fig. 7). From analysis of deep sequencing data, some genes related to wing dimorphism had higher copy numbers in the AMA library when compared with the AFA library, such as flightin, MHC and laminin. In the AMA cDNA library (Fig. 7a), among the 45 copies of flightin genes and 34 copies of KLHL10 genes as identified in comparison with the AFA library, our chosen single flightin gene and three KLHL10 genes showed up-regulation that concurred with our differential gene expression (DGE) analysis. Down-regulation of the titin gene as determined from DGE analysis (Fig. 7a) was also validated by qRT-PCR (Fig. 7b).

Discussion
In this study, we report on transcriptomic analyses derived from targeted A. suturalis life developmental stages, thereby significantly increased molecular resources available for this arthropod pest of Bt cotton and food crops, and provides a framework for understanding changes in gene expression during development. This transcriptional information not only serves as a valuable resource to better understand underlying biological and physiological mechanisms governing A. suturalis life cycle, but may also lead to the identification of novel targets for bio-rationally designed strategies for the control of this insect pest.
Analysis from the RNA sequence data, amongst the generated 88,614 unigenes and predicted against the NCBI nr protein database, 29,971 sequences of the A. suturalis transcriptome were annotated and showed specific Swiss-Prot matches. Considering the above results, the lower-bound of genes for the A. suturalis transcriptome should be more than 88,614. In fact, many sequences assembled did not match Scientific RepoRts | 5:11042 | DOi: 10.1038/srep11042 significantly to DNA/protein databases due to their general short sequence length, or that they represented significantly different genes. For example, 13.1% of genes in the red flour beetle T. castaneum, the closest relative of A. suturalis with a sequenced and annotated genome, showed no homology to other metazoan genes 24 . Although lacking the full genome information of A. suturalis, our collection of unique transcripts for different developmental stages possibly represents significant proportion of functional genes from this insect.
An interesting finding from this study is that identification of many differentially expressed genes was achieved through comparative transcriptomic analyses across the pest insect's developmental stages, thereby greatly enriching current knowledge of A. suturalis gene expression profiles. This will contribute to future research relating to molecular characterization of olfactory systems, chemical targets, developmental mechanisms and sex difference of this and other mirid plant bugs. During the early developmental stage of A. suturalis, most of the differentially expressed genes were up-regulated in the egg stage when compared with the 2 nd instar nymph stage. In contrast, most of these genes were found to be down-regulated in the AFA and eggs life stages comparison. When comparing between different developmental libraries (e.g., eggs and 2 nd instar nymphs; AFA and eggs), a large number of genes showed specific life stage-related expression profiles (e.g., odorant-binding protein, cuticular protein, cysteine proteinase, vitellogenin) that were likely involved in developmental differentiations.
From the comparison between AMA and AFA libraries, 8,185 up-regulated genes in the AFA library and 6,584 down-regulated genes in the AMA library were identified. 23 genes encoding catalase-like proteins were found to be up-regulated, while various genes such as the U4/U6 small nuclear ribonucleoprotein Prp31-like protein gene and the zinc finger protein RTS2 gene supporter of activation of yellow protein gene were down-regulated. Results obtained for the gene set enrichments in the study that compared AMA and AFA suggested that A. suturalis female adults likely to possess more active transcriptional and translational processes than did adult males, a phenomenon also found in the brown planthopper Nilaparvata lugens (Stål) 25 . Additionally, 1,351 sex biased genes (e.g. transformer-2 (tra-2) gene, odorant binding protein 6, cuticular protein 58) were investigated by transcriptomic analysis in the whitefly Bemisia tabaci (Hemiptera: Aleyrodidae) 26 , many of these genes were also identified and characterized in our transcriptome data, thereby further enriching understanding of sex differentiation in non-model pest organisms. Our transcriptome data also indicated that some A. suturalis genes were likely to have multiple roles, such as to both be associated with sex difference and participated in mechanisms related to cellular components, such as the masticatory epithelia keratin 2p protein, which was also involved in mechanisms associated with cytochrome P450 enzyme precursors. Interestingly, in the related L. hesperus (Hemiptera: Miridae), homologous protein domains were also found in transcriptome analysis of adults 27 . Our study showed that 191 cytochrome P450 genes in nymphs and 75 in adults had significantly high expression levels; almost ten times higher than the 28 cytochrome P450 genes that showed high relative expression level in pre-adults and adults stages of the hemipteran bug Halyomorpha halys 28 , and 65 times more than the transcriptome analysis of cytochrome P450 genes in the bed bug Cimex lectularius 29 . Combining this and previous studies will significantly contribute to understanding the functional roles of cytochrome P450 gene family in hemipteran insects.
In the hemipteran insect Diaphorina citri, 13 juvenile hormone genes were previously identified by transcriptome sequencing 30 . In this study, however, our transcriptomic data have identified only three hormone related neurofibromin genes, of which we validated the gene expression profiles of only one of these three genes in both AFA and AMA. Other highly expressed genes in A. suturalis adult life stage (e.g., masticatory epithelia keratin formation, zinc finger protein RTS2, supporter of activation of yellow protein) were identified which were potentially relevant to various key life developmental processes (e.g., morphogenesis, sex difference). Validating the roles of these genes in morphogenesis or sex differences will be necessary to confirm and characterize gene functions, and may offer novel gene-targeted pest management strategies for A. suturalis. We have validated expression profiles of selected groups of genes in this study based on qRT-PCR but have not validated their gene functional profiles. Our qRT-PCR results, although limited due to the low number of genes being tested, nevertheless offered support for the fold change as detected by Illumina sequencing. Furthermore, a total of 2,368 genes were found to show differential expression profiles across the whole life cycle of this insect, and these included candidate genes that underpinned homeostasis functionalities or feeding behavior (e.g., heat shock proteins, calcium-binding proteins, salivary secreted proteins) based on DGE analyses, and will require future empirical expression profile studies and gene characterization, with novel molecular gene-knockdown techniques (e.g., the CRSPR method 31 or miRNA technique 32 ) offering promising insights to characterizing gene functions in this insect pest.
This study has generated much needed genetic and genomic resources for the development of targeted control of A. suturalis. We have also identified 2,682 microsatellite DNA/simple sequence repeat units and 791,390 putative SNPs from different life stages of A. suturalis. We showed that the most common SSR types in the reduced genome of A. suturalis being dinucleotide repeats (i.e., GC/TG/AT), while transition SNPs far out numbered transversion SNPs across all life stages studied. As with other biological systems [33][34][35][36] , these potential molecular markers identified in this study will serve as useful tools to facilitate future construction of genetic maps, infer gene flow, identify genetic signatures of selection, and infer parentage studies in this emergence pest of Bt cotton and food crops in China.

Materials and Methods
Insect materials and rearing conditions. A. suturalis adults (male and female) were originally collected from a Bt cotton field at the Henan Research and Experiment Base for Modern Agriculture (35°0'13.2372"N, 113°42'28.8792"E) of Henan Academy of Agricultural Sciences (HAAS), Yuanyang, Henan province, China. Eggs, nymphs and adults were reared on green beans at 26 ± 1 °C, 70 ± 10% relative humidity, and a photoperiod cycle of 14 hr light/10 hr dark according to Feng et al. 8 . All insects of different life stages collected in this study originated from a single A. suturalis (monandrous) mating pair. Due to the limitation in female lifetime egg-laying capacity when reared with artificial diet (range 35-172 eggs/generation 37 ), we obtained sufficient eggs, 2 nd and 5 th instar nymphs, male and female adult individuals through a second-generation population-wide mating between siblings.
Samples collection and RNA isolation. Eggs were collected with damp filter papers within 24 h of oviposition 8 , 350 eggs were used in RNA extraction. The 2 nd and 5 th instar nymphs were distinguished based on developmental time since emergence from eggs. AFA individuals were identified by the black line in the abdomen of female adults. Male and female adults were collected four days after emergence from final (5 th ) nymph stage. 100 2 nd instar nymphs and 70 5 th instar nymphs, and 50 individuals of each female and male adults were collected into a glass tube using an aspirator. RNA isolation from all samples was performed as described in Dong et al. 38 , and the integrity of all RNA samples was ascertained using a 2100 Bioanalyzer (Agilent Technologies) with a minimum RNA integration value of 6. Total RNA was isolated from eggs (OD260/280 = 1.926), 2 nd instar nymphs (OD260/280 = 2.017), 5 th instar nymphs (5th) (OD260/280 = 2.024), female (AFA) (OD260/280 = 2.115) and male (AMA) adults (OD260/280 = 2.022) by the SV Total RNA Isolation System (Promega, USA).
Scientific RepoRts | 5:11042 | DOi: 10.1038/srep11042 cDNA library construction and Illumina sequencing. To acquire the entire gene expression information, Illumina sequencing was performed for all five A. suturalis cDNA libraries as described in Xie et al. 39 that represented the various developmental stages. All cDNA libraries were generated according to the protocol as supplied by Illumina (Illumina Inc., San Diego, CA, USA). Briefly, 20 μ g of pooled total RNA samples of individual A. suturalis developmental stages were digested by DNase I (Sigma, USA), and mRNA purified by oligo(dT) magnetic beads and fragmented into fragments of ~100-400 bp. The cDNA was synthetized by reverse transcriptase (Invitrogen, Carlsbad, CA, USA) using random hexamer-primers with the mRNA fragments as templates. Agilent 2100 Bioanaylzer and ABI StepOnePlus Real-Time PCR System were used to quantify and qualify all libraries. We selected five cDNA libraries with 200 bp insert size for sequencing by the Illumina HiSeq2000 TM (Illumina Inc., San Diego, CA, USA).
Bioinformatics analysis of sequencing results. Short sequence reads were assembled for individual samples from the targeted developmental stages of A. suturalis using the Trinity software. Unigenes from each of developmental stages were combined to create a lifetime unigene database of A. suturalis by the BLAST-Like Alignment Tool < http://genome.ucsc.edu/cgi-bin/hgBlat> . All raw sequence data has been uploaded to the NCBI Sequence Read Archive (SRA; accession number SRP041523).
Unigenes longer than 350 bp were aligned against the NCBI nucleotide database (non-redundant (nr)/ non-redundant nucleotide (nt)), Swiss-Prot, COG, KEGG and TrEMBL databases. Gene ontology (GO) annotations of unigenes were achieved by the Blast2GO software < http://www.geneontology.org> . For coding sequences (CDS) annotation, BLASTx alignments were carried out between unigenes against nr, KEGG, Swiss-Prot and COG protein databases, and the transcriptional directions and coding frame of unigenes were predicted from BLASTx results. Unigene CDS with no specific BLASTx matches were predicted by ESTScan 40 . All the above searches were carried out with a cut-off e-value of < 10 −5 .

Identification of simple sequence repeats (SSRs) and single nucleotide polymorphisms (SNPs).
To assess the assembled unigene quality and to identify potential new molecular markers, potential SSR markers were identified among all unigenes using the MISA tool < http://pgrc. ipk-gatersleben.de/misa/> for detecting SSR sequences 41 and possible PCR primers were designed using PRIMER v.3 42 . SNPs within unigenes were evaluated using SOAP2 and SOAPsnp 43 across the various A. suturalis developmental stages (i.e., eggs, 2 nd instar nymphs, 5 th instar nymphs, AFA, and AMA), using all unigenes as reference for SNPs discovery. All common and unique SNP sites among these samples were identified as previously described 36 .
Differential gene expression analyses. Fragments Per Kilobase of transcripts per Million mapped reads (FPKM) method 44 was used to obtain the unigene expression levels of different developmental stages. Base on FPKM values, differentially expressed genes were identified using the IDEG6 software < http://telethon.bio.unipd.it/bioinfo/IDEG6/> 23 . Gene expression variations were analyzed for specific life-stages comparisons that included the following three categories: (i) transitional developmental stages (i.e., eggs and 2 nd instar nymphs (eggs and 2 nd ), 2 nd instar nymphs and 5 th instar nymphs (2 nd and 5 th ), 5 th instar nymphs and female adults (5 th and AFA), 5 th instar nymphs and male adults (5 th and AMA)); (ii) sex differences (i.e., female adults and male adults (AFA and AMA)); and (iii) reproductive stage (i.e., female adults and eggs (AFA and eggs)).

Quantitative Real-time PCR (qRT-PCR) validation.
A total of 50 male and 50 female adults were collected four days after emergence from the final (i.e., 5 th ) nymph stage, and their total RNA purified for cDNA libraries preparation as described above. For generating the first-strand cDNA, 2 μ g total RNA of female and male A. suturalis adults was reverse transcribed in a 20 μ l reaction volume using the PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa). To confirm single amplicon product for qRT-PCR and bio-specificity for A. suturalis, gradient PCR was performed for all primer pairs and optimized for a 60 °C annealing temperature for all primer pairs. PCR amplification products were run on 1.5% 1xTAE agarose gel electrophoresis and Sanger sequenced (Biological Engineering (Shanghai) Co., Ltd.).
Quantitative Real-time PCR was performed using SYBR Premix Ex Taq™ (Tli RNase H Plus) (TaKaRa) according to the manufacturer's protocol. In the StepOnePlus™ Real-Time PCR system (Applied Biosystems) with the following PCR-cycling conditions: 95 °C for 10 s, 40 cycles of 95 °C for 30 s, 60 °C for 30 s with the concentration at 200 nM of each primer pair. The mean threshold cycle (Ct) was measured by three replicates of each gene, to test sample contamination and dimer formation, nuclease free water as template was included in all qRT-PCR experiments as the negative control. The 2 −ΔΔ Ct method 45 was used to calculate individual gene's relative expression levels. The 18 S rRNA gene (unigene 62750) was chosen as the endogenous reference gene in the qRT-PCR experiment as previously demonstrated 25,46,47 , and was expressed at the same level in the AMA (FPKM: 3.7717) and AFA (FPKM: 3.336) in FPKM analysis. All primers pairs for qRT-PCR are listed in Table S11.