Henan province in China, which contributes a quarter of China’s annual wheat harvest, is the core wheat production area. Notably, the first suction trap in Henan was constructed in Yuanyang county, Xinxiang city (35.01 N, 113.69 E). The suction trap is 8.8 m tall and collected the tiny flying insects weekly from April to June 2018–2020. Morphological characteristics were used to identify the three main wheat aphids, S. miscanthi, R. padi, and S. graminum, among the trapped samples. All samples were deposited in the Institute of Plant Protection, Henan Academy of Agricultural Sciences, Zhengzhou, China. They were stored in 75% ethanol and maintained under −20°C until DNA extraction.

The population dynamics of migration wheat aphids was summarized, and their migrating peaks were revealed. The HYSPLIT (HYbrid Single-Particle Lagrangian Integrated Trajectory) model is designed to simulate the trajectory of substances transported and dispersed through atmosphere using gridded meteorological data (Stein et al., 2015; Rolph et al., 2017). It is widely used to simulate the global dust distribution, volcanic ash dispersion, and pollutant transport (Stein et al., 2007; Stunder et al., 2007; Wang et al., 2011). Tiny insects like aphids have limited flight capacity. The flight speed of aphids is approximately 0.9 m s⁻¹ (Robert, 1987), and the aphid flights are controlled by the wind when they move above at approximately 1 m from the ground (Parry, 2013). Hence, the HYSPLIT model is also suitable to calculate the migration trajectory of tiny insects with meteorological data.

