Introduction

Influenza A viruses (IAVs) are infectious agents that cause many epidemics and seasonal influenza pandemics (Chen et al. 2018; Zhou et al. 2015). The natural reservoirs of IAVs are aquatic wild birds and mammalians such as pigs, horses, swine (Koutsakos et al. 2019; Slusher et al. 2014; Urbaniak et al. 2014). Apart from H18N11 and H17N10, almost all subtypes of IAVs are grouped in avian influenza viruses (AIVs). The H1N1, H2N2 and H3N2 subtypes of IAVs are human-adapted subtypes due to their ability to α-2,6 sialylated glycan receptors in human epithelial cells in the upper respiratory tract (Cauldwell et al. 2014; Ge and Wang 2011; Imai et al. 2012; Jayaraman et al. 2011; Long et al. 2019; Mair et al. 2014; Makarova et al. 1999; Schäfer et al. 1993; Smith et al. 2009; Viswanathan et al. 2010). Since 1968, the H2N2 has stopped circulating in human but, like other H2 subtypes, continues to circulate in birds and sporadically in pigs. The H5N1, H7N7, H7N2, H7N3 and H9N2 subtypes, circulating in domestic poultry, are not adapted to human but have been reported to cause infections without sustained interhuman transmissions (Alonso and Diambra 2020). The H5, H6, H7, H9 and H10 AIV subtypes have crossed the species barrier to infect humans via infected poultry (Kessler et al. 2021) via acquiring the ability to infect human due mutations in hemagglutinin leading to increased affinity for α-2,6 sialylated glycan receptors in human (Mair et al. 2014). (Kessler et al. 2021)The IAV subtypes that have caused pandemics in humans are H1N1, H2N2, H3N2, H5N1, H7N9, H7N7 and H9N2, in descending order of the number of pandemic deaths in human (Koutsakos et al. 2019; Long et al. 2019; Philippon et al. 2020; Song and Qin 2020; Sutton 2018; Tanner et al. 2015; Wu et al. 2020; Yamaji et al. 2020). The IAV subtypes are potential causal agents of pandemic in the human population in the future. However, (Gao et al. 2013; Koopmans et al. 2004; Marinova-Petkova et al. 2017; Meineke et al. 2019; Skowronski et al. 2006; Song and Qin 2020; Terebuh et al. 2018; Tong et al. 2018; Tweed et al. 2004; Um et al. 2021; Westenius et al. 2018) current management of human influenza A is limited to vaccines and antiviral treatments targeting viral proteins (Uyeki et al. 2019). Due to their rapid mutation capacity, IAVs have a strong evasion mechanism to the host immune system (Chang et al. 2020; Huang et al. 2021; Li et al. 2012; Peacock et al. 2017; Peng et al. 2014; Sitaras et al. 2014), thus reducing the effectiveness of antivirals targeting viral proteins. To counteract these problems, the development of drugs targeting the virus-host interactions that promote IAV infection and replication is necessary. Thus, analyzing the implications of virulence factors on host processes such as gene expression is crucial.

Among the promising targets of the treatment and prevention of influenza, the non-structural protein 1 (NS1) is a multifaceted protein majorly known for its role in evasion of the host immune system and RNA binding and is considered to be the main viral weapon for combating the innate immune response and promoting viral multiplication in the hostile environment of the host (Han et al. 2019; Zhang et al. 2019). NS1 is a key regulator in the replication and pathogenesis of avian influenza (Rosário-Ferreira et al. 2020). Studies have shown that NS1 interacts with a variety of host proteins to participate in viral replication, gene transcription regulation and induction of cell pathological changes, etc. (Nemeroff et al. 1998; Nogales et al. 2018). In host cells, NS1 of H3N2 can interact with nucleoli to regulate the stress related to cell cycle progression based on the rRNA gene promoter (Yan et al. 2017). Moreover, the NS1 protein of IAV subtype H5N1 interacts with the heat shock protein Hsp90 and is involved in the apoptosis induced by IAV via caspase cascade (Zhang et al. 2011). It was also reported that mammalian adaptation of avian influenza A virus implicates gradual shift in NS1 and that interplay among immune evasion and virus evolution is fundamental in IAV mammalian adaptation (Chauché et al. 2018).

