Transcriptome and Micro-CT analysis unravels the cuticle modification in phosphine-resistant stored grain insect pest, Tribolium castaneum (Herbst)

Phosphine (PH3) resistance in stored grain insect pests poses a significant challenge to effective pest control strategies worldwide. This study delved into understanding PH3-resistant mechanism, with the objective of informing robust and sustainable pest management strategies that could mitigate the impacts of PH3 resistance. In this regard, the transcriptomic analysis identified 23 genes associated with chitin synthesis and cuticle formation, which showed significant expression in PH3-resistant (R) strains compared to susceptible strains. Micro-computed tomography (Micro-CT) revealed an extended and tighter cuticular structure in the PH3-R Tribolium castaneum than PH3-susceptible strains but with no changes in the cuticle thickness. This altered cuticle structure may reduce PH3 penetration through cuticles rather than completely closing spiracles during fumigation. It is also hypothesized to prevent water loss from the insect body, as water production decreased in PH3-R T. castaneum due to the down-regulation of the electron transport chain function. Validation of several chitin synthesis gene expression levels revealed consistent results with those of transcriptomic analysis. Overall, integrating physical treatments using synthetic amorphous silicates, water absorbents, and cuticle-damaging materials during PH3 fumigation is recommended for its prolonged and controlled usage in the field.


Background
Pesticides and fertilizers play an instrumental role in bolstering agricultural cereal production [44,49].Their absence would likely trigger severe food supply issues across the agricultural industry [17,43].Unlike fertilizers, pesticides play a part in every stage of agriculture, from production and harvesting to storage and distribution, requiring diverse formulations for each process.Importantly, these substances must not linger in the final products, and their residual toxicity should not surpass the maximum residue levels established by individual countries [9].All consumer-bound agricultural products must strive to be free of pesticide residues while preserving the quality of these goods [13,23].
Phosphine (PH 3 ), a gaseous fumigant, is the most important fumigant for controlling insect pest not only in worldwide quarantine as an alternative to methyl bromide, but also in stored grain storages.Enormous amounts of PH 3 have been used due to their strong penetration ability and high efficacy without residue issues in grains [10,29].However, because of this, the resistance of many stored grain insect pests to PH 3 is causing a large loss of grains during the storage of agricultural products, and it has reached the point where it is considered whether to use PH 3 depending on the degree of resistance to PH 3 [29].Mechanisms of PH 3 resistance caused by insect pests have been introduced in association with changes in the electron transport chain and changes in the dihydrolipoamide dehydrogenase (dld) gene [21,25,35].In particular, many researchers claim that strong phosphine resistance is mainly related to changes in the dld gene [21,35].However, it is questionable whether PH 3 resistance is limited to changes in these metabolismrelated genes.
Notably, the use of pesticides causes two major problems: first, the development of resistance to pesticides Graphical abstract by the target organisms, which leads to increased pesticide use and their unavoidable necessity [46].Second, the excessive use of pesticides leads to environmental pollution, which results from the temporary increase in their levels entering the environment and the effect of pesticides combined with various matrices, which can cause acute and persistent problems [36].Numerous studies have reported the underlying resistance mechanisms associated with pesticide-induced resistance, such as overexpression of detoxifying enzymes, changes in target sites, and the reduction of pesticide penetration rate into the target site [7,14,16,48].Typical examples of overexpressed detoxifying enzymes include cytochrome P450s monooxygenases (P450s) and carboxylesterases [7,14].The overexpression of the enzymes has been demonstrated in several insect pests and a variety of pesticides are associated with them.P450s influence the addition of oxygen to xenobiotics.For example, they alter the phosphorothioate (P = S) moiety to the P = O form.Changes in the target sites are associated with modifications in sodium channels and acetylcholinesterases [1,50].The changes have been reported to reduce pesticidal activities, and even though the structures are altered, the original activity is maintained [1].Contact insecticides can pass through the cuticle layer and a reduced penetration rate has been observed for insecticide resistance measured based on the absorption rate of radiolabeled pesticides in an insect body, as well as the structural changes of the cuticle layers using transmission scanning electron microscope (TEM) and scanning electron microscope (SEM) [7].
PH 3 resistance mechanisms in insect pests differ from those caused by contact insecticides.Resistance mechanisms to PH 3 among insect pests are strongly associated with changes in the electron transport chain (ETC) and dld gene, the supplier of reducing equivalents to the ETC [25,34,35].In particular, several researchers have claimed that strong PH 3 resistance in stored product insect pests is primarily associated with changes in the dld expression [5,35].However, whether PH 3 resistance is limited to changes in such metabolism-related genes remains unknown because fumigant toxicity may occur under anoxic conditions rather than an aerobic environment.
Insects maintain a specific respiratory system through which O 2 and CO 2 gases are exchanged, which is referred to as the discontinuous gas exchange cycle (DGC).The insect respiratory system consists of the following phases based on the spiracle behavior: closed (C-phase), fluttered (F-phase), and open (O-phase) phases to avoid oxygen toxicity [20].Therefore, insects resistant to PH 3 treatment may develop different types of mechanisms to overcome its toxicity under anoxic conditions when they are in the C-phase.To effectively verify that PH 3 resistance is associated with changes in the cuticle layer [4], the following hypothesis has been proposed.If a certain PH 3 concentration passes through the insect respiratory system, the insect will close all spiracles or the tracheal system, and a low PH 3 concentration passes through the cuticle layer into the insect body.Subsequently, the insect cuticle structure is altered, in turn, reducing the amount of PH 3 passing through the cuticle layer.
This study used transcriptomic analysis to identify differentially expressed genes between the PH 3 -susceptible and -resistant strains of the red flour beetle, Tribolium castaneum (Herbst).After comparing the transcriptomic data, structural differences between the cuticles of PH 3 -susceptible and PH 3 -resistant strains of T. castaneum adults were examined using Micro-computed tomography (Micro-CT).In addition, changes in the cuticle layers were further analyzed to determine variations in the expression of genes associated with cuticle formation.The analysis was performed to determine the possibility of accurately describing all changes observed in different cuticle parts using single transmission or scanning electron microscope (TEM or SEM) sections as representative samples.

