Summary: The Italian government has been one of the most responsive to COVID-2019 emergency, through the adoption of quick and increasingly stringent measures to contain the outbreak. Despite this, Italy has suffered a huge human and social cost, especially in Lombardy. The aim of this paper is dual: i) first, to investigate the reasons of the case fatality rate (CFR) differences across Italian 20 regions and 107 provinces, using a multivariate OLS regression approach; and ii) second, to build a “taxonomy” of provinces with similar mortality risk of COVID-19, by using the Ward's hierarchical agglomerative clustering method. I considered health system metrics, environmental pollution, climatic conditions, demographic variables, and three ad hoc indexes that represent the health system saturation. The results showed that overall health care efficiency, physician density, and average temperature helped to reduce the CFR. By the contrary, population aged 70 and above, car and firm density, air pollutants concentrations (NO,, O,, PM,, and PM,), relative average humidity, COVID-19 prevalence, and all three indexes of health system saturation were positively associated with the CFR. Population density, social vertical integration, and altitude were not statistically significant. In particular, the risk of dying increases with age, as 90 years old and above had a three-fold greater risk than the 80–to–89 years old and four-fold greater risk than 70–to–79 years old. Moreover, the cluster analysis showed that the highest mortality risk was concentrated in the north of the country, while the lowest risk was associated with southern provinces. Finally, since prevalence and health system saturation indexes played the most important role in explaining the CFR variability, a significant part of the latter may have been caused by the massive stress of the Italian health system.