Evidence for a link between the Atlantic Multidecadal Oscillation and annual asthma mortality rates in the US

An association between climatic conditions and asthma mortality has been widely assumed. However, it is unclear whether climatic variations have a fingerprint on asthma dynamics over long time intervals. The aim of this study is to detect a possible correlation between climatic indices, namely the Atlantic Multidecadal Oscillation and Pacific Decadal Oscillation, and asthma mortality rates over the period from 1950 to 2015 in the contiguous US. To this aim, an analysis of non-stationary and non-linear signals was performed on time series of US annual asthma mortality rates, AMO and PDO indices to search for characteristic periodicities. Results revealed that asthma death rates evaluated for four different age groups (5–14 yr; 15–24 yr; 25–34 yr; 35–44 yr) share the same pattern of fluctuation throughout the 1950–2015 time interval, but different trends, i.e. a positive (negative) trend for the two youngest (oldest) categories. Annual asthma death rates turned out to be correlated with the dynamics of the AMO, and also modulated by the PDO, sharing the same averaged ∼44 year-periodicity. The results of the current study suggest that, since climate patterns have proved to influence asthma mortality rates, they could be advisable in future studies aimed at elucidating the complex relationships between climate and asthma mortality.

The atlantic multidecadal oscillation and pacific decadal oscillation. The AMO and the PDO are internal modes of climate variability. The AMO is a multi-decadal low-frequency alternation of high and low sea surface temperatures (SSTs) in the North Atlantic 21,22 with a reported oscillation period of about 65-70 years 21 . Enfield et al. 18 documented a 65-80 year cycle from 1856 to 1999 in North Atlantic SST data. In 2011, Knudsen et al. 23 explored the past 8000 years through proxy records and found a quasi-persistent ∼55to 70-year AMO signal linked to the ocean-atmosphere internal variability.
The AMO is known to influence rainfall, temperature and pressure in many regions of the Northern Hemisphere 24 , frequency of Atlantic hurricanes 25 , North American climate 26 and river flows 18 as well as the hydro-climatic conditions in other areas, like rainfall variability in Northeast Brazil 27 and occurrence of droughts in the Sahel 28,29 .
The PDO is a pattern of ocean-atmosphere climate variability in the Pacific which appears as warm or cool surface waters in the Pacific ocean poleward of 20°N 30 . Over the past century, the amplitude of this climate pattern has varied irregularly at interannual (time periods of a few years) to interdecadal (time periods of multiple decades) time scales. This climate pattern affects the Northern Hemisphere climate, the North Pacific ecosystem, North American precipitation, stream flow, surface temperature anomalies 31,32 , fluctuations of the Asian monsoon 33 , and a modulation of ENSO 34 . Minobe [35][36][37] found that PDO fluctuations were most energetic at periodicities in the 15-25 and 50-70 years bands; Chao et al. 38 found evidence for oscillatory variations at 15-20 and near to 70 years. Willmott et al. 39 jointly analysed the November-April PDO index and surface temperature and precipitation gridded data and found that the warm phases of the PDO tend to coincide with anomalously warm temperatures in northwestern North America, northern South America and northwestern Australia, and anomalously cold temperatures in Eastern China, Korea, Japan, Kamchatka, and the Southeast US and Mexico. Gedalof and Smith 40 identified 11 regime shifts in the PDO record since 1650, with the most recent occurring in 1976/77. With a 23 years average duration of a regime, they suggested that another shift is expected by around the end of this century. Many of the climate anomalies associated with the PDO are broadly similar to those associated with ENSO phases, though generally not as extreme 30,37,41 .
One recent study 42 analysed the observed temperature time-series since the beginning of the 20 th century in nine climate regions of the US and found that their time variations and oscillations were accurately reproduced by a combination of AMO and PDO oscillations with a monotonic signal associated with anthropogenic CO 2 warming. The small temperature decrease observed in the period 1938-1974 and the large temperature increase in the period 1980-2000 were thought to be caused, respectively, by the superposition of a downward trend of the oscillatory mode on an upward trend of the monotonic mode and by the superposition of an upward trend of the AMO oscillatory mode on the upward trend of the monotonic mode. Yao et al. 43 found periods of warming slowdown from 1880 to 2012, caused by the superposition of the AMO and PDO oscillations on the steady monotonic warming signal, presumably caused by increasing amounts of atmospheric CO 2 .
Subsequent studies have been focused on the idea that droughts occurred over the contiguous US (in 1996US (in and in 1999US (in -2002 were associated with North Atlantic warming (positive phase of the AMO) and northeastern and tropical Pacific cooling (negative phase of the PDO). The multidecadal oscillations in the behavior of the North Atlantic Ocean thus plays a crucial role for long-term explanation and predictability of drought frequency in the US 44 .

Results
Asthma mortality and climatic data. North America annual asthma death rates from 1950 to 2015 are shown in the first four plots of Fig. 1, divided by age groups (from the oldest to the youngest, from top to bottom). Mean death rate value over the entire time interval increases from the youngest to the oldest age group  Correlation matrix between asthma death rates, AMO, PDO, and TA. Table 1 shows the correlation coefficients (and their statistical significance) between asthma death rates for the four age groups, the AMO and PDO indices, and TA. To properly calculate correlations 45,46 , all raw data were first standardized, by subtracting the mean and dividing by the standard deviation, and detrended, by subtracting the linear trend from the original signal. The AMO index is positively and significantly (p < 0.05) correlated with asthma death rates for only the 25-34 and 35-44 age groups (25-34 yr, p = 2.06 × 10 −8 r = 0.62; 35-44 yr, p = 1.17 × 10 −10 r = 0.69), as is TA (25-34 yr, p = 0.035 r = 0.26; 35-44 yr, p = 0.006 r = 0.33). For both AMO and TA, the correlations with the other age groups turn out not be statistically-significant. The AMO index is positively and significantly correlated with TA. The PDO index is positively correlated with the AMO, TA, and the asthma death rates for the two youngest age groups, while negatively correlated with the two oldest age groups; however, none of the PDO correlations turn out to be statistically significant 47 . The four age groups are all significantly (p < 0.05) correlated with each other.
Signal analysis results. The CEEMD analysis revealed five intrinsic mode functions (IMFs) plus the trends (IMFs 6) (see Supplementary Information) for all analyzed signals shown in Fig. 1. A first visual analysis of all IMFs shows that from IMF2 to IMF4 the asthma death rate signals (all age classes), AMO, and PDO indices, recorded significative periodicities. Instead, all IMFs1 recorded predominantly noise, and all IMF5 of asthma death rate (all age classes) and TA contain an incomplete cycle. For these reasons we exclude IMFs 1 and 5 from successive analysis.
To obtain analytical results of periodicities recorded in the analyzed signals, we apply REDFIT and weighted wavelet Z-transform (WWZ) on a total of 23 IMFs (see Supplementary Information). The spectra results with periodicities above 95% Confidence Interval (CI) are reported in Table 2.
These results show us that only the IMF4 component of all signal (Fig. 2), except for the TA, have a main peaks in the ~37/~57 years periods range, all above the 95% CI and continuous throughout the investigated time window (Fig. 3).
Finally, applying a bad-pass filter to the AMO (filter 45-65 years), PDO (filter 35-55 years), and the asthma mortality raw data (filter 30-50 years), all insignificant frequencies were removed in order to obtain a more easy inspection of the relation between the asthma death rates and climatic indices. In Fig. 4, the raw and band-pass filtered data of all age groups death rates, AMO, and PDO are depicted. The visual inspection of all significant data suggests synchronicity between AMO and asthma death rates for all age groups, highlighting a slight desynchronisation effect starting from ∼1980. With regard to PDO vs AMO and asthma death rates, negative synchronicity, centered at ∼1960 and ∼2000 respectively, occurs. During the 1980-85 time interval, the high positive PDO period seems to show null relationship with other data.

Discussion
This is the first study comparing historical US annual asthma death rates for four age groups, from childhood (5-14 yr) to adulthood (35-44 yr), with relevant climatic indices, which are known to modulate drought periods in this geographic area. Droughts conditions are among the main contributors to environmental factors leading to or exacerbating asthma.
The studied period covers a long time-span (66 years), from 1950 to 2015, during which an increasing trend in surface air temperature has been observed and the historical records of precipitation, streamflow and drought indices have pointed toward increased aridity since 1950 over many land areas 48 , including the regions considered in this paper. www.nature.com/scientificreports www.nature.com/scientificreports/ time-span of asthma mortality rate data in US has been encompassed for the first time, our findings cannot be adequately compared to previous results. Nonetheless, results obtained for shorter time intervals (e.g., 1968-1987 52 ) are in agreement with our findings on increasing/decreasing sequences and trends.
The observed variations in asthma death rates cannot be an effect of the four revisions of the International Classification of Disease codes (ICD) occurred in the last 56 years (see Methods), since they do not coincide with sharp inversion or trend modification which has occurred. Instead, a combination of different factors including increasing prevalence of asthma, changing patterns of both disease severity and medical care, as well as an enhanced recognition of the disease might be the major drivers of the observed sequence of fluctuations 52 . In particular, when considering the introduction of bronchodilators from 1960's and the inhaled steroid treatment from 1990's, no evident effects on US annual asthma death rates could be recognized. Noteworthy, an increase trend of annual asthma mortality can be observed starting from 1980's and encompassing the 1990's up to 1996 in all the age groups. According to previous studies, this finding might be attributed to different factors relative to prescription patterns 53 , medication misuse, underuse, overuse, as well as toxic effects 54 .
On the other side, environmental variables may play a role that could mainly be observed when long time intervals are taken into account.
Our results pointing toward a correlation between increase/decrease sequences of asthma death rates and the AMO index oscillation may indicate a phenomenon previously unrecognized in the link between climate variability and respiratory events. Indeed, a comparable evolution of the AMO index and the fluctuations of annual asthma mortality rates for the four age groups has been observed in the 30-year period 1950-1980. Periods of maximum (minimum) in asthma death rates correspond to periods of positive (negative) values in the AMO index, thus suggesting the AMO, through its links to drought events, as a possible risk factor for asthma mortality in the contiguous US. Additionally, it is known that the 1945-1957 drought period, a European-to-North America continent event recorded by several different datasets, indices, and proxies (e.g.: Global Historical Climatic Network observation-based dataset, Palmer Drought Severity Index, tree rings 55-57 ; Reworked Coccoliths 58 and references therein), was for the first time correlated to US asthma death rates. After 1980, a combined superimposing effect of increasing TA and positive PDO occurred. This effect impairs the visual correlation between mortality data and AMO index for the last 35 years. In order to bypass this problem and to highlight the common cyclicity, the EMD analysis allowed us to find main peaks centred at ~57 and ~46 years (mean ∼51 yr) for AMO and PDO, respectively. This result is in agreement with previous literature data 23,42 . In the mortality IMF4 signal, main peaks centred at ~37 years were found from the younger to the older age groups. The comparison of all IMF4 data in respectively long-term frequency band shows that the asthma death rates signals are in phase with the AMO oscillations during the 1950-1980 time interval, where maxima in mortality correspond to positive values in AMO signal. After 1980, a slight desynchronization of about 5-7 years has been recognized, probably related to the co-occurrence of a long positive phase of the PDO index, the 1980/90 abrupt climatic shift, and a period characterized by the two largest volcanic eruption of the century, El Chichón and Pinatubo, occurred in 1983 and in 1991, respectively.   62 . Three regime shifts (1970s, 1980s and 1990s), distinguished by marked increases in temperatures or by abrupt temporal changes across different biophysical systems, have been identified in the last few decades [63][64][65][66][67][68] . Whilst documented at ocean basin or regional scales, the mechanisms behind these events, their environmental interactions, and the synchrony and scale of their effects around the globe are poorly understood. Thus, there is a considerable research gap with many disparate observations by different scientific disciplines, without a comprehensive overall assessment. In addition to its modulation of the carbon cycle 69 ; diseases (Vibrio cholerae 70 ); biotic, physical and chemical land components 71,72 ; freshwater 73 ; precipitation 74 ; marine 75,76 ; cryospheric 77 and atmospheric 59,60 systems, herein we hypothesize that the 1980 shift might have modulated also the US asthma mortality cyclicity. www.nature.com/scientificreports www.nature.com/scientificreports/ This study certainly has some limitations, such as not having considered gender classifications or having limited the analysis to a very large area without focusing on small regions, potentially useful to highlight the contribution of local climate or other local effects.  www.nature.com/scientificreports www.nature.com/scientificreports/ supported by the finding of a common mean periodicity of about 44 years among the analysed variables, which introduces some new elements in the research on asthma epidemiology, though cyclical fluctuations have been already documented in the epidemiology of other diseases (e.g., pneumoniae 78 ; influenza 79 ).
Furthermore, since drought variations in US have been shown to be controlled by AMO and PDO indices, they may be hypothesized as an emerging risk factor for asthma mortality. This is in agreement with the recent statement of the Global Asthma Network: "Environmental factors are much more likely than genetic factors to have caused the large increase in the number of people in the world with asthma, but we still do not know all the factors and how they interact with each other and with genes" 80 .
At last, the results of the current study suggest that patterns of climate variability may be emerging risk factor for asthma mortality rates. Thus, it would be advisable to develop future correlative studies including the aforementioned climatic indices in order to elucidate the complex relationships between climatic factors and asthma mortality. Finally, further studies are envisaged to test the applicability of our methodology also in different geographic areas and confirm or extend our findings.

Methods
Asthma mortality data. The annual asthma mortality dataset was downloaded from the US Centers for Disease Control and Prevention -National Center for Health Statistics -National Vital Statistics System website (1950 to 1998 data downloaded from https://www.cdc.gov/nchs/nvss/mortality/hist290.htm -1999 to 2015 data downloaded from https://www.cdc.gov/nchs/data_access/Vitalstatsonline.htm#Mortality_Multiple). This dataset is a compilation of mortality data by 10 yr-age groups, race, gender, and cause of death (according to the International Statistical Classification of Disease and Related Health Problems -ICD code), as reported annually.
Asthma deaths were identified using the following ICD codes: ICD-7 241 (from 1950 to 1967), ICD-8/9 493 (from 1968 to 1998), and ICD-10 J45/46 (from 1999 to 2015). Population age classes were categorized into children (5-14 years), youth (15-24 years), young adults (25-34 years), and middle-aged adults (35-44 years), following the World Health Organization classification. Older age groups were not taken into account in view of the possible lack of precision of death certificates for asthma in such ages. To study long-term trends, the crude asthma mortality rates (n/100000) for the 4 age groups were calculated 49 . climatic indices. The detrended monthly mean AMO and of PDO data from 1900 to 2015 were downloaded from the National Oceanic and Atmospheric Administration (https://www.esrl.noaa.gov/psd/data/timeseries/ AMO/ -https://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/PDO/). The AMO and PDO annual mean values were calculated by applying a 12 points running average. Annual mean data of the North America land temperature anomalies (TA, 1950(TA, -2015 were downloaded from https://www.ncdc.noaa.gov/cag/global/time-series/ northAmerica/land/ytd /12/1950-2015. Signal analysis. In order to single out characteristic periodicities in the time-series which we analysed, the analysis of the non-stationary (frequency changes with time) and non-linear signals was performed by applying the Empirical Mode Decomposition algorithm (EMD) by Huang et al. 81 .
The Ensemble Empirical Mode Decomposition (EEMD) and its complete variant (CEEMD) are adaptive, noise-assisted data analysis methods that improve on the ordinary Empirical Mode Decomposition (EMD) by Huang et al. 81 . This decomposition provides a powerful method to look into the different processes behind a given time series data, and provides a way to separate short time-scale events from a general trend.
Empirical mode decomposition is a form of adaptive time series decomposition method where the basis functions are derived from the signal itself, while in some standard forms of spectral analysis methods like Fourier and wavelet analysis, the basis functions are fixed as sine and cosine for the first and as mother wavelet functions for the second. The decomposition process produce IMFs that are singular function representing an oscillatory mode with one instantaneous frequency that needs to satisfy two criteria: • In the whole time series, the number of extrema and the number of zero crossings must be either equal or differ at most by one; • At any point in the time series, the mean value of the envelopes which is defined by local maxima (upper envelope) and local minima (lower envelope) is equal to zero.
This decomposition technique rests on the assumption that any complicated signal can be decomposed into a finite and often small number of components called "Intrinsic Mode Functions" (IMF) 81 , each of them representing an embedded characteristic simple oscillation on a separated time-scale. The data were detrended prior to analysis. We have added in the Supplementary Materials a complete discussion of EMD analysis and its upgrades, as well as of REDFIT and Foster's wavelet spectral technique.
With the purpose of comparing the dominant periodicities recorded in the asthma death rates with the same order periodicities documented in the reference global signal (AMO and PDO), we applied a band-pass filter using the wavelet multi-level decomposition and reconstruction technique, which is invertible and thus suitable for filtering data. In particular, we used the multiresolution analysis (MRA) algorithm to decompose a signal into scales with different time and frequency resolution organized according to a hierarchical scheme 86 .