Age influences the olfactory profiles of the migratory oriental armyworm mythimna separate at the molecular level

Background The oriental armyworm Mythimna separata (Walk) is a serious migratory pest; however, studies on its olfactory response and its underlying molecular mechanism are limited. To gain insights to the olfactory mechanism of migration, olfactory genes were identified using antennal transcriptome analysis. The olfactory response and the expression of olfactory genes for 1-day and 5-day-old moths were respectively investigated by EAG and RT-qPCR analyses. Results Putative 126 olfactory genes were identified in M. separata, which included 43 ORs, 13 GRs, 16 IRs, 37 OBPs, 14 CSPs, and 3 SNMPs. RPKM values of IR75d and 10 ORs were larger than co-receptors IR25a and ORco, and the RPKM value of PR2 was larger than that of other ORs. Expression of GR1 (sweet receptor) was higher than that of other GRs. Several sex pheromones activated evident EAG responses where the responses of 5-day-old male moths to the sex pheromones were significantly greater than those of female and 1-day old male moths. In accordance with the EAG response, 11 pheromone genes, including 6 PRs and 5 PBPs were identified in M. separate, and the expression levels of 7 pheromone genes in 5-day-old moths were significantly higher than those of females and 1-day-old moths. PR2 and PBP2 might be used in identifying Z11-16: Ald, which is the main sex pheromone component of M. separata. EAG responses to 16 plant volatiles and the expression levels of 43 olfactory genes in 1-day-old moths were significantly greater than that observed in the 5-day-old moths. Heptanal, Z6-nonenal, and benzaldehyde might be very important floral volatiles for host searching and recognized by several olfactory genes with high expression. Some plant volatiles might be important to male moths because the EAG response to 16 plant volatiles and the expression of 43 olfactory genes were significantly larger in males than in females. Conclusions The findings of the present study show the effect of adult age on olfactory responses and expression profile of olfactory genes in the migratory pest M. separate. Electronic supplementary material The online version of this article (doi:10.1186/s12864-016-3427-2) contains supplementary material, which is available to authorized users.

. In addition, host plant volatiles but not light conditions significantly affect the level of nocturnal activity in M. separata caterpillars [12,13]. Only a few studies on M. separata olfaction have been conducted, which include the antennal sensilla type [14], the EAG response to sex pheromone [15], and the remote-sensing sex pheromone trap for real-time monitoring of M. separata [16].
Recently, there has been significant progress in identifying olfactory genes in moths such as Bombyx mori [17][18][19], Manduca sexta [20], Helicoverpa armigera and H. assulta [21], and Spodoptera frugiperda [22] and S. litura [23]. To date, only three ORs (two PRs and one ORco) and one PBP were identified in M. separata antennae [24]. To gain more insights into the molecular mechanism underlying M. separata olfaction and into the olfactory mechanism of insect migration, the present study performed antennal de novo transcriptome analysis to identify olfactory genes and RT-qPCR to compare the expression profiles between sexes and different premating statuses (1-day-and 5-day-old). The link between OR gene expression and chemosensory responses as measured by electroantennography is also discussed.

Insects
Eggs of M. separata were purchased from the Chinese Academy of Agricultural Sciences and were raised in odorless insect incubators with a temperature of 25 ± 1°C, humidity of 70 ± 7%, and photoperiod (L: D) of 14 h:10 h (from 18:00). Larvae were fed with fresh maize leaves until pupation. Adults were checked at 17:00 daily. Male and female adults were divided into different insect rearing cages (30 cm × 30 cm × 30 cm) and were fed daily with fresh 10% glucose water.

EAG analysis
Recordings of whole-antenna electrical activity in response to chemicals were performed according to the standard technique described elsewhere [25]. Each chemical (Additional file 1: Table S1) (10 μL) diluted with paraffin oil to concentrations of 10 −2 and 10 −4 (V/ V) was added onto a filter paper (30 mm × 4 mm) as stimulus, and paraffin oil was used as control. Six antennae of 1-day-and 5-day-old moths were tested after exposure to each chemical for 3.5 h in the dark. As the response of antennae declined during the course of experiment, response to the 1% concentration (V/V) of (Z)-3-hexenyl acetate was used as reference and responses to all tested chemicals were standardized.