The uneven usage of the synonymous codons encoding the same amino acids is also called codon usage bias (CUB), which is widespread in biology (Deng et al. 2020). CUB is mainly determined by codon in the third position of nucleic acid and CUB analysis has been widely used in virus genome studies and contributes to the understanding of the shifts of the evolution and genetics in the virus (Cutter et al. 2006). Liu et al. (2011) found that the synonymous codon usage pattern is a mixture of coincidence and antagonism to that of the host cells. Viruses tend to use different codon usage patterns to improve overall virus survival and evade the host immune system (Moratorio et al. 2013). Previous studies indicated that similarities of relative synonymous codon usage (RSCU) of genes between pathogens and their host can impact gene expression in the host and that the pathogens depend on their hosts for their multiplication (Alonso and Diambra 2020; Jitobaom et al. 2020; Kula et al. 2018; Michely et al. 2013; Miller et al. 2017; Nambou and Anakpa 2020; Pepin et al. 2008). Recently, codon usage similarities between viral genes and human genes has been demonstrated and indicated that the infection (viral proteins engaging the host cell proteins) may induce pathological effects as a result of the expression and translation of RNA transcripts and the effects of proteins from these transcripts (Maldonado et al. 2021). Translational adaptation of human viruses to the tissues they infect has been also demonstrated and indicated that viruses adapt their codon usage toward their tropism and hijack the host defense machinery (Hernandez-Alias et al. 2021). Based on this evolution of molecular compatibility among the pathogens and their hosts, we can retrieve key information for infectious diseases, especially the translation effects more likely affected by the evolution of CUB. Thus, the co-evolution analysis based on CUB similarities of viral genes and host genes represents a way to identify host genes that may have transcripts dysregulated in response to infection. Previously, based on RSCU similarity analysis, we identified key genes involved in COVID-19 and predicted potential therapeutics for this disease (Nambou and Anakpa 2020). Proceeding so will open new avenues for understanding the pathogenesis of human infection by the IAV subtypes; however, there is no report of such analysis regarding IAV. Thus, analyzing the similarity between RSCU values of virulence genes of IAV and human gene RSCU values is potential in providing novel significant information on the pathogenesis of IAV. Because of its multifarious implications in the infectious pathogenesis, revealing the CUB of the NS1 protein will help understand the evolution and adaptability of IAVs to human as well as cellular and molecular changes induced by these viruses in human (Luo et al. 2020). However, few studies have reported the codon usage pattern of the NS1 protein of the IAV subtypes and the genome-wide human genes with CUB similar to that of this viral molecule has yet to be exhaustively scrutinized.

In this study, we downloaded NS1 coding sequences of the IAV subtypes H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2 from NCBI. The codon usage was studied to explore the adaptability and CUB of NS1 coding sequences between these IAV subtypes and human, providing a basis for a better understanding of the evolution and human adaptability of IAV subtypes. Further, RSCU-based correlation analysis allowed the identification of the human genes with CUB similar to that of NS1 gene and their functional roles.

Materials and methods

Sequence data