Insect and growth conditions
Phosphine-susceptible (PH 3 -S, Aus-10) and phosphineresistant (PH 3 -R, Aus-07) strains of T. castaneum were provided by the Plant Quarantine Technology Center of the Animal and Plant Quarantine Agency (Gimcheon, Republic of Korea).Both strains were maintained under fumigant-free conditions.All T. castaneum stock colonies were successively cultured on whole wheat flour with yeast extract (20:1, g/g) at the Kyungpook National University (Daegu, Republic of Korea) under controlled conditions (temperature of 28℃ ± 1℃, relative humidity of 50% ± 2%, and a photoperiod of 16 h light/8 h dark).All experiments except the Micro-CT scanning were conducted using T. castaneum adults selected randomly from over 4 months old.The Micro-CT scanning experiments were conducted with females and males separately.

PH 3 fumigation bioassay and genotyping of P45S mutation on dld
To compare PH 3 resistance between PH 3 -S and PH 3 -R strains, bioassays were conducted using the PH 3 fumigant and dld point mutation (P45S) was identified.PH 3 fumigation was performed in a 12 L-desiccator at 20 °C under normal photoperiod conditions (16 h light/8 h dark) according to the safety standards of the Plant Quarantine Technology Center of the Animal and Plant Quarantine Agency (Gimcheon, Republic of Korea).T. castaneum adults (n = 30) from each strain were exposed to PH 3 concentrations of 0 (control, CON), 0.005, 0.008, 0.013, 2, 2.4, and 2.75 g/m 3 for 20 h based on the PH 3 quarantine standards.Furthermore, mortality was determined at three days post-fumigation (dpf ) and the lethal concentration (LC 50 ) value was calculated by probit analysis using SAS software version 9.4 (SAS Institute Inc., Cary, NC, USA).To calculate the lethal concentration time (LCT 50 ) value, gas sampling was performed at 0.5, 1, 2, 4, 18, and 20 h during fumigation using 1-L Tedlar bags (SKC Inc., Dorset, UK) [25].PH 3 concentration was determined by gas chromatography (GC) using a flame photometric detector (Agilent 7890 B GC-FPD, Agilent Technologies Inc., Santa Clara, CA, USA) equipped with an HP-PLOT/Q capillary column (30.0 m × 530 µm × 40.0 µm; Agilent Technologies Inc., Santa Clara, CA, USA).The standard curve for PH 3 was calculated in the 0.001-3.0g/m 3 range.The P45S point mutation on dld, a diagnostic molecular marker of PH 3 resistance, was identified from each strain (n = 5) in triplicate using the cleaved amplified polymorphic sequence (CAPS) method [11].

Transcriptomic analysis using DAVID
Total RNA samples were extracted from 10 individuals in the adult stage that were selected from Aus-10 and Aus-07 strains, in triplicate.To remove flour, the collected samples were placed in microtubes and washed more than twice with DEPC-treated water, followed by homogenization with 1 mL TRIzol reagent (Thermo Fisher Scientific Inc., Waltham, MA, USA), according to the manufacturer's instructions.The total RNA concentration was calculated using Quant-iT RiboGreen RNA reagent (Invitrogen, Carlsbad, CA, USA).Only high-quality RNA preparations with RNA integrity number (RIN > 7) were used for RNA library construction.A library was prepared from 1 µg of total RNA of each sample using the TruSeq Stranded Total RNA LT Sample Prep Kit (Gold; Illumina Inc., San Diego, CA, USA).The libraries were quantified by qPCR according to the qPCR Quantification Protocol Guide (KAPA Library Quantification Kits for Illumina sequencing platforms) and qualified using the TapeStation D1000 ScreenTape (Agilent Technologies, Waldbronn, Germany).Indexed libraries were then sequenced using the paired-end sequencing method (100-nt) on the Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA).The sequenced reads were mapped to the reference genome of T. castaneum using Bowtie2 and HISAT2 programs.Differentially expressed genes (DEGs) were extracted based on log 2 fold-change and a p-value < 0.05.Gene functional annotation analysis for the DEGs was performed using the DAVID tool (http:// david.abcc.ncifc rf.gov/) to understand their biological meanings.The DEG data were analyzed using R software and visualized using principal component analysis (PCA) plots, heatmaps, and volcano plots.The "ggplot2" package in R Studio (www.r-project.org) was used to visualize the DEG data.

Sample preparation for Micro-CT analysis
T. castaneum adults, which were being reared in breeding dishes, were separated into females and males and fixed in ethanol (99.5% purity) overnight.Afterward, the samples were stained overnight with I2E solution (1% iodine dissolved in 100% ethanol) and then washed with ethanol [27].Subsequently, the lower part of a 10-µL pipette tip (polypropylene) was sealed by heating and filled with ethanol.Two T. castaneum adults were placed inside the sealed pipette tip with their heads facing upward and cotton wool was used to prevent contact between the individuals [41] and the scheme of sample preparation is shown in Fig. 1a.Gas inside T. castaneum can form bubbles as it escapes, in turn, causing blurring during the reconstruction process.Therefore, the bubbles were removed and the upper part of the pipette tip was sealed with parafilm during the reconstruction process.

Micro-CT scanning and reconstruction
The 10-µL tip containing the whole body of T. castaneum was scanned by Micro-CT (Skyscan 1272, Kontich, Belgium), and all samples were scanned under the conditions presented in Additional file 1: Table S1.The images obtained were reconstructed and converted into crosssectional images using NRecon v.1.7.1.0(Skyscan, Bruker Micro-CT, Belgium).The region of interest (ROI) was set to include the whole insect body.The threshold in the density histogram was set from 0.032963 to 0.071454 (Additional file 1: Table S1 and Fig. 1b).

