To study the reassortment pattern in Asia as comprehensively as possible, we assembled all the sequences from the Influenza Virus Resources at the National Centre for Biotechnology Information (NCBI) (https://ftp.ncbi.nih.gov/genomes/INFLUENZA/). The sequences sampled in Asia were then extracted. The sequences we downloaded are coding sequences. For each segment, we only consider the longest coding sequence. We deleted sequences with illegal characters in the codon region or incomplete length. Sequences without clear epidemiological information were also deleted. We did not consider the laboratory sequences. The genome of a virus would be preserved only when eight segment sequences met the above standards. After that, we obtained 7,075 genomes of IAVs. The distribution of the number of viruses in each subtype was shown in Figure 1A. There are a large number of IAVs with H1N1 and H3N2 subtypes, both exceeding 1,000 strains. In the distribution of the host (Figure 1B), viruses that infected humans, waterfowl, domestic birds, and swine account for the vast majority.

Based on the distance measure DMk proposed by Wei et al. (16), we introduced our feature extraction method for each segment sequence. A segment sequence is composed of 64 kinds of codons arranged in order. Since a codon may appear multiple times in the sequence, we can record the position of each occurrence. All the occurrence positions can form an array, recorded as {_{p0, p1}, p₂, ⋯ , p_n}. In order to facilitate the calculation of the interval, the position at the 0th occurrence p₀ = 0 is added to the position sequence. Then we used the reciprocal of the interval between each occurrence and the previous occurrence as its occurrence density, calculated as

Where p_i represents the position of the i th occurrence of the codon, and p_i−1 represents the position of the i-1 th occurrence of the codon. d_i reflects the interval, or density, between the i th occurrence and the i-1 th occurrence of the codon. In order to describe the order of the density array, we accumulated the density array to obtain the partial sum array, calculated as:

The partial sum array is calculated directly based on the density array, for example, s₁ = d₁, s₂ = d₁+d₂, and so on. This array is only determined by the occurrence position of the codon, and the position of the codon can be reversed based on the array. Shannon entropy reflects the order relationship of the elements in the sequence. In order to calculate Shannon entropy, we constructed a discrete probability distribution Q = {q₁, q₂, ⋯ , q_n}, where each probability is calculated as qi=si/∑j=1nsj. Then the Shannon entropy can be calculated as:

For each codon, we can calculate a Shannon entropy, and the Shannon entropy of all codons finally constitutes a feature vector. In the actual process, we removed the feature of 3 stop codons, and finally got a 61-dimensional vector for each segment sequence.

In this section, we proposed a sequence clustering method based on the self-organizing map (SOM). SOM is an unsupervised artificial neural network. Different from the general neural network training based on the reverse transfer of the loss function, it uses a competitive learning strategy and relies on the competition between neurons to gradually optimize the network. The neighbor relationship function is used to maintain the topological structure of the input space. In our clustering method, the algorithm takes the 61-dimensional vectors extracted from the segment sequences as input and maps them to the nodes of the two-dimensional grid. Segment sequences with high similarity in the input space will be mapped to the same node. We used the minisom library to implement our algorithm. First, we calculated the number of nodes in the output layer according to the empirical formula of clustering n=5N, where N is the number of genomes. The size of the output layer is size=n. The calculations above were rounded up. Then we initialize the SOM algorithm according to the following parameters: the output layer size was set to 61; the sigma parameter was set to 2; the learning rate was set to 1, and the neighborhood function was set to triangle; weight initialization method was set to PCA. The feature vector of each sequence was trained iteratively. After training, each vector was input to the SOM network and the winning node was the clustering type of the virus segment. Then we calculated the contour coefficient of the clustering result and recorded the corresponding parameters. Finally, we changed the parameters and repeated clustering. The clustering result with the highest contour coefficient was regarded as our final result for the segment clustering.

Based on our clustering for each segment, we assigned a number x to each type and recorded it as segment_x. The virus sequences contained in each type were considered to have the same segmental ancestors. Therefore, if the eight segment types of a virus can be combined from segment types of other viruses, this virus may be produced by the reassortment of these viruses. There may be a situation where different combinations of multiple viruses can produce the eight segment types of the virus. For example, the eight segment types of virus A may be combined from segment types of virus B and virus C, or combined from segment types of virus D and virus E. Then we need to consider which combination is more likely to produce A by reassortment. To this end, we considered epidemiological information of the viruses such as host, location, and time. For viruses with multiple reassortment possibilities, we thought viruses with the similar host, location, and time are more likely to reassort. Our non-phylogeny-dependent reassortment detection method was shown in Figure 2.

Based on the feature vectors extracted from each sequence, we represented the IAV genome as a 488-dimensional feature vector composed of eight segment feature vectors. In this way, the similarity of the IAV genome can be characterized by the similarity of feature vectors. We then constructed a similarity matrix of all IAVs and used it to construct a similarity network. The minimum spanning tree algorithm was used to simplify our network, and the cys layout in cystoscope was used for display, which was shown in Figure 2. In the virus similarity network, the subtypes of different virus genomes have been clearly distinguished (Figure 3A), which proves that the feature vector we extracted is reasonable. It should be noted that swine is a mixed vessel, which connected the genome of the avian influenza virus and the genome of the human infected influenza virus (Figure 3B). We also found that some human infected virus at the top of the figure. The genomes of these viruses were similar to those of domestic birds. Therefore, domestic birds also play an important role in the process of the avian influenza virus infecting humans.

Using our feature extraction and SOM-based clustering method, we divided eight segments of influenza A virus into different types, respectively, which was shown in Table 1. According to the combination of hemagglutinin (HA) and neuraminidase (NA) proteins on the surface, the influenza A virus can be divided into different subtypes. At present, HA can be divided into 18 subtypes (H1-H18), NA can be divided into 11 subtypes (N1-N11). Therefore, segments with the same subtype should be grouped when we cluster HA and NA segments. Based on this, we checked the clustering results of HA and NA segments, which were shown in Figure 4A. Each type contained only one subtype and the same subtypes were clustered together, indicating that our clustering method has no problem in distinguishing between subtypes. We also analyzed the number of types corresponding to different HA and NA subtypes, which was shown in Figure 4B. For HA subtypes, the number of types in H1, H3, H5, and H9 is relatively large, indicating that the lineages of these subtypes are complex. For NA subtypes, there are relatively more types in N1 and N2. The clustering results for all segments can be found in Supplementary Tables 1–8.

We performed reassortment detection on all the viruses in the dataset one by one in the order of their sampling time and identified 516 reassortment events. The details of each reassortment event can be found in Supplementary Table 9. According to the virus subtypes involved in reassortment, we divided reassortment events into intra-subtype reassortment and inter-subtype reassortment. Intra-subtype reassortment means that the subtype involved in reassortment is single, while inter-subtype reassortment means that the subtypes involved in reassortment are diverse. We analyzed the number of two reassortment patterns and found that the intra-subtype reassortment pattern was more common, which was shown in Figure 5A. In particular, we analyzed the situation of two reassortment patterns of eight common subtypes including H1N1, H3N2, H5N1, H5N6, H6N2, H6N6, H7N9, and H9N2, which was shown in Figure 5B. Viruses with H1N1, H3N2, H9N2 subtypes had a higher proportion of intra-subtype reassortment, while viruses with H5N1, H5N6, H6N2, H6N6, and H7N9 subtypes were more likely to be produced by reassortment of viruses with different subtypes.

Although the intra-subtype reassortment pattern was frequent, we were more concerned about inter-subtype reassortment, since such reassortment pattern is more likely to produce pandemic strains. To this end, we made statistics on the frequency of other subtypes in the inter-subtype reassortment events of each subtype. We found that there were strong reassortment signals between H1N1 and H3N2, H6N2 and H6N6, H9N2, and other subtypes, which was shown in Figure 6. H1N1 and H3N2 are common subtypes of human influenza viruses. The viruses with these two subtypes circulate seasonally every year, leading to frequent reassortment. The 2009 H1N1 pandemic is a typical example. The virus that caused this pandemic was a triple-reassorted strain whose genome is derived from H1N1 avian virus, H1N1 classical swine virus, and H3N2 human seasonal virus. We also noticed the reassortment relationship between H6N2 and H6N6. There were few subtypes related to H6N6 reassortment and more subtypes related to H6N2 reassortment. However, the reassortment association between the two subtypes was still much more than other subtypes. In 2011, the H6N6 avian influenza virus was isolated from swine in Guangdong Province, China (17), suggesting that the H6 subtype avian influenza virus has the potential to infect mammals and even humans. Therefore, the reassortment between H6N2 and H6N6 needs to be further monitored. It should be noted that we found H9N2 played an important role in the reassortment of other subtypes since its proportion was the first or second. Based on this, we analyzed the specific reassortment relationship between H9N2 and other subtypes.

H9N2 avian influenza virus spreads worldwide and is widespread in poultry in Asia (18–20). A recent study has shown that H9N2 had replaced H5N6 and H7N9 as the dominant avian influenza virus subtype in chickens and ducks in China (21). At the same time, there is also evidence that H9N2 contributed to the emergence and evolution of human infection with avian influenza viruses such as H7N9, H10N8, and H5N6 (22–24). In our study, we found that H9N2 provided internal segments for viruses with other subtypes. Among 61 inter-subtype reassortments related to H9N2, 27 events were H9N2 providing internal segments for H7N9. We also found four cases of H9N2 providing internal segments for H10N8 and four cases of H9N2 providing internal segments for H5N1. H9N2 occasionally provided internal segments for H5N2, H6N2, H6N6, H3N8, and H1N2. In addition, we found that these subtypes also provided internal segments for H9N2. These H9N2 viruses that obtained internal segments from other subtypes returned to the H9N2 gene pool and provided segments for subsequent reassortment. In general, the reassortment pattern related to H9N2 is complicated. It can provide internal segments for other subtypes, and can also obtain internal segments from other subtypes.