The overwintering eggs of G. daurica were collected from the Xianghuang Banner grasslands (44°62′N, 115°80′E) of Inner Mongolia, China on September 2021, were brought back to the Research Center for Grassland Entomology (40°48′N, 111°42′E), Inner Mongolia Agricultural University (Hohhot, Inner Mongolia, China), and were maintained in intelligent light artificial incubator (PRX-350C, Ningbo HaiShuSaiFu experimental instrument Co., Ltd.) (25°C ± 1°C, RH = 70 ± 10%, L14:D10). After hatching out, the larvae were fed on Allium mongolicum in a incubator under the same conditions just as described by Tan et al. (2018). The healthy third instar larvae in consistent physiological state were collected for further experiments.

Bioassays were conducted using a leaf-dip method based on methods of He et al. (2012) and Zhang et al. (2022) with minor modifications. 95% chlorantraniliprole (DuPont) was dissolved in acetone solution, and then diluted into five to six concentrations with distilled water containing .05% Triton X-100 at a triple gradient dose. The .05% Triton X-100 solution was used as blank control. The fresh leek cultured in the laboratory measuring 5 cm × .5 cm were immersed in the prepared various concentrations of chlorantraniliprole for 15 s, and taken out in ventilated cupboard to air dry and then placed into a Petri dish lined with filter paper. Twenty healthy 2-day third instar larvae were placed into each Petri dish with six pieces of leek leaves with four replications for toxicity assessment bioassays. Mortality was assessed after 48 h of exposure, larvae that did not move when gently pushed with a fine hair brush were considered dead. The control groups was below 5%, mortality was corrected using Abbott’s formula (Abbott, 1925). The LC₅₀ value and slope were calculated by regression probit analysis conducted with the POLO-Plus software (2002).

Sequences annotated as ryanodine receptor genes were screened from the transcriptome data published in Genbank (Bioproject No. 785282), and specific primers listed in Table 1 were designed for segmental cloning. The Total RNA was extracted from the whole body using RNAiso reagent (TaKaRa, Dalian, China) based on the manufacturer’s instructions, cDNA was synthesized by reverse transcription using the PrimeScript^® RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China) and as the PCR templates. Total PCR reaction system (50 μl) included cDNA template (2 μl), forward and reverse primers (1 μl, 0.2 μmol/L), amplification enzyme mix (25 μl), RNAse Free ddH₂O (21 μl). 2 × Vazyme LAmp^® Master Mix was purchased online from Vazyme Co., Ltd. and used for high fidelity and rapid PCR, the annealing velocity is high to 1 kb/s. PCR amplification conditions were performed as follows: 98°C for 3 min, followed by 35 cycles of 98°C for 10 s, 60°C for X s (X: each primer used a different annealing time), and 72°C for 3 min. Finally, it was extended for 5–10 min at 72°C (Table 1). The expected-size fragments were purified, ligated into pMD19-T vector (TaKaRa, code: D102A), and positive transformants were selected for plasmid isolation using MiniBEST Plasmid Purification Kit (TaKaRa, Dalian, China) and sent to Sangon Biotech company for sequencing. The fragments were overlapped and aligned with the annotated GdRyR gene identified from the transcriptome database using DNAMAN software, and the spliced complete CDS sequence of the GdRyR gene was submitted to the NCBI database, the Genbank accession number is OP828593.

ORF Finder was applied to search the open reading frame of the GdRyR gene, and sequence integrity was verified by NCBI’s Blastx tool to predict the gene isoelectric point, molecular mass, atomic composition, protein structural domain, signal peptide, and transmembrane region. Phylogenetic tree was constructed using the BLAST tool and DNAMAN V6.0 software to align amino acid sequence homology, using the neighbor-joining method in MEGA software and repeatedly running 1,000 times.

The total RNA was isolated from individuals of G. daurica population collected from the Xianghuang Banner grasslands for cDNA preparation and PCR. The two pairs of specific primers containing the potential mutation sites G4911 and I4754 were designed based on the sequence of the C-terminal transmembrane region (Table 1). The specific bands obtained by PCR amplification were purified and sent to company for sequencing, and the sequences from different individuals were compared in the point mutation region to detect the frequency of allelic mutation sites.