Post-processing for extraction of the cuticle layer
To analyze only the cuticle layer, which is the hardest part of the organism, the grayscale indexes in the reconstructed cross-sectional image datasets of T. castaneum obtained by Micro-CT scanning were set from 254 to 255 using a CTAN v.1.18.8.0 (Bruker Micro-CT, Belgium).
To analyze the same part of the sample, only 371 slice images in the posterior direction from the first abdominal sternite were used as a reference for all individuals.In addition, ROI was set accordingly to exclude the insect legs that could interfere with the analysis.To remove various noises from the body cavity of the sample, a new ROI was set and white speckles smaller than 50 pixels were removed.ROI was set to exclude the cotton wool surrounding the cuticle layer, and each procedure was visualized using CTvox v.3.3.0 (Bruker Micro-CT, Belgium) and CTAN (Fig. 1c).Morphometric measurements were conducted on each processed cuticle layer dataset.The selected parameters were analyzed according to the manufacturer's protocol: mean total cross-sectional object area (2D Obj.Ar, µm 2 ), mean eccentricity (2D ECC), object surface area to volume ratio (2D Obj.S/Obj.Reconstructed cross-sectional images were thresholded with a grayscale index (254-255) (binarization).Cutting processing was only performed for 371 slides from first abdominal sterna to the posterior direction and the parts corresponding to the legs were cut (cutting).Removing noise from the inside (removal of inside noise) and outside (removal of outside noise) of the cuticle layer was followed V, 1/µm) obtained by dividing 2D object surface area (2D Obj.S, µm 2 ) by 2D object volume (2D Obj.V, µm 3 ), mean number of objects per slice (2D Obj.N), average object area per slice (2D Av.Obj.Ar, µm 2 ), closed porosity (2D Po(cl), %) in 2D analysis, and percent object volume (3D Obj.V/TV, %) obtained by dividing 3D Obj.V by 3D total volume (3D TV, µm 3 ), object surface area (Obj.S, µm 2 ), object surface area to volume ratio (3D Obj.S/Obj.V, 1/ µm), structure thickness (3D St. Th, µm), and structure model index (3D SMI) in 3D analysis.The mentioned parameters have been grouped into basic parameters, cuticle structure and thickness parameters, and connectivity parameters.Basic parameters consist of 3D Obj.V/TV, 3D Obj.S, and 2D Obj.Ar, serving the purpose of confirming the segmentation of the cuticle layer of T. castaneum with error and loss.Cuticle structure and thickness parameters comprise 2D ECC, 2D Obj.S/ Obj.V, 3D obj.V, and 3D St.Th, designed for comparing the properties of the cuticle and plate thickness.Connectivity parameters, on the other hand, are composed of 2D Obj.N, 2D Av.Obj.Ar, 3D SMI, and 2D Po(cl), serving the purpose of comparing the continuity or tightness of the cuticle structure.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR)
Total RNA was extracted from 10 individuals in the adult stage that were selected from Aus-10 and Aus-07 strains, respectively, in triplicate.The extraction was performed using the same methods described in "Transcriptomic analysis using DAVID" section.The absorbance of the total RNA extracted was measured at 260/280 nm using a μDrop plate system (Thermo Fisher Scientific, Waltham, MA, USA), and the concentration of RNA was normalized to 50 ng/µL.Complementary DNA (cDNA) was synthesized using RNA and Maxima First Strand cDNA Synthesis Kit with dsDNase (Thermo Fisher Scientific, Waltham, MA, USA) and all samples were diluted to the same concentration.Luna Universal qPCR Master Mix (New England Biolabs, Ipswich, MA, USA) was used for qRT-PCR and verification was carried out using the QuantStudio 3 Real-Time PCR system (Applied Biosystems, Foster City, CA, USA).The experimental procedures were performed three times.The primers used to amplify genes associated with chitin biosynthesis are listed in Additional file 2: Table S2.

Measurement of chitin contents
The procedure of cuticle content measurement using Calcofluor (F3543, Sigma-Aldrich Co., St. Louis, MO) was slightly modified following the reference [19].Briefly, adults of T. castaneum were randomly selected from PH 3 -S or PH 3 -R strain group and starved for 5 days before the experiment.To measure chitin contents, colloidal chitin suspension was prepared using 0.5 g of adult T. castaneum (approximately 250 individuals) and commercial chitin from shrimp shells (C7170, Sigma-Aldrich Co., St. Louis, MO) for the standard curve.Samples were homogenized using a mortar and pestle with liquid nitrogen and were dissolved in 10 mL hydrochloric acid (37% HCl, extra pure grade, Duksan Pure Chemical Company, Ansan, Republic of Korea) with stirring for 50 min at room temperature (RT).Then, 25 mL of distilled water (DW) was added for the precipitation of chitin.The colloidal suspension was centrifuged at 3500 rpm for 5 min, and the pellet of chitin was washed with DW until the pH of the solution was above 3.5.Finally, the pellet was resuspended with 100 mL DW.For the standard curve, 5 µg/µL of colloidal chitin suspension from commercial chitin was used for the stock solution, and a standard curve was plotted at the range of 0, 1, 2.5, 5, 10, 25, and 30 µg chitin.Each colloidal chitin suspension was incubated with 100 µL of Calcofluor solution at the dark condition for 15 min, then centrifuged at 13,000 rpm for 5 min.Pellets were washed with DW for two times and resuspended with 200 µL of DW in a black 96-well plate.The fluorescent measurement was conducted using Spec-traMax iD3 (Molecular Devices, San Jose, CA) with excitation at 350 nm and emission at 430 nm.

Statistical analysis
Data are expressed as means ± standard deviation.Statistically significant differences between PH 3 -S and PH 3 -R strains were evaluated using a t-test.One-way ANOVA followed by Tukey's test for multiple comparisons was performed to evaluate statistically significant differences between susceptible and resistant strains of T. castaneum based on sex.All statistical analyses were performed using Graphpad Prism 9.0 (Graphpad Software Inc., San Diego, CA, USA) and the means were considered significant at p < 0.05.