The complete human genome was obtained from GENCODE (https://www.gencodegenes.org/human/). The nucleotides sequences of IAV NS1 protein sequences were downloaded from the NCBI database after searching with specific keywords. After manual curation, only NS1 protein sequences were included. The IAV types included were those having been reported to cause human infections with number of NS1 sequences available higher than 1000. Thus, H1N1 (45,217 NS1 sequences), H3N2 (52,887 NS1 sequences), H3N8 (4553 NS1 sequences), H5N1 (7486 NS1 sequences), H5N2 (2894 NS1 sequences), H5N8 (1958 NS1 sequences), H7N9 (3957 NS1 sequences) and H9N2 (9375 NS1 sequences) were finally included in this study.

Relative synonymous codon usage (RSCU) analysis

According to the equation by Sharp and Li (1986), we calculated the relative synonymous codon usage (RSCU) values among different codons of IAV nonstructural protein 1 and human gene sequences using the R packages “vhcub” (Anwar et al. 2019). RSCU refers to the probability for a codon to be used between synonymous codons of corresponding amino acids. When RSCU < 1, it means that the codon is used less; conversely, when RSCU > 1, it means that the codon is used more. When RSCU = 1, there is no bias in the use of the codon, with RSCU value of > 1.6 and < 0.6 representing "overrepresented" and "underrepresented" codons, respectively. For calculation of RSCU, please refer to the work by He et al. (2019).

Principal component analysis (PCA)

PCA is a common multivariate statistical method, which can transform the dimension of multiple indicators into several comprehensive indicators. In our study, PCA was performed and visualized using the R packages “ade4”, “ggplot2” and “ggrepel” for analyzing the codon usage of NS1 sequences in different IAV subtypes based on the RSCU values.

Effective number of codons (ENC) analysis

ENC reflects the degree of preference for the unbalanced use of synonymous codons in the codon family. ENC was used to assess the CUB of the IAV NS1 coding sequences using the R package “vhcub” (Anwar et al. 2019). The codon with high expression usually has high codon bias and relatively small ENC. The gene codon with low expression has low codon bias and relatively high ENC. Therefore, ENC can be used to measure the expression level of endogenous genes. The ENC measures CUB for a gene, ranging from 20 (high CUB; a gene uses only one codon) to 61 (no bias; synonymous codons are used equally).

ENC-plot analysis

ENC-plot is based on ENC as the abscissa and GC3 (GC content of the third codon) as the ordinate of a scatter diagram, and can be used to evaluate the role of mutation pressure in codon bias. If mutation pressure only acts on the choice of codon bias, the points in the graph will appear on the top of the expected ENC curve. On the contrary, there may be other factors, such as natural selection pressure, involved. ENC-plot was performed using the R package “vhcub” (Anwar et al. 2019).

Parity rule 2 (PR2) analysis

The PR2 plot was used to explore the effects of mutation pressure and natural selection pressure on the usage of a single gene codon. The figure was drawn with A3/(A3 + U3) as the vertical coordinate and G3/(G3 + C3) as the horizontal coordinate. The central point of the figure is where A = U (0.5) and C = G (0.5), representing the equal role of mutation pressure and natural selection pressure. PR2 plot was performed using the R package “vhcub” (Anwar et al. 2019).

Codon adaptation index (CAI) analysis

The codon adaptation index (CAI) is a parameter, which is used to measure the relative adaptation of viruses and their potential hosts. CAI values range from 0 to 1. A sequence with high CAI is believed to have a strong adaptability to its hosts. CAI analysis was performed by R package “vhcub” (Anwar et al. 2019).

Neutrality analysis

Neutrality plot (GC12 values against GC3 values) was used to evaluate the effects of mutation pressure and natural selection on codon usage. The slope of the regression curve between GC12 and GC3 content represents the mutation pressure (Jia et al. 2015; Sueoka 1988). When the slope = 1.0, there is no natural selection pressure or the natural selection pressure is weak (Jia et al. 2015; Sueoka 1988). On the contrary, when the regression curve deviates from the diagonal, natural selection has a significant impact on the usage of codons (Jia et al. 2015; Sueoka 1988). The GC12 was calculated as the average of GC1 and GC2 from each sequence. The visualization of the neutrality plots was performed using the GraphPad Prism software (GraphPad Software, San Diego, California USA, www.graphpad.com).

Relative codon deoptimization index (RCDI) analysis

The relative codon deoptimization index (RCDI) value is a parameter of codon deoptimization to measure the similarity of codons used by a given gene over the reference genome. The higher genetic similarity between the virus and the host is, the higher the translation ratio is, and the closer the RCDI value is to 1 (Mueller et al. 2006). The lower RCDI value indicates higher adaptability between the virus and the host. However, high RCDI may indicate that some genes are expressed during the latency phase or show low replication levels. RCDI analysis was performed with the “vhcub” (Anwar et al. 2019) library in R.

Similarity index

SiD was used to measure the effect of host CUB on IAV NS1. The value of SiD is between 0 and 1. Higher SiD values indicates that the host has a decisive influence on codon usage. SiD value was calculated using the R package “vhcub” (Anwar et al. 2019).

Synonymous codon usage order (SCUO) analysis

The SCUO (0 ⩽ SCUO ⩽ 1) is the synonymous CUB for the whole sequence. A higher SCUO indicates a higher CUB in the sequence. SCUO was analyzed using the R package “vhcub” (Anwar et al. 2019).

Identification of human genes with CUB similar to NS1 CUB

To disclose human genes probably affected by the transcription of NS1, we calculated the Pearson correlation between the human genes and the NS1 sequences from the IAV subtypes based on the RSCU values. The correlation analysis was performed using the Hmisc library in R. The human genes which were correlated with NS1 sequences with |r|> 0.5 and p < 0.05 were considered as those potentially affected by the NS1 transcription.

Protein–protein interaction network of human genes with CUB similar to that of NS1 gene of IAV subtypes

To construct the protein–protein interaction (PPI) network of the human interactome of NS1 gene of the IAV subtypes, we selected the human genes correlated with IAV NS1 from each subtype as the input list for analysis in Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/). The basic settings of analysis chosen were: full network for the network type; confidence for the meaning of network edges; the active interaction sources included textmining, experiments, databases, co-expression, neighborhood, gene fusion and co-occurrence; and the minimum required interaction score at high confidence was 0.7. The obtained PPI networks were exported for visualization in Cytoscape 3.4.0 software (Cytoscape Consortium, SanDiego, CA, USA) and the MCODE plugin was used for the extraction of PPI clusters with hub genes based on the following settings: degree cut-off ≥ 2, node score cut-off ≥ 0.2, K-core ≥ 2, and max depth = 100.

Functional enrichment analysis

The functional analysis of the genes was performed using the TCGAanalyze_EAcomplete function in the TCGAbiolinks library in R. The enrichment terms with adjusted p value of < 0.05 were retained as significantly enriched terms. The visualization was performed with the same library using the function TCGAvisualize_EAbarplot.

Statistical analysis

The GraphPad Prism version 9.0.0 for Windows (GraphPad Software, San Diego, California USA, www.graphpad.com) was used for statistical processing and the one-way ANOVA followed by Tukeys multiple comparison posttest was used to compare the differences in ENC, RCDI, and CAI values between IAV subtypes. The adjusted p-values < 0.05 indicated statistically significant differences among subtypes. The p-values among groups were represented by an inset network on the figures.

Results

A- and U-ending codons are preferred in NS1 coding sequences of IAV subtypes

As shown in Fig. 1, the results of RSCU analysis were used to estimate the codon usage patterns of NS1 coding sequences from different subtypes of IAV. The NS1 coding sequences tended to use codons with A-/U-ending. For H1N1, 27 preferred codons among which 19 had A-/U-ending and9 were highly used with A-U-ending were identified. For H3N2, 28 preferred codons were identified among which 22 had A-/U-ending and 12 were highly used with 10 A-/U-ending. For H3N8, 27 codons were preferred, among which 19 had A-/U-ending and 10 were highly used with 9 of them having A-U-ending codons. For H5N1, 27 preferred codons, with 16 of them having A-/U-ending were identified; 13 of them were highly used and 9 had A-/U-ending. For H5N2, 28 preferred codons among which 19 had A-/U-endingand 10 were highly used with 9 having A-U-ending were identified. For H5N8, 29 preferred codons among which 19 had A-/U-ending and 10 were highly used were identified. For H7N9, 26 preferred codons among which 21 had A-/U-ending and 12 (10 codons with A-/U-ending) were highly used were identified. For H7N9, 27 codons were preferred, and among them, 20 codons had A-/U-endingwhile 11 were highly used codons.