The amino acid sequence of GdRyR was compared with domains containing mutation loci for P. xylostella RyR (G4946E, I4790M), S. frugiperda for (G4900E, I4743M) and S. exigua (G4891E, I4743M) (Zou et al., 2019), which correspond to G4900E, I4754M mutation sites of GdRyR. Mutation of “Ile” was replaced by “Met” in I4754 in GdRyR using AlphaFold2, the mutation model GdRyR-C-M4754 was constructed. Similarly, “Gly” was substituted into “Glu” in G4900 in GdRyR, and the mutation model GdRyR-C-E4900 was constructed. The models WT GdRyR-C, GdRyR-C-M4754 and GdRyR-C-E4911were processed separately using Autodock Tool-1.5.7 for energy minimization and dehydrogenation as molecular docking receptors.

The 3D structures of chlorantraniliprole and cyanobacteriamide were downloaded from the PubChem website, and converted into a suitable format for molecular docking using Open Babel-2.4.1, and used as ligands for molecular docking using Autodock Tool-1.5.6 for energy minimization etc.

Based on the reported binding sites of the rabbit RyR complex crystal structure with chlorantraniliprole, the active pockets of the constructed 3D model were predicted and the docking-box was constructed. The box parameters used for the docking of GdRyRs were as follows: center-x = 8.18, center-y = −1.56, center-z = .877; size-x = 99.75, size-y = 99.75, size-z = 99.75. The processed protein receptors were docked to the ligands using Autodock vina1.1.2 software. Computing platform was as follows: Microsoft-PC: Intel(R) Core(TM) i5-1035G4CPU@1.10GHz 1.50 GHz.

By molecular docking, the two ligands were selected for optimal conformation to construct complexes with the ryanodine receptor, and the binding modes were analyzed using the PLIP online tool, including: hydrogen bonding, hydrophobic interactions, π-π, π-cation, halogen bonding and non-covalent interaction forces such as interactions. Pymol-2.1.0 was used to visualize the binding modes.

Homologous sequence alignment: https://blast.ncbi.nlm.nih.gov/Blast.cgi; Predicting protein properties such as gene isoelectric point, molecular mass and atomic composition: ProtParam online tool; Protein domain prediction: https://www.ebi.ac.uk/interpro/result/InterProScan/(InterProscan); Signal peptide prediction: http://www.cbs.dtu.dk/services/SigRyRlP/; Homologous modeling: https://swissmodel.expasy.org/(SWISS-MODEL) (Waterhouse et al., 2018); Protein quality evaluation: https://Main page-MolProbity (duke.edu) (Molprobity) (Williams et al., 2018); Protein crystal structure PDB data: http://www1.rcsb.org/; NCBI databse: https://www.ncbi.nlm.nih.gov/; Transmembrane region prediction (TMHMM): http://www.cbs.dtu.dk/services/TMHMM/(Möller et al., 2001); Action force analysis tools (PLIP): https://projects.biotec.tu-dresden.de/plip-web/plip (Melissa et al., 2021); Small molecular crystal structure: PubChem (nih.gov ).

MEGA X was applied to construct phylogenetic tree (Kumar et al., 2018); Open Babel-2.4.1 was used to convert the file format of the chemical structure type (O’Boyle et al., 2011); The 3D protein crystal structure was visualized by Pymol-2.1.0 (Delano, 2002); Modeller (Version-10.2) was used for homologous modeling (Chen et al., 2021); Autodock Tool-1.5.7 (Morris et al., 2009) and Autodock vina 1.1.2 (Trott and Olson., 2010) were used for molecular docking; AlphaFold2 was applied for deep learning modeling (Jumper et al., 2021).

Linear regression of the dose–mortality relationship (Y = −2.557 ± 2.291 X) was fitted to the observed data (i.e. no significant deviation between the observed and the expected data; χ ² = 1.81, df = 3) and LC₅₀ was considered valid (inf lim < LC₅₀ < sup lim: 21.25 < 44.10 mg/L < 80.33).