Transcriptomic analysis revealed the relationship between PH 3 resistance and chitin-related biological process
The LC 50 and LCT 50 values of the PH 3 -R strain were 235-fold and 96.3-fold, respectively, higher than those of the PH 3 -S strain (Fig. 2a and Table 1).The diagnostic gene marker cleavage, a P45S point mutation on dld, was observed in the PH 3 -R strain but not in the PH 3 -S strain (Fig. 2c).Furthermore, the homogeneous pool of P45S point mutations was dominant in the PH 3 -R strain, accounting for over 60%.In the PCA plot with 12,898 identified transcripts in T. castaneum, two groups were distinctly separated, with 82.62% in the first principal component (PC1) and 10.4% in the second principal  component (PC2) (Fig. 2d).A total of 508 transcripts were selected based on a p-value < 0.05, and upregulated and downregulated gene expressions were grouped by hierarchical clustering based on the three replicates for each group (Fig. 2e and Additional file 2: Table S3).
To determine significant differences between the PH 3 -S and PH 3 -R strains based on the DEGs, functional annotation analyses were conducted and significant gene ontology (GO) terms at p < 0.05 were divided into four categories: cellular component (CC), molecular function (MF), biological process (BP), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (Fig. 3a).Notably, chitin-related terms were ranked in every category: chitin binding (p < 0.01) in MF, chitin metabolic process (p < 0.01) and chitin catabolic process (p < 0.05) in BP, and amino sugar and nucleotide sugar metabolism (p < 0.05) in the KEGG pathway (Fig. 3a).In addition, 23 genes associated with chitin synthesis and cuticle  S3 formation were selectively represented in the volcano plot, including genes that were functionally annotated by DAVID (Fig. 3b and Additional file 2: Table S4).Chitin-related genes in the PH 3 -R strain were significantly upregulated when compared to those in the PH 3 -S strain (Fig. 3b).Chitinase 5, 10, and 13 (cht5, cht10, and cht13), which are associated with chitin synthesis, were significantly upregulated (p < 0.01) by 5.2, 6.8, and 8.6-fold, respectively, in the PH 3 -R strain, indicating that the chitin catabolic process in the PH 3 -R strain was more activated than in the PH 3 -S strain (Fig. 3b).A significantly high level (p < 0.05) of the facilitated trehalose transporter Tret1 (tret1) gene expression was observed on the right side of the volcano plot, indicating that trehalose transporters, which are associated with the regulation of trehalose levels, were more highly activated in the PH 3 -R strain than in the PH 3 -S strain.Trehalose is a precursor of chitin [37].Moreover, the expression of genes encoding peritrophic matrix proteins 1-A, 1-B, and 1-C (pmp1-a, pmp1-b, and pmp1-c) was consistently upregulated in the PH 3 -R strain with gene counts of the chitin-binding term in MF (Fig. 3b and Additional file 2: Table S4).Overall, the PH 3 -R strain exhibited a remarkably distinct gene expression pattern in chitin metabolic and catabolic processes.Therefore, we speculate that the cuticle layer in the PH 3 -R strain was influenced by the significant differences in the expression patterns of chitin-related genes.

PH 3 -induced morphometric changes in the cuticle of T. castaneum Basic parameters
No significant differences were observed in the cuticle appearance in the 3D modeling of T. castaneum based on strain and sex (Fig. 4).The results of morphometric analysis conducted on the segmented cuticle layers of both PH 3 -S and PH 3 -R (Union) strains and female and male (Subset) T. castaneum are presented in Fig. 4. In the basic 3D morphometric parameters, the total volume of interest (TV, µm 3 ) represented the total volume of the space used to analyze the cuticle layer (Fig. 5a, black space and white space) and object volume (Obj.V, µm 3 ) was the volume of the binarized cuticle layer within the TV (Fig. 5a, white space).The percent object volume in 3D analysis (3D Obj.V/TV, %) represented the volume of the cuticle layer as a percentage of the total volume.The value was calculated by summing the continuous 2D cross-sectional images of T. castaneum.The object surface area in 3D analysis (3D Obj.S, µm 2 ) was the total surface area of the cuticle layer occupying the 3D space (Fig. 5, white surface).The mean total cross-sectional object area in 2D Fig. 4. 3D modeling of PH 3 -susceptible (S) and resistant (R) Tribolium castaneum whole body and segmentation of the conserved cuticle layer.In 3D modeling of whole body, each conserved cuticle layer (red box) was segmented and analyzed.SF, susceptible female; RF, resistant female; SM, susceptible male; RM, resistant male analysis (Obj.Ar µm 2 ) was the mean area occupied by the cuticle layer in the cross-sectional image in the 2D plane (Fig. 5c).Regarding the basic parameters, no significant differences were observed between PH 3 -S and PH 3 -R strains in the Union group and among Subset groups.The results indicate that all parts of the cuticle layer of T. castaneum were segmented with minimal error and no loss.

