A. lucorum protein sequences were retrieved from the NCBI Genome database (https://www.ncbi.nlm.nih.gov/assembly/GCA_009739505.2/) (32). Based on previous studies and records in Flybase (http://www.flybase.org/) (33), the GPCRs of Drosophila melanogaster (34), A. pisum (22), Bombyx mori (17), Tribolium castaneum (16), and Pediculus humanus humanus (18) were collected. By using D. melanogaster GPCRs as references and A. lucorum protein sequences as queries, BLASTP searches (35) were performed with a cut-off e-value of 1e−5 to look for all GPCR candidates. Then, seven-transmembrane (7TM) domain and annotation information was adopted as the basic criteria for all GPCR candidates. The GPCR candidates in which the number of 7TM domains was more than four or the annotation information indicated it was a GPCR were retained. The remaining GPCR candidates were also confirmed by means of BLASTX analysis in the UniProtKB/Swiss-Prot database. Using all GPCRs that we collected, pre-phylogenetic analysis with the maximum likelihood method was the final criteria to remove non-GPCRs from candidate pools. If a candidate showed fewer genetic relationships with known GPCRs by phylogenetic analysis and the hit sequences in BLASTX analysis indicated they were not GPCRs, they were classified as a non-GPCR and removed from our analysis.

The 7TM domains for all GPCR candidates were predicted with the server TMHMM (v2.0) (36) from the Centre for Biological Sequence Analysis (http://www.cbs.dtu.dk/services/TMHMM/). Functional annotations of the target proteins were done using InterProScan (37). In addition, the chromosomal location of each GPCR candidate was extracted from the genome annotation file of A. lucorum.

Partial GPCRs of R. prolixus and C. lectularius that also belonged to Heteroptera were also obtain based on previous study (24, 27). GPCRs from D. melanogaster, A. pisum, C. lectularius, and R. prolixus were assigned to a family/subfamily according to previous results (22, 24, 27, 34). Putative A. lucorum GPCRs were classified into different families/subfamilies according to the families to which their orthologous proteins were assigned. Amino acid sequences of the putative A. lucorum GPCRs in each family/subfamily were aligned with receptors of the same family/subfamily in D. melanogaster, A. pisum, C. lectularius, and R. prolixus using MAFFT v7 (38). Phylogeny tests were accomplished using the bootstrap method with 1,000 replications to reconstruct maximum likelihood (ML) trees using IQ-TREE (39) and the best-fit tree model was determined with ModelFinder (40). It should be noted that the GPCRs of R. prolixus and C. lectularius were uncompleted, which were composed of opsins, biogenic amine receptors, and neuropeptide GPCRs. For the Drosophila sequences, the name of the GPCRs were used, while for A. pisum and R. prolixus, the protein names were same as in previous work (22, 27), and for C. lectularius, the accession numbers in NCBI were used. The GPCRs of A. lucorum identified in this work were numbered according to their families.

To study the expression profiles of the GPCRs, a total of 39 transcriptome data of A. lucorum were downloaded from the genome project of A. lucorum (Accession: PRJNA526332) in the NCBI Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra/) (30, 41), which included egg and different tissues (leg, head, body, mouthpart, wing, and gut) of nymphs and adults. Each tissue contained three biological replicates. In detail, we downloaded the SRA data first and then we used an SRA-Toolkit to split the paired-end reads. Clean reads were obtained from the raw data using Trimmomatic (42) to remove reads with quality scores lower than 10 and adapter sequences. To analyze gene expression profiles, clean reads of each sample were mapped to A. lucorum gene sets using hisat2 (43), and then the TPM value (44) of each putative GPCR gene was calculated with featureCounts (45). These TPM expression values were scaled and served to generate a cross-sample normalized trimmed mean of the M-values (TMM) gene expression matrix (46). Finally, the heatmap was drawn in ITOL (https://itol.embl.de) (47) using the normalized matrix. The value used for each sample was the mean of three independent biological replicates.

MCScanX (48) was used to classified the duplication types of different duplicate GPCR genes. First, the homology with different genes in the genome of A. lucorum was determined by a whole-genome BLASTP analysis with a max target seqs of 5 and a cut-off e-value of 1e−5. Then, the homology with different genes and the chromosomal location were combined and all genes were classified into various types, including the segmental duplication, and tandem duplication. Finally, the duplication types of GPCRs were extracted based on these results. All visualized works were accomplished in TBtools (49).

Insect family-A GPCRs include opsins, biogenic amine receptors, neuropeptide and protein hormone receptors, and purine GPCRs (17, 18, 22). In this study, 98 family-A GPCRs were identified in the genome of A. lucorum, and these receptors were composed of seven opsins, 30 biogenic amine receptors, 58 neuropeptide and protein hormone receptors, and three purine GPCRs ( Tables 1 , 2 ).

Family-B GPCRs play vital roles in many biological processes, including growth, development, and reproduction. They are characterized by long N-terminal domains, and they form a small group of receptors that are structurally and functionally divergent from other groups of GPCRs (56). Within this family, family-B GPCRs can be further subdivided into three subfamilies: B1–B3, which are greatly divergent in both function and structure. In total, 21 family-B GPCRs were identified from the genome of A. lucorum in this study, which consisted of nine B1 subfamily members, three B2 subfamily members, and nine B3 subfamily members ( Table 3 and Figure 3 ).

The B1 subfamily is made of largely classical hormone receptors. It comprises three types of hormone receptors in D. melanogaster: diuretic hormone 31 receptor (DH31-R/hector), CRF-like diuretic hormone 44 (DH44-R), and pigment dispersing factor receptor (Pdfr) (57–60). In our study, all three types of hormone receptors were identified. The parathyroid hormone receptor (PTHR), which is involved in the calcium and phosphate homeostasis and bone growth in vertebrates, is also a subfamily-B1 GPCR (58). There are two PTHRs in mammals, which are involved in calcium and phosphate homeostasis and bone growth (61, 62). In insects, PTHR-like (PTHRL) have been identified from T. castaneum, A. mellifera, P. h. humanus, and N. lugens (18, 26, 57), but its counter-parts in D. melanogaster, B. mori, A. pisum, and A. gambiae are not found (17, 22, 34, 63). T. castaneum has two distinct PTHRLs (57), N. lugens possesses a pair of homologous PTHRLs (26), and A. mellifera only has one PTHRL (57). In our study, we also identified one PTHRL, B9, which shared a low e-value (1e−104 and 1e−111) with two PTHRLs in N. lugens. These results showed that genes coding for PTHR are divergent among insects.

The B2 subfamily is characterized by a long extracellular N-terminal domain and a GPCR proteolytic site (57, 58, 64). Based on phylogenetic analysis and sequence similarity, three receptors (B10–12) were classified in the B2 subfamily, which correspond to a calcium-independent receptor for α-latrotoxin (Cirl), starry night (stan), and CG15744, respectively. However, the orthologs for CG11318 and CG15556 were not identified in the genome of A. lucorum.

There is only one group of receptors in the B3 subfamily; i.e., Methuselah (mth)/Methuselah-like (mthl) (57, 58). This gene family is involved in the modulation of life span and stress responses. No counterpart for the mth gene family has been identified in vertebrates (57). In insects, the number in the B3 subfamily is highly variable (18, 57). In D. melanogaster, this subfamily can be divided into two groups based on their structure, 12 mth ectodomain‐positive members (mth, mthl2, mthl3, mthl4, mthl6, mthl7, mthl8, mthl9, mthl10, mthl11, mthl12, and mthl13) and four mth ectodomain‐negative members (mthl1, mthl5, mthl14, and mthl15) (65). In our study, nine receptors (B13–21) were identified in this family. Based on phylogenetic analysis, B13 and B20 may belong to the mth ectodomain‐positive group, and others may be members of the mth ectodomain‐negative group. At the mRNA level, B13, B18, and B21 showed a higher expression level than other mthl members, which indicated these three members of the B3 subfamily may play a more important role in A. lucorum.

Family-C GPCRs possess a large ligand-binding extracellular domain and form constitutive dimers (34, 66–68). There are three types of GPCRs in family-C, the glutamate and γ-amino butyric acid (GABA-B) receptors, the bride of sevenless (boss-type) receptors, and the metabotropic glutamate (mGlu) receptors. Until now, nine family-C GPCRs from D. melanogaster and seven from A. pisum have been reported. By using these reference sequences, 10 family-C members ( Table S1 nad Figure S2 ) of A. lucorum were identified here.

Like A. pisum (22), there are two GABA-B receptors (C1 and C2) in A. lucorum. C1 shares a 78% sequence similarity with D. melanogaster GABA-B-R1, while C2 has a 44% sequence identity with D. melanogaster GABA-B-R2. The orthologous gene to D. melanogaster GABA-B-R3 has not been found in two Hemiptera insects. The boss-type receptor was first identified as a ligand for sevenless tyrosine kinase, which was involved in eye differentiation in D. melanogaster. Subsequently, boss has been implicated in the glucose-response (69, 70). It has been reported in D. melanogaster (34), A. gambiae (63), A. pisum (22), and T. castaneum (16), but not in B. mori (17), A. mellifera, N. vitripennis, and P. humanus corporis. Here, C3 was identified as the orthologous gene to boss. These results indicated that boss has been randomly lost in insects during their evolutionary process. Moreover, three mGlu receptors (C4, C5, and C6) were found in A. lucorum, whereas there are only two mGlu receptors in D. melanogaster (34) and one mGlu receptor in A. pisum (22). Small expansions of A. lucorum mGlu receptors have been observed. There are some unclassified receptors in this family. C7 was the orthologous gene to smog. C8 and C9 showed 53 and 38% sequence identities with their counterparts in D. melanogaster, respectively. As shown in Figure S2 , C9 showed a high expressed level in egg, adult head, and all nymph tissues, except leg, while the function of its orthologous gene is unclear. C10 was an orphan receptor that has not been reported in D. melanogaster and A. pisum (22, 34).

Family-F GPCRs comprise the frizzled gene family and the smoothened gene (34, 71). In this study, four putative A. lucorum GPCRs were identified in the frizzled/smoothened GPCR family, which are orthologous to D. melanogaster fz, fz2, fz3, and smo ( Table S1 and Figure S3 ). Our results indicated that orthologs for D. melanogaster fz4 were missing in A. lucorum. Among family-F, F4 (smo) showed the highest expression in the egg with a TPM of 61.

Compared with other well-studied insects, we noticed that the number of genes coding for GPCRs is obviously larger than for other insects, especially expanded among the biogenic amine receptors, neuropeptide and protein hormone receptors, and B1 subfamily ( Table 4 ). There were 26 GPCRs duplicated in A. lucorum. Twenty-three of them were classified into the three subfamilies mentioned above. By MCScanX analysis, we found six tandem duplication events occurred among Rh6, Rh7, AstA-R, ETHR, FMRFaR, and AdoR, while most GPCR genes duplicated dispersedly. Considering the location of GPCR genes on chromosomes, it suggested that the duplication of GPCR genes mainly occurred as independent duplications and transitions ( Figure 4 ).

The duplicate biogenic amine receptors in A. lucorum included 5-HT1A, 5-HT1B, 5-HT2A, 5-HT2B, 5-HT7, Dop2R, mAChR-B, Oamb, Octbeta2R, and Octbeta3R. These biogenic amine receptors can regulate many behaviors including flight and fight, learning and memory, sleep and wakefulness, feeding, and social and reproductive behaviors (72–74). For example, 5-HT1A was related to locomotor activity in B. mori. Injecting the antagonist of the Bm5-HT1A receptor into larvae caused slow or weak motility, and adults had lowered courtship vitality or moving speed (75). mAChRs have also been reported to be critical in regulation of locomotory behavior in Drosophila (76). In addition, 5-HT1B mediates hemocyte phagocytosis and serotonergic signaling performs critical modulatory functions in immune systems. Moreover, 5-HT7 and Dop2R were shown to be associated with learning ability (77, 78) and octopamine receptors were required for ovulation in D. melanogaster (79). The ligands of neuropeptide and protein hormone receptors and the B1 subfamily belong to neuropeptides, which also play an important role in the regulation of development, reproduction, feeding, courtship, aggression, olfaction, locomotor activity, circadian rhythm, and many other physiological processes in insects (21, 80, 81). The gene expansion of these three subfamilies indicated that A. lucorum had a more complex peptidergic signaling system.

Expansions of genes associated with omnivorousness and mesophyll feeding, such as those related to digestion, chemosensory perception, and detoxification, were also observed in A. lucorum (32). Gustatory receptors (Grs) and odorant receptors (Ors) are thought to be the most important chemosensory receptors. Like GPCRs, Ors, and Grs are seven-transmembrane domain receptors but belong to the chemosensory 7tm receptor superfamily (82, 83). It has been suggested that the cause of gene expansion in GPCRs might be similar to that of chemosensory receptors, also to better adapt to the environment (32). A. lucorum is found in natural and agricultural ecosystems throughout the world (30), and many of them are generalists, exhibiting diverse feeding habits or preferences (e.g., feeding on leaf, stem, inflorescences, nectar, pollen, and fruit) (32). These results indicated the complex peptidergic signaling system is more favorable for A. lucorum to adapt to multiple living environments and multiple hosts.

FMRFamide (FMRFa) is a cardioexcitatory peptide that was first isolated from the nervous system of the clam, Macrocallista nimbosa (84), and is active as a tetrapeptide only in mollusks and annelids. Since the discovery of FMRFa, peptides with extended length at the N-terminal portion have been reported, such as myosuppressin (Ms) (85). Here, nine receptors (A54–62) displayed a certainly evolutionary kinship with the FMRFaR of D. melanogaster, A. pisum, C. lectularius, and R. prolixus, while the MsR is missing in A. lucorum ( Figure 2 ). By reconstruction of the phylogenetic tree of MsR and FMRFaR with more species (17, 86), we found these receptors were closer to the known FMRFaRs ( Figure 7 ). Among these, A54 and A55 were clustered in a single clade with the FMRFaRs that had been identified in other insects, whereas the other seven receptors are clustered on the other single clade. The most similar proteins in UniProtKB/Swiss-Prot of these receptors are both D. melanogaster FMRFaRs (CG2114). Among these receptors, A54 was the ortholog to the insect FMRFaRs with the smallest e-value of 1.75E−130, and A55, which was adjacent located on chromosome 9, is a tandem duplication that occurred at the beginning of this receptor expansion. However, the other seven receptors (A56–62) were scattered across six chromosomes, which indicated these seven receptors might have arisen from transposition.

By searching in the genome of other heteropteran insects, we found there are only one or two FMRFaRs in each heteropteran species. We suggest A54 and A55 should be classified as FMRFaRs, while the others (A56–62) were named as FMRFaR-like for the moment. This branch might be another unknown GPCR or even contain the MsR. Recent research had found that FMRFaR stimulates intracellular calcium signaling through the IP3R and helps maintain neuronal excitability in a subset of dopaminergic neurons for positive modulation of flight bout durations (87), and the ligand can reduce spontaneous muscle contractions, such as in the intestinal muscle and the heart rate, which also have an impact on movements (88, 89). A. lucorum has great flight capacity and its adults can fly 151.3 km within 48 h (90). A large FMRFaR-like branch evolved in A. lucorum may help it maintain strong flight capability.

In our research, only 23 GPCR genes (17.3% of all GPCR genes) were expressed highly (TPM >10) in at least one tissue, while most GPCR genes showed a low expression level in A. lucorum. This result is consistent with previous studies in which most GPCRs showed a low endogenous expression level, even at the mRNA level (91–93). Certainly, a low expression level of GPCRs does not necessarily equate to functional insignificance (94).

Here, through transcriptome analysis, we found the expression levels of opsins were higher than in other subfamilies. Opsins that originated early in metazoan evolution mediate the response to visual stimuli primarily. When stimulated by light, opsins can activate a downstream signaling cascade by conformational change (95). Most opsin genes are expressed in photoreceptors, but there are opsins expressed in other tissues, suggesting some nonvisual functions (96, 97). In A. lucorum, opsin genes (A1–3) were expressed highly not only in head but also in leg, wing, and mouthpart, indicating these opsins may execute some nonvisual functions ( Figure 1 ). Among these, the LW opsin (A1) showed the highest expression levels in the head tissue of adults (TPM = 28,787, at least 20 times more than other GPCR; Figure 5 ). The peak absorbance of the LW opsin is 500–600 nm, which corresponds to yellow-green light. As night sets in, the natural ambient light is increasingly dominated by longer wavelengths (98, 99). The importance of LW opsin had been reported in many nocturnal insects (100, 101). The adults of the A. lucorum were mainly active from dusk to early morning (102). High expression levels of LW opsin may help the organism adapt to a low light environment. We found B10, orthologous to Cirl, is the most widely expressed GPCR gene, which can be tested in all tissues ( Figure 5 ). Cirl belongs to a unique branch of GPCRs and, specifically, is an adhesion GPCR (103, 104). The orthologs of Cirl have been discovered in almost all animals from invertebrates to vertebrates, including humans (105). There are three homologs of Cirl in most vertebrates (Cirl‐1, Cirl‐2, and Cirl‐3) and two in birds and worms, whereas there is only one homolog in insects—which is most homologous to vertebrate Cirl-2 (103, 106). The expression pattern of insects Cirl, which had been reported to be expressed in multiple tissues (103, 107), was also like vertebrate Cirl-2 (108). Although there is only one Cirl member in insect species, Cirl is still involved in multiple physiological processes, which can regulate sensory, developmental, reproductive, and immune functions in insects (104, 109). Here, B10 had been detected in all transcriptomic samples and its distribution range is wider than in T. castaneum and D. melanogaster (103, 106). This kind of expression pattern suggested B10 is crucial in the development of A. lucorum.

Within the neuropeptide and protein hormone receptor subfamily, LGRs are a distinct subgroup with important functions in development and reproduction (110). Three distinct types of LGRs have been defined based on their structural characteristics and they are distinguished by the number of leucine-rich repeat (LRR) motifs, the absence or presence of a low density lipoprotein receptor domain class A (LDLa) motif, and their type-specific hinge region. Generally, type B LGRs have about twice the number of LRRs compared to the other two types. An exclusive feature of the type C LGRs is the presence of at least one LDLa motif in the ectodomain. The more general type containing only one LDLa will be referred to as type C1, whereas type C2 contained more than one LDLa (111).

Type C2 LGRs were first discovered in echinoderms, mollusks, and in one insect species (Pediculus humanis corporis). In our study, we found it existed in all hemipteran insects that we studied. Combining recent work, we mentioned that type C2 LGRs are reported in many hemipteran insects and P. h. humanus (18, 23, 25, 26), and are lost in other orders of insects (16, 17, 63, 86) ( Figure 6 ). Until now, the presence of type C2 LGRs have been found in all Hemiptera insects in which their GPCRs have been identified (22, 23, 25, 26). Among insects, Phthiraptera is one of the orders most closely related to Hemiptera (112). Type C2 LGRs may be present in the common ancestor of these two orders. To clarify the distribution of type C2 LGRs in insects, we checked all protein sequences in the non-redundant protein sequences database (nr) of NCBI. The result showed, except for Hemiptera and Phthiraptera, type C2 LGRs were also founded in Zootermopsis nevadensis of Blattodea. In terms of functionality, LGRs have important functions in development and reproduction. Type A LGRs and type B LGRs are stimulated by large dimeric protein hormones (110), regulating the adult eclosion of insects, and cuticle tanning (113, 114), while type C1 LGRs are the receptors of insulin-like peptide 7 and insulin-like peptide 8 and they coordinate organ growth in D. melanogaster (115–117). At present, the function of type C2 LGRs is undetermined. The function of type C2 LGRs and the existence of type C2 LGRs in Phthiraptera need to be explored in future work.