Sample collection and total RNA extraction
Antennae of 1-day-and 5-day-old 25 adult M. separata (males and females separately) were collected after 3.5 h in the dark. Samples from each group were immediately homogenized in TRNzol-A + (TIANGEN Biotech, Beijing, China) on ice, and total RNA was extracted according to the manual. The concentration and purity of total RNA were determined by using a spectrophotometer NanoDrop2000 (ThermoFisher, Waltham, MA, USA). RNA with an A 260 /A 280 ratio between 1.75-2.05, an A 260 /A 230 ratio > 1, and a concentration > 400 ng/μL were used in the following experiments. Extracts were treated with DNase I (Takara, Kusatsu, Shiga, Japan) to remove any DNA. RNA extractions were replicated three times.

Analysis of de novo transcriptome
The cDNA libraries for transcriptome analysis were prepared using TruSeq SBS Kit v3-HS (Illumina, San Diego, CA, USA) following the manufacturer's recommendations. Briefly, oligo (dT) magnetic beads were used to isolate poly (A) mRNA from collected total RNA. A fragmentation buffer was added to break the mRNA into short fragments. A random hexamer primer was used to synthesize the first-strand cDNAs using short RNA fragments as templates. The second-strand cDNAs were synthesized using a buffer, dNTPs, RNase H, and DNA polymerase I. After purification, the short cDNAs were linked to sequencing adapters. The libraries were sequenced using a paired-end transcriptome platform Illumina NextSeq500 (Illumina, San Diego, CA, USA) with 2 × 150 bp. Adapters and low-quality bases (base quality < 20) containing reads were discarded from the raw reads to obtain clean reads for analysis. De novo transcriptome assembly was conducted with the short reads assembling program, Trinity [26]. BLASTx alignment (E value < 0.00001) between unigenes and protein databases, including NCBI non-redundant protein sequences, Gene Ontology [27], Kyoto Encyclopedia of Genes and Genome [28], eggNOG [29], and SWISS-PROT was successively performed.

Profiling analysis of gene expression based on antennal transcriptome
The gene expression level was calculated using the RPKM method based on the results of antennal transcriptome analysis [41], in which the number of mapped reads per million reads for a gene is divided by the length of that gene. When there were multiple transcript variants for a gene, the longest one was used in calculating its expression level. The identification of differentially expressed genes between male and female antennae was conducted by using DESeq (version 1.18.0) [42]. P-values (p-value < 0.05) and fold-changes (|fold change| > 2) were calculated to determine the differentially expressed genes in the present study.

RT-qPCR of olfactory gene expression
Single-stranded cDNAs were synthesized from 1 μg of total RNA with ReverTra Ace qPCR RT Kit (Toyobo, Kita-ku, Osaka, Japan) following the manufacturer's recommendations. qRT-PCR was performed with SsoFast™ EvaGreen® Supermix (Bio-Rad, Hercules, California, USA), following the manufacturer's protocols, in a CFX-96™ PCR Detection System (Bio-Rad, Hercules, California, USA). The PCR primers used are listed in Additional file 1: Table S2. Actin and AK were used as reference genes. The difference in gene expression was measured by the 2 -ΔΔCt algorithm [43]. For every gene, expressions among all studied tissues were measured with 1-dayold female antennae as control (IR1 with 1-day-old male antennae as control). Then, the logarithm of all data to base 10 was used in generating a heat map. RNA extraction was repeated three times for each sample and two or more RT-qPCR replicates were prepared for each sample.

Data analysis
Data were analyzed using SPSS 17.0 (SPSS Inc., Chicago, IL, USA). Significance between two samples was determined by using an independent-samples ttest. The critical P value for each test was set at 0.05. A heat map was generated by using the software Per-mutMatrixEN [44].

Sequencing and unigene assembly
Using the Illumina NextSeq 500 sequencing system, a total of 75,705,108 (11,355,766,200 Figure S1). For the evolutionary genealogy of genes, eggNOG annotation was used to assess the transcriptome (Fig. 1), and 5105 (10.75%) unigenes were assigned to signal transduction mechanisms. Among  unigenes, male antennae showed 28 genes with a higher expression, whereas 118 genes were upregulated in the female antennae (Fig. 2).