Cuticle structure and thickness parameters
The mean eccentricity in 2D analysis (2D ECC) was the average eccentricity of the shape of the cuticle layer in the segmented cross-sectional image of T. castaneum.ECC increased as the shape of the cuticle layer became more elliptical, with a perfect sphere having a value of 0 (Additional file 1: Fig. S1a).The ratio of solid surface to volume in 2D analysis (2D Obj.S/Obj.V, 1/µm) and in 3D analysis (3D Obj.S/Obj.V, 1/µm) represented the uniformity or unevenness of the surface area of the cuticle layer in relation to its volume, with lower values indicating greater uniformity and higher values indicating greater unevenness (Additional file 1: Fig. S1b and c).Based on the structural parameters, no statistically significant differences were observed between the union and subset groups.Thickness in 3D analysis (3D St. Th, µm) refers to the average thickness of the continuous cuticle layer and Fig. 5 Basic structural parameters of PH 3 -susceptible (S) and resistant (R) strain of Tribolium castaneum using Micro-CT analysis.Data were obtained by analyzing the same part, divided into PH 3 -S and R (Union), and male and female (Subset).Schematic illustrations are representing for a percent object volume (3D Obj.V/TV, %) obtained by dividing 3D object volume (3D Obj.V, µm 3 ) by 3D total volume (TV, µm 3 ), b object surface (3D Obj.S, µm 2 ) and c mean total cross-sectional object area (2D Obj.Ar, µm 2 ).All data are presented as mean ± standard deviation.Union group was analyzed using t-test and subset group was analyzed using one-way ANOVA, followed by Tukey's multiple comparisons.Asterisk represents the statistic mark in Union group (*p < 0.05; **p < 0.01; ***p < 0.001).SF, susceptible female; RF, resistant female; SM, susceptible male; RM, resistant male is a parameter equivalent to trabecular thickness in bone research (Additional file 1: Fig. S1d).No significant differences were observed in the cuticle thickness between the Union and Subset groups.Overall, the results indicate that all segmented cuticle layers in PH 3 -S and PH 3 -R strains of T. castaneum have similar properties and continuous plate thickness.

Connectivity parameters
Considerable differences were observed in connectivity parameters between the Union and Subset groups.The mean number of objects per slice in 2D analysis (2D Obj.N) indicated the mean number of cuticle parts distinguished in a single cross-sectional image of T. castaneum in 2D analysis.The number of cuticle parts decreased and exhibited a low value when the parts were merged, while the number of cuticle parts increased and exhibited a high value when the parts were disconnected (Fig. 6a).The average object area per slice in 2D analysis (Av.Obj.Ar, µm 2 ) is the mean area of cuticle parts distinguished in a single cross-sectional image of T. castaneum in 2D analysis.Unlike the "2D Obj.N" parameter, the "Av.Obj.Ar" parameter had higher value when the cuticle parts were interconnected and merged, and lower value when they were disconnected (Fig. 6b).These two parameters facilitate the comparison of connectivity levels within the interconnected structure of cuticle layers.The "2D Obj.N" parameter showed that the PH 3 -S strain had a higher average number of cuticle parts distinguished in a single cross-sectional image in 2D analysis than the PH 3 -R strain (p < 0.001), while the "2D Av.Obj.Ar" parameter showed that the PH 3 -S strain had a lower average area of cuticle parts distinguished in a single cross-sectional image in 2D analysis than the PH 3 -R strain (p < 0.001).The structure model index in 3D analysis (3D SMI) indicated a dominant shape (plate-like or cylindrical) in the 3D structure.The ideal plate-like, cylindrical, and spherical shapes had 0, 3, and 4 values, respectively, while cylindrical and spherical cavities had − 3 and − 4 values, respectively (Fig. 6c).The plate-like structure became dominant and the value decreased as the structure became more densely packed, while the cylindrical structure became more dominant and the value increased as the empty spaces in the structure increased.The PH 3 -S strain had a significantly higher SMI value than the PH 3 -R strain (p < 0.05).The "closed porosity (percent) in 2D analysis" (Po (cl), %) of the 2D cross-sectional image of the cuticle layer in T. castaneum was determined.Processed cuticle layers were not completely filled with cuticle layer parts, but rather consisted of enclosed spaces.Enclosed pores that were completely surrounded by a white cuticle layer decreased as the cuticle layers became more disconnected and increased as the layers became more connected (Fig. 6d).Regarding the "Po(cl)" parameter, the PH 3 -S strain had a lower closed porosity than the PH 3 -R strain (p < 0.01).The connectivity parameters indicated that the PH 3 -R strain had a more continuous and tight cuticle structure than the PH 3 -S strain.

Summary of the transcriptome and Micro-CT data with validation using chitin measurement and qRT-PCR
In the PH 3 -R group, the number of individuals per given mass (0.5 g) was higher by 19.2% than in the PH 3 -S group.Total chitin contents in a single insect from the PH 3 -R strain was slightly lower than that from the S strain, but the difference was not statistically significant (Fig. 7a).To selectively investigate and validate gene expression levels, eight genes involved in the chitin biosynthesis pathway were analyzed by qRT-PCR.The eight genes included those encoding glucosamine-6-phosphate N-acetyltransferase (gna), phosphoacetylglucosamine mutase (pagm), UDP-N-acetylglucosamine pyrophosphorylase 1 (uap1), chitin synthase 1 and 2 (chs1 and chs2), N-acetylglucosaminidase 2 and 3 (nag2 and nag3), and chitinase (cht13) (Fig. 7b).Significant differences were observed in the expression of uap1, nag3, and cht13 between PH 3 -S and PH 3 -R strains, whereas no differences were observed in the expression of gna, pagm, chs1, and chs2 (Fig. 7c).The expression level of uap1 was lower in the PH 3 -R strain than in the PH 3 -S strain.Two genes, nag3 and cht13, which are associated with the conversion of chitin to N-acetylglucosamine (GlcNAc), exhibited higher expression levels in the PH 3 -R strain than in the PH 3 -S strain.In particular, the expression of cht13 was constantly upregulated in the PH 3 -R stain, with its expression level being 4-to 8-fold higher than that of the PH 3 -S strain, based on transcriptome data and qRT-PCR validation results (Figs. 3b and 7c).
Overall, trehalose, which is a precursor of chitin synthesis, exhibited a higher concentration in the PH 3 -R strain than in the PH 3 -S strain due to the high expression of the trehalose-specific transporter 1 gene, tret1 (Figs.3b  and 8).Additionally, the expression of nag3, cht5, cht5, and cht13, which are associated with the conversion of chitin to GlcNAc, was significantly upregulated and chitin contents were lower in the PH 3 -R strain, suggesting that the strain may utilize chitin for other purposes, such as energy metabolism (Fig. 7).Consequently, the chitin may experience a reduction in the lamellae source, which consists of chitin polymers and nanofibrils in the PM (Fig. 8).To compensate for the reduction in chitin in the PM, the expression of PM protein-related genes, pmp-1a, b, c, and peritrophin-1, was significantly upregulated in the PH 3 -R strain (Figs.3b and 8).Conversely, hydrocarbons in the cuticle layers were produced in the oenocytes through fatty acid synthesis, and the expression of fatty Fig. 6 Connectivity parameter of PH 3 -susceptible (S) and resistant (R) strain of Tribolium castaneum using Micro-CT analysis.Data were obtained by analyzing the same part, divided into PH 3 -S and R (Union), and male and female (Subset).Schematic illustrations are representing for a mean number of objects per slice (2D Obj.N), b average object area per slice (Av.Obj.Ar, µm 2 ), c structure model index (3D SMI) and d closed porosity (percent) (Po(cl), %).All data are presented as mean ± standard deviation.Union group was analyzed using t-test and Subset group was analyzed using one-way ANOVA, followed by Tukey's multiple comparisons.Asterisk represents the statistic mark in Union group (*p < 0.05; **p < 0.01; ***p < 0.001).SF, susceptible female; RF, resistant female; SM, susceptible male; RM, resistant male acid synthesis-related genes, fabp and far1, was upregulated in the PH 3 -R strain.Furthermore, a more continuous and tight cuticle structure was observed in the PH 3 -R strain.Overall, PH 3 exposure led to notable changes in the cuticle layers and chitin metabolism, thereby enhancing PH 3 resistance in T. castaneum (Fig. 8).