In this study, the migration trajectories of S. miscanthi and R. padi were simulated with the online HYSPLIT model (https://www.ready.noaa.gov/HYSPLIT_traj.php). In these analyses, the meteorological data were obtained with the one-degree GDAS (Global Data Assimilation System) model, the model calculated star times were set on the 02:00 UTC time (10 a.m. in Beijing time) on the aphid migration peak days, and the flight heights were chosen as 10, 50, and 100 m AGL (above ground level). The aphid migration trajectories were simulated in 24 h.

In this study, the infections of three main aphid bacterial symbionts, Serratia symbiotica, Hamiltonella defensa, Regiella insercticola, and one common insect bacterial symbiont, Wolbachia, were detected in S. miscanthi and R. padi with specific primers listed in Table S1 . In these specific amplifications, aphid total DNA was extracted from a single aphid using an Ezup Column Animal Genomic DNA Purification Kit (Sangon Biotech, Shanghai), following the manufacturer’s recommendations. Before DNA extraction, every aphid was washed with 70% ethanol and sterile water several times to remove surface contamination. Aphid elongation factor-1α gene was used as a reference to evaluate the extracted DNA quality, and the low-quality ones were excluded in the bacterial symbiont detections. The 16S rRNA sequences of the known aphid bacterial symbiont strains were retrieved to uncover the phylogenetic positions of the identified bacterial symbiont strains of S. miscanthi and R. padi. The sequences were aligned using the MUSCLE program in MEGA 7.0 (Kumar et al., 2016). Phylogenetic analysis was performed in IQ-TREE 1.6.12 (Nguyen et al., 2015). The support for each node was assessed by resampling 5,000 ultrafast bootstraps (Hoang et al., 2018). The best substitution model was selected by Bayesian information criterion in ModelFinder (Kalyaanamoorthy et al., 2017). The phylogenetic tree was visualized in Figtree 1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

Herein, the bacterial composition of S. miscanthi and R. padi was further uncovered by amplicon sequencing. We randomly selected 15 individuals of S. miscanthi and R. padi collected in the same year and divided them into five replicates. The genomic DNA of the pooled aphids in each replicate was extracted using the DNeasy Blood and Tissue Kit (QIAGEN, Germany) according to the manufacturer’s instructions. The possible surface contamination was removed as described above. An approximately 460-bp fragment of the bacterial 16S rRNA V3–V4 region was amplified with primers 341F: CCTACGGGNGGCWGCAG and 806R: GGACTACHVGGGTATCTAAT (Qin et al., 2021). The PCR reactions were performed in triplicate in a 50-μl mixture containing 5 μl of 10 × KOD Buffer, 5 μl of 2 mM dNTPs, 3 μl of 25 mM MgSO₄, 1.5 μl of each primer (10 μM), 1 μl of KOD Polymerase, and 100 ng of template DNA. PCR amplifications were carried out with the following program: 94°C for 2 min, followed by 30 cycles at 98°C for 10 s, 62°C for 30 s, and 68°C for 30 s and a final extension at 68°C for 5 min. The amplicons were further purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer’s instructions. Subsequently, the purified amplicons were quantified using the ABI StepOnePlus Real-Time PCR System (Life Technologies, Foster City, USA), pooled in equimolar and generated the sequencing libraries. The next sequencing was performed by the Gene Denovo Biotechnology Co. (Guangzhou, China) on an Illumina HiSeq 2500 platform with a 2×250-bp paired-end method according to the standard protocols.

The reads were trimmed using Btrim with the cutoff for average quality scores higher than 20 over a 5-bp window size and a minimum length of 180 bp (Kong, 2011). The reads ranging from 242 bp to 244 bp were employed in the following analyses with QIIME2 (version 2022.11) (Bolyen et al., 2019). Primers and adaptor sequences were removed by cutadapt plugin (Martin, 2011). ASVs (amplicon sequence variants) were obtained by merging the sequences and removing chimera with the DADA2 plugin (Callahan et al., 2016). The taxonomic assignment of representative sequences was carried out using the feature-classifier plugin against the SILVA database (silva-138-99). Four alpha diversity indexes, ACE (abundance-based coverage estimator), Chao1, Shannon, and Simpson, were calculated in QIIME2; Bray–Curtis distance matrix was calculated in QIIME2, and then PCoA (principal coordinates analysis) and NMDS (non-metric multi-dimensional scaling) were visualized with Bray–Curtis distances and plotted in the ggplot2 package. The relative abundance of the top 10 bacterial taxa among the samples and the cladogram of the differential species were analyzed and visualized by the MicrobiotaProcess package (https://github.com/YuLab-SMU/MicrobiotaProcess), with the feature table and taxonomy table obtained in QIIME2.

In this study, most of the migrant wheat aphids were trapped in 2018 (4,051 samples), and 2020 had the least number of trapped migrant wheat aphids (737 samples). In 2019, 1,271 migrant wheat aphids were trapped ( Figure 1A ). R. padi was predominant among the samples of trapped wheat aphids, whereas S. miscanthi and S. graminum were less common. However, among the three aphid species, the number of trapped aphids was not significantly divergent (ANOVA test, F = 2.334, p = 0.1778). In addition, different patterns of wheat aphid migration were discovered over the years ( Figures 1B–D ). The results indicated that the migrations of R. padi had two distinct peaks during 2018–2020. Specifically, in 2018 and 2019, the first R. padi migration peak occurred in the second week of May, and the second peak occurred in the first week of June. However, in 2020, the two R. padi migration peaks happened 1 week earlier than that of 2018 and 2019. In contrast, only one peak was observed in the migration of S. graminum and S. miscanthi in 2018 and 2019, which occurred in the first and third weeks of May, respectively. The S. miscanthi migration experienced two adjacent peaks in 2020, which occurred in the fourth week of April and the second week of May, respectively. However, no obvious peak was revealed in S. graminum migration in 2020, since at most three samples were trapped in the weeks and a total of 15 samples were collected.

Here, the meteorological data reflecting the local climate conditions near the suction trap (thereafter named Yuangyang), as well as those reflecting the climate conditions of a wider area, such as Xinxiang city, were involved in the statistical analyses to comprehensively uncover the underlying roles of temperature and humidity on the migrations of wheat aphids. The results indicated that there was no significant difference in the temperatures of April and May during 2018–2020, in either the Yuanyang or Xinxiang region (ANOVA test, p > 0.05) ( Figures 2A, B ). However, a significant difference in their relative humidity was observed (ANOVA test, p < 0.05) ( Figures 2C, D ). Furthermore, the multiple comparisons using Tukey’s test revealed that the relative humidity of April 2018 was significantly different with that of 2020 either in the Yuanyang or Xinxiang region. However, the relative humidity of May 2018 was significantly different from that of 2019 and 2020 in the Yuanyang region, and in the Xinxiang region, the relative humidity of May 2018 was significantly different from that of 2019.

To identify the origin of the migrating aphids, we simulated the backward migration trajectories of immigrations. On the other hand, the forward migration trajectories of emigrations were simulated to uncover the destinations of the flying aphids. The results demonstrated similar migration trajectories at the three layers: 10 m, 50 m, and 100 m AGL ( Figures 3 , 4 ). Most S. miscanthi migrations on 7 May 2018 and 28 April 2020 originated from the west and southeast of Yuanyang county, respectively ( Figures 3A, B ). In addition, on 21 May 2019 and 14 May 2020, S. miscanthi migrated to the northeast and northwest of Yuanyang county, respectively ( Figures 3C, D ). The migration trajectories of R. padi are summarized in Figure 4 . The results indicated that most R. padi migrations on 14 May 2018 originated north of Yuanyang country ( Figure 4A ). However, most of the migrating R. padi on 14 May 2019 and 7 May 2020 originated northeast of Yuanyang country ( Figures 4B, C ). According to the simulation of forward migration trajectories, most R. padi would migrate to the south and east of Yuanyang county on 7 June 2018 ( Figure 4D ). On the contrary, most R. padi would migrate to the northeast and north of Yuanyang county on 7 June 2019 and 28 May 2020, respectively ( Figures 4E, F ).

In this study, a total of 237 and 350 samples of S. miscanthi and R. padi were subjected to the specific bacterial symbiont detections, respectively ( Table 1 ). S. symbiotica, H. defensa, and R. insercticola were detected in the aphids; however, no infection of Wolbachia was identified. In S. miscanthi, infection rates of S. symbiotica and H. defensa were 12.23% (29/237) and 1.27% (3/237), respectively. On the other side, in R. padi, infection rates of S. symbiotica and H. defensa were 7.71% (27/350) and 1.42% (5/350), respectively. The infection rates of R. insercticola in S. miscanthi and R. padi were 2.95% (7/237) and 4.29 (15/350), respectively. However, between S. miscanthi and R. padi, no significant difference was observed in the infections of S. symbiotica (χ ² = 2.28, p = 0.13), H. defensa (χ ² = 0.027, p = 0.86), and R. insercticola (χ ² = 0.69, p = 0.40). In the phylogenetic tree ( Figure 5 ), the monophyly of S. symbiotica, H. defensa, and R. insercticola was robustly supported. In these symbiont clades, the strains identified in S. miscanthi and R. padi shared close relationships with that isolated in other aphid species, which further verified their infections.

In this study, a total of 3,874,860 raw reads (average 129,162 reads per sample) were yielded in the 16S rRNA V3–V4 amplicon sequencing. After trimming, denoising, and chimera removal, 8,518 representative ASVs were obtained, and their taxonomic positions were identified against the SILVA database. The relative abundance of the top 10 bacterial taxa among the samples was revealed at four taxonomic levels ( Figure 6 ). The results showed that compositions of bacterial community harbored in the S. miscanthi 2018 population was different from that of the remaining aphid populations. For instance, at the phylum level, in the 2018 S. miscanthi population, Firmicutes was dominant with 79.60% relative abundance and Proteobacteria was the second main taxon with 8.46% relative abundance. However, in the remaining aphid populations, Proteobacteria was dominant with a relative abundance ranging from 34.49% in the 2018 R. padi population to 78.47% in the 2019 R. padi population. A similar divergence was also observed at the order and family level. At the genus level, Buchnera was identified in all samples with an average relative abundance of 11.26%. Furthermore, two aphid facultative symbionts belonging to Serratia and Rickettsiella and three common insect bacterial symbionts, Arsenophonus, Rickettsia, and Wolbachia, were identified. In S. miscanthi, Serratia was the dominant aphid facultative symbiont with an average abundance of 8.41%. On the other side, Arsenophonus was the dominant aphid facultative symbiont in R. padi with an average abundance of 23.30%. In S. miscanthi, the average abundance of Rickettsiella was 10.11% much higher than that of R. padi, which had only 0.13% average abundance. Rickettsia was mainly identified in the 2019 S. miscanthi population with a 9.49% abundance. Wolbachia was only identified in the 2019 R. padi population with a 0.05% abundance. The results of biomarker discovery indicated that Arsenophonus was the most differentially abundant bacterial taxon in R. padi ( Figure S1 ). Additionally, within S. miscanthi and R. padi, the abundance of Rickettsiella, Arsenophonus, and Serratia was varied among the years. Their lowest relative abundance was observed in the 2018 aphid populations ( Figure S2 ). In S. miscanthi, the relative abundance of Rickettsiella in the 2018 and 2019 aphid populations was significantly different. A similar result was also observed in Serratia ( Figure S2A ). In R. padi, the relative abundance of Arsenophonus in the 2018 aphid population was significantly different from that in the 2019 and 2020 aphid populations ( Figure S2B ).

In this study, the alpha diversity among the samples was varied ( Figure S3A ). Furthermore, it was found that R. padi had a relatively higher richness and evenness than S. miscanthi ( Figure S3B ). In beta diversity analyses, both PCoA and NMDS analyses indicated that the bacterial community of the 2018 S. miscanthi population was distinct from the other samples ( Figure S4 ). The samples of S. miscanthi trapped in 2019 and 2020 had similar bacterial communities. In R. padi, the bacterial communities were separated by years in PCoA. However, in the NMDS analysis, bacterial communities of the R. padi samples trapped in 2018 and 2020 and that of S. miscanthi samples trapped in 2019 and 2020 were not well separated.

In this study, the amount of migration wheat aphids varied among years. Previous studies indicated that temperature and humidity are the key environmental factors influencing aphid flights (Parry, 2013). Herein, we analyzed the underlying influences of temperature and humidity on the wheat aphid migrations with two types of meteorological data. The results showed that the difference in relative humidity among the years was significant, indicating that relative humidity probably influenced the wheat aphid flight. We also found that the migrating populations of the three wheat aphids varied from each other. R. padi was the most abundant among the identified migrating aphids, followed by S. miscanthi, and S. graminum was the least abundant. This result was consistent with the findings of our aphid surveying in the wheat field (data not shown). In Yuanyang country, R. padi was the dominant wheat aphid species, particularly in the late wheat development stages characterized by higher temperatures. Previous studies demonstrated that R. padi performed better in high-temperature habitats (Asin and Pons, 2001) and showed greater ecological plasticity when exposed to harsh environments (Zhu et al., 2021). We hypothesize that the R. padi population would proliferate faster with a temperature increase than other wheat aphids; thus, more winged R. padi would be induced due to the low host plant quality and crowding.

A previous study indicated that in Chinese northern winter wheat plant regions, aphids would probably migrate into wheat fields during the wheat heading to flowering stages, and then migrate out in the wheat milk-ripe stage (Li et al., 2014a). In the Yuanyang country, the periods of wheat heading to flowering stages ranged from around late April to early May, and the milk-ripe stage of wheat ranged from around late May to early June. Hence, in this study, we roughly used the half of May as the threshold to distinguish the types of aphid migrations; those that occurred earlier were categorized as immigrations, while those that occurred later were categorized as emigrations. Although alate aphids will undertake appetitive flight in short distances (Loxdale et al., 1993), long-distance migrations are also observed in various aphid species (Kring, 1972). This study found that R. padi and S. miscanthi had different migration patterns. Specifically, two migration peaks were observed in R. padi, identified as immigrations and emigrations. However, in S. miscanthi, the migration patterns varied over the years. Only immigrations were observed in 2018, whereas only emigrations were observed in 2019. In 2020, two adjacent migration peaks were uncovered. These results suggest that the field population of S. miscanthi in Yuanyang county in 2018 likely included the migrated aphids, but few S. miscanthi migrated out in late June. However, in 2019, most S. miscanthi field specimens might have been overwintering individuals that migrated out in late May.

Wind speed and direction are the key factors influencing the initiation, path, speed, distance, and duration of aphid flight (Parry, 2013). HYSPLIT is widely used to simulate the migration trajectories of insects, such as Helicoverpa armigera (Feng et al., 2009), Spodoptera frugiperda (Westbrook et al., 2016), Pantala flavescens (Cao et al., 2018), and Anopheles mosquitos (Huestis et al., 2019). Herein, HYSPLIT was used to simulate the trajectories of R. padi and S. miscanthi on their peak migration days. The results indicated that the immigrations typically originated from the south, whereas the north was the primary destination of the emigrations. These findings were consistent with the previous conclusions that wheat aphid migration trajectories in Chinese wheat fields were typically from south to north in line with the monsoons (Li et al., 2014a). However, there were also some exceptions. In 2018, it was predicted that the S. miscanthi would migrate from the west, whereas the predicted direction of the R. padi in emigrations was south.

Previously, using specific PCR, H. defensa and R. insercticola have been detected in S. miscanthi (Li et al., 2014b). Recently, the bacterial communities of S. miscanthi are investigated with 16S amplicon sequencing (Wang et al., 2022a). However, the examined aphids in these studies are collected in wheat fields; the bacterial community of the migrant wheat aphids is less known. In this study, using specific PCR and phylogenetic analysis, the infections of S. symbiotica, H. defensa, and R. insercticola were identified in S. miscanthi and R. padi migrants. However, no significant difference was observed in their infection rates. Amplicon sequencing indicated that Proteobacteria and Firmicutes were the dominant bacterial taxa in the wheat aphid migrants. Moreover, we identified two aphid facultative symbionts (S. symbiotica and R. viridis), and three common insect bacterial symbionts (Arsenophonus, Rickettsia, and Wolbachia). Arsenophonus has been detected in diverse aphid groups (Jousselin et al., 2012), whereas R. viridis has been identified in pea aphid (Tsuchida et al., 2010). To our knowledge, it is the first report of them in S. miscanthi and R. padi. Furthermore, Arsenophonus was identified to be the most differentially abundant bacterial taxon in R. padi. Additionally, to uncover the potential interactions between the infections of bacterial symbionts and aphid flight, the relative abundance of Rickettsiella, Arsenophonus, and Serratia, was used in the statistical analyses. In both S. miscanthi and R. padi, the lowest relative abundance of the bacterial symbionts was observed in the 2018 aphid populations. Furthermore, the specific PCR results revealed that the 2018 populations of S. miscanthi and R. padi had the lowest rates of bacterial symbiont infection. Since most of the aphid migrants were trapped in 2018, it indicated that aphid bacterial symbionts did not probably promote the wheat aphid migrations. However, this hypothesis needs to be further verified by comparing the aphid flight capabilities of bacterial symbiont infected and uninfected populations.

In conclusion, in this study, we used the data obtained from the Yuanyang suction trap site, uncovered the population dynamics of migrant wheat aphids in the main wheat planting region in China, and simulated their migration trajectories. Furthermore, we revealed the potential roles of climate conditions on the aphid migrations. Additionally, the interactions between wheat aphid migrants and bacteria were investigated with specific PCR and amplicon sequencing.