Fig. 1
figure 1

The relative synonymous codon usage (RSCU) values of the NS1 coding sequences of IAV subtypes

Full size image

ENC of the NS1 coding sequences of IAV subtypes

ENC was used to describe the CUB of NS1 coding sequences of IAV subtypes. A smaller ENC value indicates a stronger CUB. The boxplot showing the ENC values of NS1 coding sequences of IAV subtypes was depicted in Fig. 2. Statistical analysis indicated significant differences among the mean ENC values of pairwise-compared IAV subtypes (Additional File S1). The mean ENC value of NS1 sequences of H1N1 (56.07 ± 2.11) was the highest, followed by those of H3N2 (52.01 ± 1.54), H3N8 (51.56 ± 1.27), H5N2 (51.34 ± 1.97), H5N8 (51.02 ± 1.79), H5N1 (49.98 ± 2.35), H9N2 (49.61 ± 2.11) and H7N9 (48.81 ± 2.23). These results indicated that the CUB of NS1 from H7N9 was the highest, followed by H7N9, H5N1, H5N8, H5N2, H3N8, H3N2 and H1N1, in this order.

Fig. 2
figure 2

The effective number of codons (ENC) values of NS1 of IAV subtypes. The egdes of the inset network indicates the distribution of adjusted p-values among the compared node groups of IAV subtypes

Full size image

Trends in codon usage variations

PCA, a common multivariate statistical method, is used to estimate the relationship between samples and variables. In this study, PCA was used to evaluate the usage of synonymous codons in NS1 coding sequences of IAV subtypes. The first two principal components of the NS1 coding sequences accounted for 30.08% and 28.26%, respectively (Fig. 3A). These values reflected that the first two principal components are the main influencing factors for the codon usage of NS1 coding sequences of IAV subtypes. In addition, PCA was clustered into two categories based on NS1 coding sequences of IAV subtypes, among which some individual points were overlapped, suggesting that these points may have the same codon usage trend. Moreover, Fig. 3B demonstrated that the points of NS1 coding sequences of IAV subtypes were covered by the point of gene sequences from Homo sapiens, indicating that NS1 coding sequences of these viruses may share similar codon usage patterns with their human host. The individual overlap of each virus with H. sapiens was as depicted in Supplementary Figure S1.

Fig. 3
figure 3

Principal component analysis (PCA) for NS1 of IAV subtypes based on the RSCU values. A PCA based on the RSCU values of NS1 of IAV subtypes; B PCA based on the RSCU values of NS1 of IAV subtypes and RSCU values of all genes of H. sapiens

Full size image

ENC-plot analysis

ENC-plots (ENC vs GC3) were generated to explore factors that affect the codon usage of the NS1 coding sequences of IAV subtypes. If codon bias is only affected by natural selection pressure, all points fall below the standard curve. However, points above the standard curve indicate that codon usage is affected by mutation pressure. The points representing the ENC values of NS1 coding sequences of IAV subtypes H7N9, H5N8, and H5N1 were exclusively below the expected ENC curve, suggesting that natural selection pressure is the main factor affecting the codon usage of NS1 in these viruses (Fig. 4). For H5N2, H3N8, H3N2, H9N2 and H1N1, the points were on both sides of the expected ENC curve (Fig. 4), demonstrating that the NS1 coding sequences of H5N2, H3N8, H3N2, H9N2 and H1N1 may be affected by both natural selection and mutation pressure.

Fig. 4
figure 4

ENC-plot analysis of NS1 of IAV subtypes. ENC-plot analysis of the NS1, with ENC, plotted against GC3s of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9, and H9N2 was generated. The red curve represents the expected ENC curve when the CUB is determined by the GC3s composition only. The black dots represent the distribution of ENC values of NS1of AIV subtypes

Full size image

Neutrality plot

In order to determine the levels of natural selection pressure and mutation pressure in codon usage of NS1 coding sequences of IAV subtypes, we performed neutrality plot analysis. As shown in Fig. 5A, Pearson correlation analysis indicated significant negative correlations between the GC12 and the GC3 values of NS1 coding sequences in the IAV subtypes H1N1 (r = -0.28, p = 0.00), H5N2 (r = -0.24, p = 3.22 \(\times {10}^{-20}\)) and H5N8 (r = -0.08, p = 0.008). On the contrary significant positive correlations were recorded for H3N2 (r = 0.30, p = 0.00), H3N8 (0.12, p = 3.18\(\times {10}^{-8}\)), H5N1 (r = 0.05, p = 0.004), H7N9 (r = 0.51, 4.48\(\times {10}^{-133}\)) and H9N2 (0.21, 1.60\(\times {10}^{-0.49}\)). Moreover, as shown in Fig. 5B, the slopes of simple linear regression analysis were -0.11, 0.15, 0.09, 0.03, -0.19, -0.04, 0.27 and 1.19 for H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2, respectively. These results indicated that the levels of natural selection in the codon usage of NS1 coding sequences in H1N1, H3N2, H3N8, H5N1, H5N2, H5N8 and H7N9 were respectively 89%, 85%, 91%, 97%, 81%, 96% and 73%. For H9N2, the slope was higher than 1, indicating exclusive neutrality. These results indicated that natural selection was the main factor affecting the codon usage of NS1 coding sequences in H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2 while mutation pressure was a minor influencing factor.