General target sites of fumigants including phosphine
Insects maintain a unique respiratory system in which O 2 and CO 2 gases are exchanged and are referred to as the DGC, which begins with all spiracles closed.Spiracles are the gateways to the tracheal system and are usually tightly closed (C-phase or closed-spiracles).During the C-phase, there is minimal gas exchange and O 2 concentration in the tracheal system decreases, whereas CO 2 concentration increases [20].When O 2 concentration reaches 4%-5%, the spiracles flutter, O 2 flows into the tracheal system and reaches the same concentration as that of mitochondrial respiration.This is called the flutter (F)-phase, in which O 2 introduced by diffusion is removed from the tracheal spaces, and approximately 15%-25% of CO 2 generated by mitochondria is released into the tracheal spaces by diffusion.In this phase, water loss is minimized.CO 2 produced by the mitochondria is partially released to obtain sufficient O 2 [20].The insect maintains O 2 concentration at approximately 4%-5% during the F-phase and the gas concentration gradient in the spiracle is adjusted to 16%-17%.In addition, CO 2 concentration is maintained at 4%, which is fourfold lower than that in the mitochondria, where CO 2 is produced.Therefore, CO 2 accumulates during the F-phase, in turn, causing the spiracles to open (O-phase) and a high amount of O 2 is introduced into the tracheae [20].The unique respiratory system of insects makes them an excellent target site for fumigants they are exposed to.Target pests exposed to fumigants suppress the inflow of fumigants by closing their spiracles [24], which alters the gas concentration gradient in the tracheal system and suppresses the emission of CO 2 generated from the mitochondria, thereby affecting the overall metabolism in the insect pest and causing imbalances.Kim et al. [24] found that the opening rate of spiracles in Plodia interpunctella (Hübner) larvae was significantly reduced upon fumigation with ClO 2 [24].The metabolic imbalance in insect pests has been demonstrated by several studies.Sitophilus oryzae (Linnaeus), which is resistant to PH 3 , regulates energy production by reducing the use of total O 2 and the production of CO 2 and H 2 O [25].The PH 3 -R strain of S. oryzae exhibits reduced oxidative phosphorylation even in the absence of PH 3 treatment.Therefore, nonsynonymous gene mutations associated with the subunit of complex 1 of the mitochondrial ETC are induced, which in turn, suppresses the overall production of H 2 O and CO 2 [25].