Based on the GdRyR sequence information in the transcriptome, a total of four target fragments were amplified. After overlapped, a complete open reading frame (ORF) of 15,399 bp, encoding 5,133 amino acids, was obtained by BlastX alignment and sequence splicing. The molecular weight of the protein is 582.53 kDa, and the isoelectric point (pI) is 5.51. The GdRyR consists of five elements: carbon (C), hydrogen (H), Nitrogen (N), Oxygen (O) and Sulphur (S), and the total number of atoms is 81,490, the chemical formula is C₂₅₈₉₂₅H₄₀₅₃₈N₆₉₉₂O₇₇₉₉S_236, it has six transmembrane structure regions and no signal peptide. As shown in Figure 1, the N-terminal part of the gene contains one MIR domain (Mannosyltransferase, IP3R and RyR) located between amino acid sites 220-397, two RIH domains (RyR and IP3R Homology) located at amino acid sites 447-642 and 2254-2484, respectively, and three SPRY domains (SPla and RyR) located at amino acid sites 664-802, 1091-1214 and 1563-1706, respectively. The gene also has four RyR domains located at amino acid sites 854-944, 967-1056, 2861-2952, 2979-3062, and one highly conserved RIH-associated domain located at amino acid sites 4021–4137 before the transmembrane helix; The C-terminus has six transmembrane helices (TM1 to TM6) located at amino acid sites 4474-4496, 4658-4680, 4739-4761, 4881-4903, 4929-4951, and 5009-5028, respectively.

The amino acid sequences encoded by the GdRyR gene were blasted with RyR of other insect species researched in the NCBI database, and a phylogenetic tree was constructed using the neighbour-joining (NJ) method. Figure 2 showed that GdRyR has the highest homology with RyR of potato beetle Leptinotarsa decemlineata (GenBank Accession No: QZZ63290.1), and amino acid sequence identity is 91.33%. The amino acid homology with Hypothenemus hampei RyR (QEE14187.1) and Tribolium castaneum RyR (NP 001308588.1) were 88.05% and 86.70%, respectively, and the RyR of the four Coleopteran species clustered as one clade. The clustering results indicated that the GdRyR gene is highly similar to taxonomically similar insect taxa.

The two potential allelic mutation sites and the occurrence frequency of Gly4911Glu and Ile4754Met were examined by individual extraction of RNA, using reverse-transcribed cDNA as a template, respectively. Fifty beetles collected from the grasslands of Xilingol Xianghuang Banner, Inner Mongolia, were conducted to gene sequence testing, the mutation loci and allelic frequencies of the GdRyR are summarized in Table 2. The results showed that only the heterozygous genotype G4911E was present in the test population, accounting for approximately 12% of the population; the heterozygous and homozygous genotypes of I4754M were present in 32% and 2% of the population, respectively.

The Modeller and AlphaFold2 are applied to construct the GdRyR 3D crystal structure (WT GdRyR), the C-terminal wild-type 3D structure of GdRyR (WT GdRyR-C) and the mutant structure (GdRyR-C-M4754) (GdRyR-C-E4911), respectively. Figure 3A showed the 3D structure of GdRyR consists of multiple α-helices, β-folds, β-turns, irregular curls, and some extended extension structures. In lateral view, the protein crystal structure resembles the letter “Y,” with a wide N-terminus and a narrow C-terminus, separated by a small angle at the C-terminus. As it can be seen from Figures 3C, E, G that the 3D structure of the C-terminal region of GdRyR is mainly composed of several parallel sets of α-helices, exhibiting obvious transmembrane protein properties. The 3D structural model of Ramachandran Plot diagram and the test results are shown in Figure 3 and Table 3, respectively. The percentages of residues in the optimal regions for WT GdRyR, WT GdRyR-C and GdRyR-C-M4754 and GdRyR-C-E4911 were 86.9%, 88.3%, 87.8%, and 88.1%, respectively. The proportions of amino acid residues in the acceptable region were 95.8%, 94.6%, 94.9%, and 94.6% respectively, all are greater than 90%, indicating that the constructed models were all reasonable.

The binding affinities of WT GdRyR-C to the optimal conformation of chlorantraniliprole and cyantraniliprole were −31.38, −33. 89 kJ/mol, respectively. The binding pattern was analyzed using PLIP and drawn by Pymol. As Figure 4 shows: Y4660 (Tyr-220), K4663 (Lys-223) form a hydrophobic interaction with the benzene and pyrazole rings of chlorantraniliprole, K4663 forms a hydrogen bond with the oxygen atom on the carbonyl group, K4762 (Lys-322) forms a hydrogen bond with the nitrogen atom of the pyrazole ring, G4911 (Gly-471) forms a hydrogen bond with the NH of the methylamino structure on the benzene ring, R4916 (Arg-476) forms a hydrogen bond with the nitrogen atom of the pyridine ring (Figure 4A). Y4660 forms a hydrophobic interaction with the methyl group on the benzene ring of cyantraniliprole, K4762 forms a hydrogen bond with the nitrogen atom of the pyrazole ring, G4911 forms a hydrogen bond with the NH of the amide structure on the side of the pyrazole ring, the nitrogen atom of R4916 forms a π-cation with the pyrazole ring of cyantraniliprole (Figure 4B).