Fig. 5
figure 5

Neutrality plot analysis of GC3s against GC12s for all the coding sequences based on NS1 of IAV subtypes. A Pearson correlation analysis between the GC3s and GC12s of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9, and H9N2. B Neutrality plots of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9, and H9N2

Full size image

Parity analysis

PR2 bias plot was used to analyze whether CUB exists in NS1 coding sequences of IAV subtypes. When there is no deviation from natural selection and mutation, the points fall in the middle of the graph, and the coordinates of this middle point is (0.5, 0.5). Figure 6 showed that most of the points fall in the upper right quadrant of the PR2 plots from each of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2 viruses, indicating that A3s were more than U3s while G3s were more than C3s, which demonstrated that the codon usage pattern of NS1 in the IAV subtypes was affected by mutation and natural selection pressure.

Fig. 6
figure 6

The PR2 plots of NS1 of IAV subtypes. The AT [A3%/(A3% + T3%)] and GC [G3%/ (G3% + C3%)] bias of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9, and H9N2 were plotted. The center of the plot (both the coordinates is 0.5) indicates a position where there is no bias

Full size image

Codon usage adaptation of NS1 of the subtypes of IAV to human

In order to explore the adaptability and codon usage optimization of NS1 of the IAV subtypes and their human host, we calculated CAI, RCDI and SiD to analyze the adaptability of codons. The CAI values of NS1 in the IAV subtypes were 0.57 ± 0.015, 0.59 ± 0.011, 0.58 ± 0.015, 0.6 ± 0.024, 0.58 ± 0.012, 0.59 ± 0.011, 0.59 ± 0.012, 0.59 ± 0.014 for H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2, respectively (Fig. 7A). Significant differences in CAI were found among the pairwise comparisons of the eight viruses (Additional File S1). These values indicated a high degree of adaptation between the IAV subtypes and human. Furthermore, we analyzed the role of CUB in gene expression. The average value of RCDI was the highest for H7N9 (3.25 ± 0.17), followed by H9N2 (3.22 ± 0.19), H5N8 (3.15 ± 0.13), H3N2 (3.1 ± 0.14), H1N1 (3.06 ± 0.17), H5N1 (3.01 ± 0.19), H5N2 (2.95 ± 0.18) and H3N8(2.92 ± 0.21), indicating that the codon usage deoptimization occurred in this decreasing order (Fig. 7B). Significant differences in RCDI were found among the pairwise comparisons of the eight viruses (Additional File S1). Significant differences in SCUO values were also found among the eight IAV subtypes (Fig. 7C, Additional File S1). In addition, SiD analysis was used to assess the influence of the human on codon usage patterns. We found that SiD values were the highest in H1N1 (0.4919), suggesting that the human exerts more effects on this virus (Fig. 7D). However, the RCDI and SiD values were not significantly different among the studied IAV subtypes. In addition, a negative correlation was recorded between the CAI values of NS1 of the eight studied IAV subtypes and the values of ENC and GC3 (Fig. 7E), which indicated that the highly expressed genes had a superior scale of CUB and preferred the codons with A or U at the synonymous position. Positive correlations were also recorded between CAI, RCDI and SCUO values (Fig. 7E).

Fig. 7
figure 7

Adaptation analysis of NS1 of IAV subtypes. A The codon adaptation index (CAI) analysis of IAV NS1. B The RCDI values of IAV NS1. C The SCUO analysis of IAV subtypes. D The similarity index (SiD) analysis of IAV NS1. E Heatmap showing the Pearson correlation between the ENC, CAI, RCDI and SCUO values of NS1 coding sequences of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9, and H9N2. The egdes of the inset network indicates the distribution of adjusted p-values among the compared node groups of IAV subtypes in A–C. Statistical analysis was not applicable to D

Full size image

Identification of human genes with CUB similar to NS1 CUB and their functional roles

In order to identify human genes with CUB similar to CUB of NS1 gene of the IAV viral subtypes, Pearson correlation based on RSCU values was calculated among human genes and NS1 sequences. We obtained 2147 human genes correlated to NS1 with |r|> 0.5 and p < 0.05 (Additional File S2), which were considered as human genes with CUB similar to CUB of NS1 gene; RSCU values of all of these genes were positively correlated with RSCU values of NS1 gene (Additional File S2). Functional enrichment (Fig. 8) results showed that the human genes with CUB similar to CUB of NS1 were enriched in the biological processes including regulation of transcription (77 genes), DNA-dependent regulation of transcription (26 genes), regulation of RNA metabolic process (27 genes), cellular response to stress (66 genes), defense response (56 genes), response to virus (17 genes), activation of immune response (13 genes), oxidation reduction (76 genes) cell cycle phase (37 genes) and M phase of mitotic cell cycle ( 35 genes). For cellular component, extracellular region (171 genes) was the most enriched term whereas in the category of molecular function, adenyl nucleotide binding (165 genes) was the most represented term. In pathway enrichment analysis, we found that pathways such as “activation of IRF by cytosolic pattern recognition receptors” (22 genes), “role of RIG1-like receptors in antiviral innate immunity” (17 genes), “role of Wnt/GSK-3”_signaling in the pathogenesis of influenza (22 genes). The pathways related to influenza such as “role of PI3K/AKT signaling in the pathogenesis of influenza” (18 genes) and “role hypercytokenemia/hyperchemokinemia in the pathogenesis of influenza” (12 genes) were also enriched. Immune-related pathways such as “complement system”, “PPAR/RXR activation”, “role of cytokines in mediating communication between immune cells” (14 genes), “CD27 signaling in lymphocytes”, and “interferon signaling” (10 genes) were also significantly enriched.