Transcriptomic analysis of genes expressed in PH 3 -susceptible and PH 3 -resistant T. castaneum strains
Changes in the expression of genes in the T. castaneum PH 3 -S and PH 3 -R strains are illustrated in Fig. 2e.Notably, considerable changes in gene expression were observed in the chitin synthesis pathways when compared to other metabolic pathways (Fig. 3a).A transcriptomic analysis of PH 3 -R T. castaneum was performed post-PH 3 exposure and compared to untreated specimens, and these results identified 44 significantly upregulated cytochrome P450 (CYP) genes, such as CYP6A1-like, CYP345A1, CYP4C1, and CYP346B1, in both susceptible and resistant groups following exposure [31].In contrast, our study focused on transcriptomic analyses of both PH 3 -S and PH 3 -R T. castaneum strains without PH 3 fumigation.Emphasis was placed on functional annotation analysis of innate  S1 and Additional file 2: Table S3 properties in the PH 3 -S and PH 3 -R strains in this study.Notably, even in the absence of PH 3 exposure, CYP6BQ9, a brain-specific gene, was significantly upregulated in the PH 3 -R strain compared to the S strain (Additional file 2: Table S3).In this regard, there is a difference of significant genes between our study and Oppert et al. [31] due to the difference of PH 3 exposure.
Insect cuticles are composed of structural cuticle-binding proteins and chitin, which form strong exoskeletons to protect insects from various environmental stresses [30].Cuticle layers contribute to insecticide resistance in various insect pests with different types of cuticle modifications.Some insecticide-resistant insect pests have thicker epicuticles to inhibit insecticide penetration to the target sites [7,45].A 50% reduction in the radiolabeled insecticide ( 14 C-deltamethrin) penetration was observed in the deltamethrin-resistant Anopheles gambiae (s.l) strain compared to the deltamethrin-susceptible strain based on an insecticide-impregnated paper assay, and relatively high amounts of cuticular hydrocarbons were detected in relation to the upregulation of cytochrome P450 4G16 (CYP4G16) and CYP4G17 in abdominal oenocytes, which catalyze hydrocarbon synthesis [7].CYP4G16 is associated with decarboxylase activity that transforms [9,10] 3 H-octadecanal to n-heptacosane, whereas P4504G17 does not perform such an activity [7].
The epicuticle consists of various hydrocarbons, polysaccharide-binding proteins, free fatty acids, and wax esters; therefore, certain epicuticle components have been modified based on the type of pesticide an insect pest is exposed to.For example, the PH 3 -R strain of T. castaneum has been reported to have different concentrations of glycerolipids (1.13-53.1-foldvariance) and phospholipids (1.05-20.0-foldvariance) [4].The observation was also made in a stored grain insect pest Rhyzopertha dominica (Fabricius), which is resistant to PH 3 , with higher glycerolipid and phospholipid contents being observed in the resistant strain than in the control strain [4].The authors suggested that the variation in glycerolipid and phospholipid contents could be associated with the long-term energy reservoirs and barriers that protect insect mitochondria from PH 3 toxicity, respectively [4].Furthermore, significant differences have been observed in the content of the cuticular hydrocarbon, 3-methylnonacosane between the PH 3 -R and PH 3 -S strains, while 2-methylheptacosane content significantly differs between the PH 3 -resistant and PH 3 -susceptible strains of R. dominica [3].Alnajim et al. [3] demonstrated that high concentrations of cuticular hydrocarbons could lead to PH 3 -resistant insect strains to reduce PH 3 penetration into their bodies, in turn, protecting them from fumigant toxicity [3].
However, the two studies did not measure the thickness of the epicuticles of PH 3 -R strains of T. castaneum and R. dominica.
Conversely, some insecticide-resistant insect pests, such as Drosophila melanogaster (Meigen), harden their endocuticles by forming more laminated structures [40].DDT penetration in the 91-R strain of D. melanogaster decreased by approximately 1.3-fold when compared to the control strain of Canton-S flies, and the reduced penetration was attributed to the high amounts of cuticular hydrocarbons and relatively thick laminated cuticles [40].In the TEM analysis, endocuticles of the 91-R strain of D. melanogaster had a denser laminated structure than those of the Canton-S populations.Similarly, hardening of endocuticles with an increase in the number of chitin layers (98 layers in resistant strains) has been reported in β-cypermethrinresistant strain of Bactrocera dorsalis (Hendel) when compared to the susceptible strain (75.67 layers) [26].The phenomenon contributed to an increase in cuticle sizes among insecticide-resistant oriental fruit flies, which inhibited β-cypermethrin penetration, when compared to the susceptible flies.In this study, cuticle thickness in the PH 3 -R strain of T. castaneum adults remained unchanged when compared to that in the PH 3 -S strain.Therefore, our results and the findings of Alnajim et al. (3,4) suggest that PH 3 resistance in T. castaneum is associated with reduced PH 3 penetration because of altered cuticle composition [3,4,6].However, demonstrating the hypothesis is a challenge because of high amounts of certain hydrocarbons.Generally, high hydrocarbon concentrations in relation to insecticide resistance are associated with thick cuticles [6,7,45].In another study, transcriptomic findings for PH 3 -resistant Cryptolestes ferrugineus (Stephens) revealed significant enrichment of mitochondrial and cuticular protein genes [12], mirroring our results (Fig. 3).These results supported the evidence of a relationship between PH 3 resistance and cuticle.Nevertheless, previous studies have yet to produce evidence for alterations in the cuticle layer itself with only gene expression evidence.In this context, our combined Micro-CT and transcriptomic data robustly suggest that cuticle layer modifications correlate with PH 3 resistance (Figs. 3 and 6).Chen et al. [12] also proposed a potential molecular mechanism for PH 3 resistance based on RNA-seq data, speculating about potential changes in cuticle thickness or composition [12].However, our findings suggested no cuticle thickness alterations but significant variations in the connectivity of cuticle layers (Fig. 6).This discovery of cuticle layer modifications represents a pioneering insight into the mechanisms underlying PH 3 resistance.

Modifications of cuticles in the PH 3 -resistant strain of T. castaneum
Micro-CT analysis revealed the development of new cuticle forms in the PH 3 -R strain of T. castaneum (Figs. 6 and 8), which is a unique feature when compared to other insecticide-resistant insect pests.The cuticles of the PH 3 -R strain were more anti-osteoporotic than those of the PH 3 -S strain (Fig. 6c) which implies that the PH 3 -S strain of T. castaneum is not osteoporotic.The PH 3 -S strain of T. castaneum formed a cuticle that minimizes water loss and protects them from other environmental stresses.However, the PH 3 -R strain of T. castaneum exhibited increased production of hydrocarbons from abdominal oenocytes and development of a new cuticle matrix with a more connected and tight structure using different building blocks (Fig. 8).Consequently, the cuticle thickness of the PH 3 -R strain of T. castaneum remained unchanged, although more pores were observed in the modified cuticles (Fig. 6d and Additional file 1: Fig. S1).The observation was notable because PH 3 does not easily penetrate the insect body through pores that serve as reservoirs for penetrated PH 3 gas.
According to previous studies, fumigants can pass through the cuticle layers as well as through the spiracles [8,15,28,39,42].Tsao et al. [42] observed that CO 2 expiration in the American cockroach, Periplaneta americana (Linnaeus) increased during fumigation using glucosinolate breakdown products as their mode of fumigant toxicity in the insect respiratory system [42].Price [32] also found that the uptake of 32 PH 3 was considerably lower in the resistant strain of R. dominica and the absorbed PH 3 was rapidly metabolized in the PH 3 -R strain [32].However, a study conducted by Price a year later [33] revealed that PH 3 resistance in R. dominica was acquired from the active exclusion of PH 3 rather than metabolic detoxification [33].A recent study based on transcriptomic analysis revealed the upregulation of P450s in the PH 3 -R strain of R. dominica [47], however, the role of P450s in the metabolism of PH 3 remains unclear.
The observed changes in the cuticle structure of T. castaneum can be explained by the following two reasons.First, low absorption of PH 3 , which may be the primary reason for the changes in cuticle structure in T. castaneum as trace amounts of PH 3 penetrated under anaerobic conditions [28].Second, water loss in the PH 3 -R strain of T. castaneum due to insufficient water production caused by the respiratory system shut down during PH 3 fumigation.The observation has also been made in PH 3 -R strains that are not exposed to PH 3 .