Expression of olfactory genes between male and female antennae
The expression levels of 43 olfactory genes, including 10 CSPs, 6 GRs, 2 IRs, 16 OBPs and 9 ORs in male antennae were significantly higher than that in female antennae, and the difference among 3 OBPs (OBP13, OBP15, and OBP3) and IR1 was > 10-fold (Fig. 10). The expression level of 38 olfactory genes, including CSP9, 4 GRs, 8 IRs, 9 OBPs and 16 ORs in female antennae was significantly higher than that in male antennae and the difference in OBP9 was higher by 10-fold. The expression level of 4 IRs, 4 OBPs and 6 ORs was the same between male and female antennae. The difference in expression levels of 16 olfactory genes, including 3 GRs, IR2, 3 CSPs, 3 OBPs, and 6 ORs between male and female antennae was influenced by their pre-mating status (Fig. 10).
The expression levels of 4 PBPs (except for PBP3) and 5 PRs (except for PR6) were significantly higher in male antennae than that in female antennae, and the difference between PBP2 and the 3 PRs (PR1, PR2, and PR4) was 10-fold higher. The expression level of PBP3 in 5day-old females was 10-fold higher than that in males. The expression level of PR6 in 1-day-old males was significantly higher than that in females, but that in 5day-old males was significantly smaller than in the female counterpart (Fig. 10). ORs (OR12, OR3, and OR21) in 5-day-old moths were significantly higher than that observed in 1-day-old moths, and the difference of OBP15 was 10-fold higher (Fig. 10). The expression levels of 43 olfactory genes, including 6 CSPs (CSP3, CSP4, and CSP7), 9 IRs (IR25a and IR2), 7 GRs, 8 OBPs (OBP14), and 13 ORs (OR6) in 5-day-old moths were significantly smaller than the 1day-old moths, and difference among CSP9, OBP9, and OR22 was 10-fold higher. The expression levels of 14 olfactory genes, including GR12, 8 OBPs (GOBP1, OBP16, and OBP21), and 5 ORs (OR10, OR4, and ORco) in 1-day-old moths were similar to those of 5-day-old moths. The differential expression of the 28 olfactory genes, including 5 CSPs (CSP1), 3 GRs, 3 IRs (IR75d), 7 OBPs (GOBP2), and 10 ORs between different premating statuses was influenced by gender, particularly for CSP11, CSP12, and OBP3 (Fig. 10). Three PBPs (PBP2, PBP3, and PBP4) and four PRs (except for PR4 and PR5) were upregulated in 5-day-old moths relative to that in 1-day-old moths, and a 10-fold difference in expression between PBP2 and PR3 was observed (Fig. 10). The expression level of PR4 in 1-dayold females was significantly higher than that in 5-dayold females. The expression level of PBP5 was the same between 1-day-old and 5-day-old moths. The difference in the expression level of PBP1 and PR5 between different pre-mating statuses was influenced by gender.
In 6 of 7 sex pheromones (except Z9E11-14:OAc) with evident EAG responses (>0.4 mV), the male responses were significantly greater than that in females, and the responses to the 10 −2 concentration were significantly larger than that observed with the 10 −4 concentration (Fig. 11). The EAG responses of 5-day-old males to Z11-16:Ald, Z9-16:Ald, and Z13-18:OAc were significantly greater than 1-day-old males, but the EAG responses of 5-day-old males to Z11-16:OH and E\Z12-14:OAc were significantly smaller than 1-day-old males. Female EAG responses to Z9E11-14:OAc were significantly greater than males, and the EAG responses of 5-day-old females were significantly larger than 1-dayold females (Fig. 11).
The EAG responses at 10 −2 concentration were significantly greater than that using 10 −4 . Furthermore, the EAG responses of males were significantly greater than that in females, and 1-day-old males were significantly greater than that in 5-day-old males. At the same time, the EAG responses of 5-day-old females to E2-5:Ald and heptanal were significantly larger than that in 1-day-old females at 10 −2 concentration (Fig. 12).