The binding affinities of GdRyR-M4754 to the optimal conformation of chlorantraniliprole and cyantraniliprole were −32.23, −33. 91 kJ/mol, respectively. Combined pattern analysis diagram was seen in Figure 5: Y4660 (Tyr-220) and K4663 (Lys-223) form hydrophobic interactions with the pyridine ring of chlorantraniliprole, K4913 (Lys-473) forms hydrophobic interactions with the benzene ring as well as the methyl group, Y4660 forms a hydrogen bond with the oxygen atom on the carbonyl group of the benzene ring, K4663 forms a hydrogen bond with the nitrogen atom of the pyrazole ring, R4916 (Arg-476) forms a hydrogen bond with the NH of the amide structure on the pyrazole side (Figure 5A). Y4660, V4910 (Val-470) and K4663 form hydrophobic interactions with the pyridine ring, pyrazole ring and methyl group of cyantraniliprole, respectively. K4663 forms a hydrogen bond with nitrogen on the pyridine ring, K4762 (Lys-332) forms a hydrogen bond with the nitrogen atom of the pyrazole ring, G4911 (Gly-471) forms a hydrogen bond with NH of the amide structure on the benzene ring side, R4916 forms a hydrogen bond with the oxygen atom of the carbonyl group on the pyrazole ring side, the nitrogen atom of K4663 forms a π-cation interaction with the pyridine ring of cyantraniliprole (Figure 5B).

The binding affinities of GdRyR-E4911 to the optimal conformation of chlorantraniliprole and cyantraniliprole were −25.89, −27.54 kJ/mol, respectively. Combined pattern analysis diagram was seen in Figure 6: K4762 (Lys-329) forms hydrophobic interactions with methyl of chlorantraniliprole, K4769 (Lys-322) and R4770 (Arg-330) form hydrophobic interactions with the pyridine ring and the chlorine ion on the pyridine ring, respectively. D4878 (Asp-438) forms hydrophobic interactions and hydrogen bonds with the pyrazole ring and the NH of the pyrazole side amide structure, respectively (Figure 6A). L4497 (Leu-60) and L4500 (Leu-57) form hydrophobic interactions with the methyl and benzene rings of cyantraniliprole, respectively. M4501 (Met-61) and Q4868 (Gln-428) form hydrophobic interactions with the pyridine ring, Y4759 (Tyr-341) and F4881 (Phe-441) form hydrophobic interactions with the benzene ring, Y4759 forms a hydrogen bond with the NH of the amide structure on the side of the benzene ring, Q4868 (Gln-428) forms a hydrogen bond with the nitrogen atom of the pyrazole ring (Figure 6B).

As shown in Table 4: The mutation of isoleucine (Ile) to methionine (Met) at position 4754 in GdRyR showed no decrease in binding affinity to diamide insecticides; the mutation of glycine (Gly) to glutamic acid (Glu) at position 4911 in GdRyR showed a significant decrease in binding affinity to diamide insecticides. The binding modes of two diamide insecticides with WT GdRyR, GDRYR-M4754 and GDRYR-E4911 are shown in Figures 4–6: After mutating isoleucine (Ile, position 4754) of GdRyR to methionine (Met), the mode of interaction with diamide insecticides, the residues involved in the formation of the action force and WT GdRyR are essentially the same, mainly because the isoleucine (Ile) at position 4754 of WT GdRyR is not involved in the formation of the force action between ligands and receptors. The mutation to glutamic acid (Glu) in WT GdRyR, which has a hydrogen bonding interaction between glycine (Gly) at position 4911 and the diamide insecticide, prevents the formation of hydrogen bonds with the diamide insecticide, resulting in a reduced receptor-ligand interaction and decreased affinity. The mutations of amino acid from Gly to Glu at position 4911 in WT GdRyR may lead to resistance to diamide insecticides on G. daurica, while further studies are needed to determine whether the I4734M mutation is associated with resistance.