Miscellaneous changes in the PH 3 -resistant strain of T. castaneum
Some of the PM protein-encoding genes in the PH 3 -R strain of T. castaneum were upregulated (Figs.3b and  8).The PM consists of chitin and proteins [2], and it is critical in T. castaneum for the compartmentalization of digestive enzymes and ingested foods, as well as providing protection from abrasive food particles and enteric pathogens [18].PM proteins are resistance factors that are known for their role in promoting Cry1Ac resistance in the cotton bollworm, Helicoverpa armigera (Hübner) [22].However, further studies are required to determine whether the expressional changes in PM proteinencoding genes are associated with PH 3 resistance in T. castaneum or are merely follow-up to changes in chitin synthesis.Gut microbes can degrade the PH 3 absorbed by insects [38].

Conclusion
The present study revealed significant upregulation of chitin-related genes, including chs1, cht5, cht10, and cht13 in the chitin catabolic process term, and PM protein-encoding genes pmp-1a, pmp-1b, and pmp-1c in the chitin-binding term based on functional annotation analysis.A significant differential expression of genes was observed in the insect cuticle structure, with chitin being the primary cuticle component.The cuticle structure of the PH 3 -R strain was examined and compared to that of the PH 3 -S strain of T. castaneum using Micro-CT analysis.More extended and tighter cuticle structure was observed in PH 3 -R T. castaneum than in PH 3 -S strains, with no changes in cuticle thickness.Insect pests that are resistant to contact insecticides are known for their mechanism of rapidly metabolizing and excreting the penetrating insecticides using detoxifying enzymes, such as CYP450s and carboxylesterases, and decreasing affinity to the target sites due to structural changes induced by gene mutations.However, as a fumigant, PH 3 affects the respiratory system of insect pests by penetrating insect bodies through spiracles in the respiratory tract.Therefore, the resistance acquired by insects following exposure to fumigants differs from conventional resistance mechanisms associated with contact insecticides.Micro-CT analysis revealed a unique cuticle structure that was more connected and tight to prevent PH 3 penetration into the insect body during fumigation when the spiracles are completely closed and to potentially reduce water loss due to insufficient water production.In conclusion, T. castaneum resistance to PH 3 toxicity induced significant changes in the insect cuticles and modifications of metabolic processes that balance of O 2 use and CO 2 production.Future studies should focus on establishing strategies to reduce the resistance acquired by insect pests toward PH 3 toxicity, which could facilitate the reduction of PH 3 content in atmospheric environments, including agricultural environments.

Fig. 1
Fig. 1 Cuticle layer structure analysis using Micro-CT.a Scheme of sample preparation.b Summary of Micro-CT scanning and reconstruction process.c Post-processing of 3D rendering structure data for comparing parameters between phosphine-susceptible and resistant strain.Reconstructed cross-sectional images were thresholded with a grayscale index (254-255) (binarization).Cutting processing was only performed for 371 slides from first abdominal sterna to the posterior direction and the parts corresponding to the legs were cut (cutting).Removing noise from the inside (removal of inside noise) and outside (removal of outside noise) of the cuticle layer was followed

Fig. 2
Fig. 2 Transcriptomic analysis of phosphine (PH 3 ) susceptible (S) and resistant (R) strain of Tribolium castaneum.a, b Mortality and comparison of PH 3 -resistant parameter based on lethal concentration (LC) 50 (a) or lethal concentration-time (LCT) 50 .FC, fold-change.c Genotyping of P45S mutation on dihydrolipoamide dehydrogenase.rr, not detected.;Rr, heterogeneity; RR, homogeneity of P45S.d Principal component analysis (PCA) plot of transcriptome of PH 3 -S and PH 3 -R group.e Heatmap of significant genes (p < 0.05) with hierarchical clustering.Color key means Z-score

Fig. 3
Fig. 3 Functional annotation analysis of differentially expressed genes (DEGs, twofold-change, p < 0.05) compared to PH 3 -susceptible (S) with resistant (R) Tribolium castaneum.a Top terms of DAVID functional annotation (*p < 0.05; **p < 0.01; ***p < 0.001) in the cellular component (CC), molecular function (MF), biological process (BP), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.b Volcano plot of DEGs compared to R/S stain.Genes were selected by chitin-related genes based on DAVID analysis.Gene's full name with details is shown in Additional file 2: TableS3

Fig. 7
Fig. 7 Comparisons of total chitin contents and quantitative reverse transcription polymerase chain reaction (qRT-PCR) in phosphine (PH 3 )-susceptible (S) and resistant (R) Tribolium castaneum.a The number of adults per equivalent mass (0.5 g) for Calcofluor staining and the ratio of total chitin mass per single adult T. castaneum.b Scheme of chitin biosynthesis pathway representing differences between PH 3 -S and PH 3 -R strain of T. castaneum.c Relative gene expression level involved in chitin biosynthesis pathway using qRT-PCR.All data for chitin contents and relative mRNA level presented as mean ± standard deviation and analyzed using a t-test.Asterisk represents the statistic mark (*p < 0.05; **p < 0.01)

Fig. 8
Fig. 8 Summary illustration of transcriptome and Micro-CT data of phosphine (PH 3 )-susceptible (S) and resistant (R) strain of Tribolium castaneum.Significant gene expressions were represented using red color (upregulated in R) and blue color (downregulated in R) in transcriptomic analysis.Gene symbols are shown in Additional file 1: TableS1and Additional file 2: TableS3