Olfactory genes of M. separata
We identified a total of 126 putative olfactory genes in the antennae of M. separata, including 43 OR genes, 13 GR genes, 16 IR genes, 37 OBPs, 14 CSPs, and 3 SNMPs. The number of ORs, IRs, OBPs, and CSPs in M. separate is smaller than that in B. mori [19,45], which was conducted using the whole genome, but was the same as that of other moths that were studied using the same protocol (antennal transcriptome), which included Agrotis ipsilon [46], C. suppressalis [37], H. armigera [47], O. furnacalis [35,36]. However, the number of GRs and SNMPs was larger than those of most moths studied.
The maxillary palps, which harbor considerably fewer sensilla than antennae, are believed to specialize in taste reception in various moth species [48]. However, 14 gustatory receptor genes were identified in M. separate antennae, indicating that the antennae might be an important taste organ for M. separata. As shown in a former study, sensilla chaetica were distributed around each antennal segment of M. separata [14], which was with a terminal pore and with a shape that was higher and wider than that of other sensilla and have been suggested to have contact/chemoreceptor functions [49]. Sugar as the supplementary nutrient is essential for egg- ripening [8] and as the energy supplement influences the duration of migration [9]. Five sweet receptors were identified in M. separata that specifically recognize glucose, galactose, fructose, mannose, sucrose, maltose, trehalose, raffinose, glycerol, and mannitol [50]. GR1, a sweet receptor with higher RPKM than other GRs might be used to recognize the most important sugar for M. separata. In addition, the expressions of GR1 in 1-dayold males were larger than that in females and 5-day-old males. Two conserved members of SNMP (SNMP1 and SNMP2) were detected in various moth species [51]. A new SNMP member SNMP3 that has been identified in S. exigua recently and other four moths [52] was also reported in M. separata. SNMP3 shares only less than 30% identity with SNMP1 and 2 [52]. Different insects have evolved various feeding behaviors in host plants and the evolution of olfactory acuity and discriminatory power of insect must be consistent with the difference in plant odor [53]. Approximately 5 ORs, 2 GRs, 2 IRs, 7 OBPs, and 4 CSPs of M. separata were not effectively clustered with those of Lepidoptera (bootstrap values <50). These olfactory genes might be correlated with the special habitat of M. separata, including attacking rice, maize, sorghum, wheat, and other  gramineous plants. M. separata larvae have evolved a specific detoxification system for benzoxazinoid, which is part of the chemical defense system of graminaceous plants [54].