Fig. 8
figure 8

Functional enrichment analysis of human genes with CUB similar to that of NS1 gene of IAV subtypes. The GO and pathway functional enrichment analysis was performed and the top 10 terms were visualized for each GO categories and pathway category

Full size image

Protein–protein interaction network prediction and identification of hub clusters of the human genes with CUB similar to that of NS1 of IAV subtypes and their functional roles

In order to explore whether proteins corresponding to human genes with CUB similar to NS1 gene of IAV subtypes interact with each other, we selected the top-2000 genes from these genes for PPI network construction (Supplementary Figure S2). The network was composed of 1877 nodes with 12,707 edges, and the average number of neighbors was 13.578 (Fig. 9A). To uncover hub genes from this network, we performed network clustering analysis using the MCODE plugin and identified 46 clusters. The hub clusters with the highest score (cluster1, score = 31.273) was composed of 34 nodes and 516 edges (Fig. 9A). Further analysis indicated that the hub genes in cluster1 were enriched in the biological process of nuclear division (12 genes), mitosis (12 genes), M phase of mitotic cell cycle (12 genes), M phase (11 genes), regulation of cell cycle (11 genes) and organelle fission (9 genes) (Fig. 9B). The most enriched cellular components of the hub genes (cluster 1) were spindle (10 genes), microtubule cytoskeleton (11 genes) and condensed chromosome (6 genes) while the most significant molecular functions were ATP binding (9 genes) and adenyl ribonucleotide binding (9 genes) (Fig. 9B). The most significant pathways were “role of CHK proteins in cell cycle checkpoint control” (5 genes), ATM signaling (5 genes) and mitotic roles of polo-like kinase (5 genes) (Fig. 9B).

Fig. 9
figure 9

Hub gene subnetwork extracted from the protein–protein interaction (PPI) network of the human genes with CUB similar to that of NS1 of IAV subtypes. A The PPI network based on the top-2000 (based on Pearson correlation) human genes with CUB similar to that of NS1 of IAV subtypes was constructed and the MCODE plugin was used for extracting the hub gene subnetworks and the hub gene cluster with the highest score was visualized. B Functional enrichment analysis of hub genes extracted from the human genes with CUB similar to that of NS1 of IAV subtypes. The GO and pathway functional enrichment analysis was performed and the top 10 terms were visualized for each GO categories and pathway category

Full size image

Similarities and differences of the human protein–protein interaction networks of NS1 coding sequences from different IAV subtypes

In order to discover the dissimilarities and the commonalities between the IAV subtypes H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2, we retrieved the human genes interacting with each of the IAV subtypes. The intersection analysis indicated important intersection genes between the list of genes from the eight IAV subtypes (Fig. 10A). Genes specific to each of the IAV subtypes were also identified (Fig. 10A). The hub subnetworks from the PPI networks constructed from human genes with CUB similar to CUB of NS1 gene of each IAV subtype were identified and depicted in Fig. 10B. These results indicated that though similarities were found in the interaction networks of IAV subtypes, some specific genes were associated to each IAV subtype, which may explain the differences in the response of human to each type of viral infections. The biological processes and pathways enriched from the list of genes from each subtype was summarized in Additional File S3 and indicated differences and similarities among the IAV subtypes.

Fig. 10
figure 10

Similarities and differences of the PPI networks of human genes with CUB similar to that of NS1 coding sequences from different IAV subtypes. A Intersection analysis of the different PPI networks of human genes with CUB similar to that of NS1 coding sequences from each IAV subtype. B Networks of hub genes of the different PPI networks of human genes with CUB similar to that of NS1 coding sequences from each IAV subtype

Full size image

Discussion

The influenza A virus NS1 protein is a key viral protein that plays a key role in immune escape and the weakening of immune system in the host (Fernandez-Sesma 2007). Herein, starting from the CUB analysis, we aimed to study the factors influencing the expression of NS1 gene of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2 during adaptation to the human host and reveal the human genes whose expression would be regulated by the expression of NS1 in the human host. We first computed the CUB indices and subsequently analyzed the correlation between the RSCU values of human genes and those of different NS1 coding sequences of IAV subtypes H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2. The NS1 gene of IAV subtypes had high CAI values and codons ending with A/U were preferentially used. In addition, natural selection was found as the most determinant factor for CUB. The adaptation indices based on the NS1 coding sequences indicated a high degree of adaptability of different IAV subtypes to the host. Most importantly, we identified the human genes with CUB similar to that of NS1 gene and found that they were majorly involved in biological processes related to transcription and cellular response to stress; these genes were also involved in multiple pathways including those related to influenza, response to virus, defense response and immune response. This is the first study to reveal human genes with CUB similar to CUB of NS1 gene of IAV subtypes. The present findings may help develop vaccines and therapeutics for influenza A.