Co-receptor for odorant detection in M. separata
IRs act in combinations of up to three subunits, comprising individual odor-specific receptors and one or two broadly expressed co-receptors, and heteromeric IR complex formation is necessary and sufficient for mediating odor-evoked electrophysiological responses in vivo and in vitro [55]. IR25a and IR8a function as coreceptors because of their higher sequence identity compared to the other IRs [32], expression in all or most OSN cell bodies, and co-expression with divergent IRs in former studies [56][57][58][59]. However, in the present study, the expression levels of IR8a (RPKM = 1.9) and IR25a (RPKM = 13.2) were evidently smaller than that of IR75d (RPKM = 509.1). IR75d has been identified in various insects such as Grapholita molesta [60], Athetis dissimilis [61], Rhynchophorus ferrugineus [62], B. mori, Acyrthosiphon pisum, Pediculus humanus humanus, and dipterans [32]. In Heliconius melpomene, the IR75d gene is highly expressed across various tissues and sexes [63]. In D. melanogaster, IR75d is expressed in three different OSNs that are housed in three coeloconic sensilla (ac1, ac2, and ac4) [64] and is the only ORN class that occurs in more than one sensilla subtype [65]. In addition, IR75d is co-expressed with IR84a and IR76a in ac4 sensilla [66]. These findings suggest that IR75d might be a co-expressed receptor in M. separate.
Orco is the common olfactory receptor co-receptor in insects and is present in apparently all OR-expressing ORNs [67]. The expression level of Orco is always higher than that of other ORs [23]. It has been proposed that Orco is a membrane localization protein which  stabilizes ORs in the dendritic membranes, as well as a chaperon molecule that facilitates in the correct protein folding of ORs and in forming heteromers with the particular OR as a cation channel [68,69]. Downregulation of the Orco transcript not only results in weaker EAG responses [58,[70][71][72], but also diminishes the significant behavior preference [72][73][74] for host volatiles and pheromones. However, the RPKM value of ORco (5.5) in M. separata was less than 9 general ORs and 2 PRs. The low expression level of ORco in M. separata provides allows odorant receptors to mediate odorant responses also in the absence of Orco. When ORs from different species were expressed without Orco in heterologous expression assays such as Drosophila Schneider 2 (S2) and S. frugiperda 9 (SF9) cells, these also evoked ligandspecific responses [75]. In addition, the functional expression of BmOR-1 and BmOR-3 in modified HEK 293 cells was possible without co-receptor BmOR-2 and BmOR-3 being heterologously expressed. Furthermore, modified HEK 293 cells are highly sensitive and selective responsiveness to bombykal [76]. The Orco agonist VUAA1 did not collectively alter mating efficiency, the number of eggs produced, and the number of hatched nymphs in the bed bug, Cimex lectularius L [77]. In fact, Orco is not essential for odorant detection when ORs are successfully inserted into the plasma membrane [75].

Pheromone recognition of M. separata
The sex pheromone of M. separata isolated from Japan and Taiwan was identified as a blend of (Z)-11-hexadecenyl acetate (Z11-16:Ac) and (Z)-11-hexadecenol (Z11-16:OH) [78,79]. However, Z11-16Ald is the main female sex pheromone in mainland China [10,11], and male moths had the highest peak response to Z11-16:Ald [15]. In the present study, the EAG response of male oriental armyworm to Z11-16:Ald was significantly greater than that to other pheromones. In addition, the male EAG response to Z11-16:Ald was significantly greater than that in females, and the EAG response to the 10 −2 concentration was significantly greater than that to the 10 −4 concentration. However, the EAG responses to Z11-16:Ac was very small (<0.2 mV) and were not significantly influenced by concentration, sex, and pre- Fig. 10 Expression levels of olfactory genes in male and female antennae with different day time measured by RT-qPCR. Gene expression was calculated relative to the Actin and AK as reference genes and expression in 1-day-old female antennae was arbitrarily defined as control for all genes. Gene expression in other tissue were normalized to 1-day-old female antennae. Log scale of gene expression was used to generate heatmap. F1, 1-day-old female, F5, 5-day-old female, M1, 1-day-old male, M5, 5-day-old male mating status. Z11-16:Ald is the main female sex pheromone in M. separata. The calling behavior of females is the most frequent and the sex pheromone titer in the sex gland is the highest in 4-5 day old M. separata [80]. In the present study, the EAG responses to Z11-16:Ald of 5-day-old males were significantly greater than that of 1-day-old males, which was consistent with the female calling period of M. separata. In all candidate genes for sex pheromone recognition, PR2 and PBP2 were the most abundant PR and PBP, respectively, and were specifically expressed in male antennae. In addition, the expression level of PR2 and PBP2 in 5-day-old males was 10-fold higher than those found in 1-day-old males. The consistency between the EAG response to Z11-16:Ald and the expression profile of PR2 and PBP2 suggests that these two genes might be used in recognizing Z11-16:Ald, the main female sex pheromone.
In addition to Z11-16:Ald, there were also EAG responses to other sex pheromones such as Z9-16:Ald, Z11-16:OH, and E\Z12-14:OAc. Male EAG responses to these pheromones were significantly greater than female EAG responses, and the responses to the 10 −2 concentration were significantly greater than that to the 10 −4 concentration. These sex pheromones might be identified by male M. separata. Besides PR2 and PBP2, there were also other 5 PRs and 4 PBPs in the antennae and the expression levels of 3 PBPs and 4 PRs in male antennae were significantly higher than the females. These genes might be used in recognizing former sex pheromones. For example, MsepPR1 is used as a specific pheromone receptor for Z11-16:Ac, which is the main sex pheromone component of the Japanese M. separata [24]. EAG responses of 5-day-old males to Z9-16:Ald was significantly greater than that of 1-day-old males in the present study, and Z9-16:Ald might be a minor sex pheromone component in M. separata female glands [15]. In accordance with the EAG response, the expression levels of 2 PBPs (PBP1 and PBP3) and PR3 were higher in 5-day-old male moths than in the 1-day-old moths. E\Z12-14:OAc is the main female sex pheromone of O. furnacalis, which is a major pest of maize [81]. M. separata males might be capable of distinguishing conspecific females from other species by recognizing their sex pheromone components or might search host plants with female sex pheromone of moths that infest the same host as that of M. separata. For example, M. separata might find maize plants by recognizing E\Z12-14:OAc.

Host recognition in M. separata
At least four distinct large-scale and long-distance migration events of M. separata between overwintering Fig. 11 EAG response of male and female moth with different day time to sex pheromone with 10 −4 (left) and 10 −2 (right) concentration. F1, 1-day-old female, F5, 5-day-old female, M1, 1-day-old male, M5, 5-day-old male sites in southern China and northern temperate zones occur annually, including two northward displacements in the spring and early summer and two southward displacements in summer and fall [3]. The migration takes place in 1-2 days after eclosion for the orient armyworm [5]. Adults need a sufficient amount of sugar as the energy supplement for flight, and sugar is a crucial factor that influences the duration of migration [9]. In addition, adults also need a sufficient amount of sugar for survival; otherwise, these soon perish [8]. Under natural conditions, the adult armyworm feeds on flower nectar or honey-dew of aphids [50]. In the present study, EAG responses of 1-day-old moths to 16 plant volatiles were significantly larger than the responses of 5-day-old moths, and expression level of the 43 olfactory genes (6 CSPs, 9 IRs, 7 GRs, 8 OBPs, and 13 ORs) in 1-day-old moths were significantly higher than those in the 5-day-old moths. These plant volatiles might be used by an adult armyworm to search for food, and might be mediated by these olfactory genes, which were differentially expressed in 1-and 5-day-old moth antennae. In particular, the EAG responses to heptanal, Z6-nonenal, and benzaldehyde were greater than 1.2 mV, which is the EAG response to sex pheromones. Heptanal, Z6-nonenal, and benzaldehyde are flower volatile compounds [82,83]. These three flower volatiles might be very important to adult armyworms as they search for food. The RPKM values of the 7 OBPs and 4 CSPs were > 1000, whereas that of the 6 ORs and 2 IRs were > 10. Therefore, these olfactory genes might play a critical role in the recognition of host volatiles.
Spermatids are observed at the time of male M. separate emergence [8], whereas oocytes ripen only after the female moth has consumed a sufficient amount of food [84]. In accordance with the larger sugar requirement of the female, the expression levels of 38 olfactory genes (CSP9, 4 GRs, 8 IRs, 9 OBPs, and 16 ORs) in females were significantly higher than those in the males (Fig. 10). However, no plant volatiles with higher female EAG responses than that in the males were observed in the present study. The inconsistency between gene expression and EAG response might be attributable to the loss of important chemicals. In fact, identifying all stimulating chemicals may be a relatively difficult task. A panel of 480 odorants, including esters, acids, aldehydes, ketones, alcohols, pyrazines, aromatics, terpenes, and sulfur compounds, against 21 larval odorant receptors (more than 10,000 odorant-receptor combinations) were conducted in fly, but no strong odorants were identified for Or2a and Or49a [85]. Other than sex pheromones, host plant volatiles might also be used by sexually mature males to find females for mating. In Lepidoptera, several cases of synergism between plant semiochemicals and pheromones have been observed both in the laboratory and field [86]. In our study, the male responses to 16 plant volatiles with high EAG responses (>0.4 mV) were significantly greater than the female responses, and the expression levels of 43 olfactory genes (10 CSPs, 6 GRs, 2 IRs, 16 OBPs, and 9 ORs) in males were significantly greater than in the females. These results prove that some plant volatiles are very important to male moths, and using of these volatiles might be a strategy to optimize mating opportunities. However, our understanding of how insect responses to plant stimuli aid in mating are limited [87].
In the present study, the EAG response to volatile organic compounds of larvae host plants was very small, which include β-myrcene, (Z)-3-hexenyl acetate, nonanal, linalool, β-caryophyllene, (E)-2-hexenal, and (Z)-3hexen-1-ol from wheat [88,89]. These low EAG responses may be attributable to the pre-mating status of the moths used in the present study. Only 1-and 5day-old unmated male and female moths were used in our study. The most important behavior of 1-and 5day-old moths involved finding their nutrition for sex maturation as well as searching for a mating partner. Therefore, volatiles of larvae host plant might not be required by 1-and 5-day-old unmated moths.

Additional files
Additional file 1: Table S1. Chemicals tested in the Electroantennogram recordings. Table S2. Primers used in the study. (DOC 267 kb) Additional file 2: Figure S1. Size distribution of all 41,056 unigenes assembled from the pooled M. separata RNA extract. Figure S2. Aligned putative full ORF of OBP gene sequences of M. separata. Six conserved cysteines were highlighted by blue color. Figure S3.