The CUB analysis suggested that the usage of codons of NS1 coding sequences of IAV subtypes is mainly influenced by natural selection. The RSCU analysis showed that the A-/U-ending codons were those preferred in the coding sequences of NS1 of IAV subtypes, suggesting that NS1 gene of IAV subtypes is less affected by mutation pressure. The results were supported by the findings of Luo et al. (2020). All IAV viruses isolated from humans tend to eliminate GC and CpG components during long-term adaptation to avoid activating the anti-viral immune response (Greenbaum and Ghedin 2015; Luo et al. 2020). The ENC-plot, neutrality plot and PR2 analyses clearly indicated that codon usage of NS1 coding sequences of IAV subtypes is influenced by natural selection in a broader extent and mutation pressure in a minor degree. We also found that the NS1 coding sequences of H1N1 and H3N2 were greatly affected by mutation pressure compared to other IAV subtypes, which may increase their pathogenicity compared to other subtypes and this was supported by the findings that mutations in the NS1 region of a virus will affect the stability and function of the protein (Solórzano et al. 2005). The NS1 protein plays an important role in the process of innate immunity against the host, allowing the IAV to replicate more effectively in the host cell (Fernandez-Sesma et al. 2006; Monteagudo et al. 2019). The virus relies on the host cellular components and metabolism to complete its own replication and assembly. The codon usage pattern of viruses reflects the evolution of viruses to improve their survival and fitness in the host (He et al. 2017; Hu et al. 2011; Luo et al. 2019; Taylor et al. 2017). Therefore, comparing the CUB of IAV helps improve our understandings on the evolutionary process of the virus. We found that the ENC value of the NS1 coding sequences of IAV subtypes was higher than 35, indicating lower CUB. Lower CUB in RNA viruses can avoid the competition on synthetic resources between the virus and the host cells so as to improve the efficiency of viral replication (Butt et al. 2016). NS1 is a nonstructural protein of IAV that plays a pivotal role in the early life cycle of the virus. Moreover, we found that the mean ENC values between the eight subtypes was significant, indicating that the codon usage pattern of these subtypes differed significantly from each other. Numerous studies have reported that the relationship between the virus and the host affects the occurrence, dynamics, and evolution of infectious diseases (Gibbs et al. 2010; Irwin et al. 2012; Puigbò et al. 2010; Rodelo-Urrego et al. 2013; Torres-Pérez et al. 2011). CAI values showed that the NS1 coding sequences of IAV subtypes had strong adaptability to human. The high adaptation of NS1 gene of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2 to human indicated the shaping of their codon usage profiles for improved adaptation to those of human, thus achieving their increased replication rate, which indicates high pathogenicity of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2 in humans. The results of CAI was supported by RCDI analysis because low RCDI may indicate high adaptability of both the virus and the host (Perez-Padilla et al. 2009). However, we found that there are differences in the RCDI values among the eight IAV subtypes, among which the RCDI values of H1N1 and H3N2 were the highest, suggesting that both viruses are better adapted to human. The low RCDI values indicated high similarity and adaptation of codon usage among NS1 gene of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2 and human, allowing these IAV subtypes to express NS1 gene efficiently with an elevated replication rate. The SiD value of NS1 coding sequences of IAV subtypes was about 0.49, indicating that the human is less affected by the usage of NS1 codons. The results were verified by neutrality plot and ENC-plot analysis. Therefore, in this study, we found that IAV NS1 showed a certain adaptability to human, which may be beneficial to IAV infection of human populations. The SiD value of NS1 coding sequences of H1N1 was the highest, suggesting the similarity of the genome between NS1 of H1N1 and human, which explains why H1N1 is reported to be a common human influenza A in the IAV subtypes. The high SiD value of H1N1 NS1 sequences may reflect the length of time the virus has circulated in human hosts (about 1918–2009) compared to the other viruses (H3N2 from the 1960s to present). The major factors shaping CUB in various species are natural selection and mutational pressure associated with genetic drift. In the present work, natural selection was noted as the main factor shaping the codon usage of NS1 coding sequences in H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2. This might reflect the length of time the viruses have evolved with the human. To the best of our knowledge, a systematic analysis of the evolutionary adaptation of NS1 gene of IAV subtypes to human has not been performed so far. Thus, our study is the first to explore this area and our results may be essential for the development of vaccines against these viruses (Baker et al. 2015). Previous studies indicated that codon usage adaptation of NS1 gene in human was involved in the spread of the pandemic Zika virus (de Melo Freire et al. 2015; Freire et al. 2018).

Up to date, the human interactome of virulence genes of IAV is not elucidated. Numerous studies have experimentally revealed proteins that can be affected by the IAV subtypes. For example, it was demonstrated that IAV NS1 hinders the activation of NF-κB, alpha/beta interferon, Jun N-terminal kinase and AP-1 transcription factors (Ludwig et al. 2002; Wang et al. 2000). However, the effect of NS1 in gene regulation in infected human remains to be explored. To date, the use of high throughput approaches in the identification of the effect of IAV subtypes would improve our knowledge on the molecular mechanism of these viruses or their virulence factors. However, studies based on those approaches are limited. Here, using the RSCU correlative approach, we uncovered human genes that could be significantly correlated with the NS1 gene of IAV subtypes H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2. This is the first study to uncover human genes with CUB similar to CUB of NS1 gene of IAV subtypes on a large scale. It is worth noting that the similarity in CUB may be the need for the host cell to rapidly produce the proteins (translation) when the transcript is present. The identification of these human genes with CUB similar to CUB of NS1 gene will be useful for diagnostic, prognostic and therapeutic purposes in human cases of infection by IAV subtypes. Moreover, PPI network construction indicated that human proteins in the network strongly interacted with each other, indicating that protein encoded by human genes with CUB similar to that of NS1 gene of H1N1, H3N2, H3N8, H5N1, H5N2, H5N8, H7N9 and H9N2 are involved in coordinated functions leading to pathological changes in the body. Thus, we analyzed the functional role of these human genes and found that human genes in the interactome were involved in the regulation of transcription. Previous studies have indicated that transcription of genes in the host is affected by pathogens (Kainulainen et al. 2014; Terasaki et al. 2016). Our results corroborated with these findings. Recent single cell RNA sequencing analysis indicated that the complexity of the expression profile of viral genes can lead to multifarious host transcriptional responses during the influenza A virus infection (Sun et al. 2020). Other studies have indicated that influenza virus infection can affect cellular transcription by altering the distribution of the polymerase II on cellular genes, which leads to host shut-off (Engelhardt and Fodor 2006). Thus, the modulation of genes involved in transcription or the development of therapeutics based on these genes will be vital for treating human infections caused by the IAV subtypes. Processes such as “cellular response to stress”, “defense response”, “response to virus”, “activation of immune response”, “oxidation reduction” and “cell cycle phase” were also predicted to be functionally altered by the NS1 of IAV subtypes. These observations confirmed the multifunctional role of NS1 gene especially its involvement in the immune and weakening of host defense following viral infection. Some of these processes have been reported to be affected by IAVs in past studies (Jang et al. 2015; Josset et al. 2008). Moreover, we found that the pathways enriched largely included those associated with influenza (eg. “role of Wnt/GSK-3″_signaling in the pathogenesis of influenza”, “role of PI3K/AKT signaling in the pathogenesis of influenza” and “role hypercytokenemia/hyperchemokinemia in the pathogenesis of influenza”), further confirming the credibility of our approach that the similarity of RSCU of parasite and host may affect the diverse processes in the host and drive the subsequent pathological conditions. The pathological conditions induced by the NS1 gene is probably driven by abnormal regulation of cell cycle as demonstrated by the enrichment of the hub genes of the human PPI network of NS1 of IAV subtypes in processes such as nuclear division, mitosis, M phase of mitotic cell cycle and the pathway of “role of CHK proteins in cell cycle checkpoint control”. Thus, our study not only revealed pathways already known to be affected by NS1 gene of IAV subtypes, but also disclosed new pathways and processes that could be modulated in future, thus enhancing our understanding of the molecular mechanism of IAV infections. The obtained data will also help develop appropriate diagnostic approaches and therapeutics for the infected patients. In our previous work, we identified the human genes whose RSCU values are correlated with those of the genes of respiratory coronaviruses (Nambou and Anakpa 2020). Overlapping of these coronavirus-related human genes (Nambou and Anakpa 2020) with the human genes correlated with NS1 gene of IAV subtypes, we identified 1723 common genes; the functional enrichment of these common genes indicated their involvement mainly in defense response (56 genes) and immune-related processes and pathways (Supplementary Figure S3). These observations indicated that the NS1 gene of IAV subtypes may induce some pathogenetic mechanisms similar to those of other respiratory viruses, namely the coronaviruses.

Our study presents some limitations. We only included IAV subtypes with the number of NS1 coding sequences higher than 1000 in order to keep statistical power; doing so, the other IAV subtypes were excluded from the study, which may lead to some information loss. In addition, we did not include the IAV subtypes that had not infected human, including these viruses may be informative for predicting their effect in case of human infection by these viruses. Furthermore, comparative analysis of the co-adaptation of NS1 gene of IAV subtypes and human to that of other natural reservoirs and intermediate hosts was not performed. Performing such analysis would help understand the differential adaptation of IAV subtypes to their hosts. Moreover, due to the huge number of analyses performed, we did not compare the results from allele A and allele B of NS1 gene, which will be addressed in a future study. In the future, we will also update our research based on experimental and public data and verify our results through in vivo and in vitro experiments.

In sum, our study showed that the CUB of the NS1 coding sequences of IAV subtypes is highly conserved and subtype-specific though similarities were found. The codon usage pattern is mainly affected by the pressure of natural selection and mutation pressure in a lesser extent. Additionally, IAV NS1 showed good adaptability to the host. The findings are vital for understanding the CUB of IAV NS1 and the factors that affect the evolutionary process of IAV, and may contribute to the development of the IAV vaccine. We also identified the human genes with CUB similar to CUB of NS1 gene of IAV subtypes and their function, which may shed light in the pathogenesis of human infection by these viruses and may help develop adequate treatments for IAV-infected patients. The functional role of the predicted human genes needs to be verified by other methods, especially at the translational regulation level.