Our aim is to model rare species (with few occurrences) but modelling the distribution of species with few occurrence data and many predictor variables leads to model overfitting. Thus, we use the recently developed ensemble of small models, which showed high predictive accuracy in modelling the distribution of rare species to estimate the current and future distribution of 56 rare (and endangered) saproxylic beetle species. Thus, we stacked predictions from individual species distribution models to derive rare species richness. We used current and five future general circulation models and three representative concentration pathways to test whether the distribution of hotspots for rare species shifts due to climate change under different future scenarios. Moreover, we verified the representativeness of existing protected area systems under future climate conditions in Italy. Specifically, we identified potential hotspots for rare species richness through a cumulative relative frequency distribution function. The current surface covered by hotspots is 50.4% of the study area corresponding to 151,223km(2) (mainly from central to northern Italy). Currently, only 35,124km(2) of rare saproxylic hotspots are covered by protected areas (PAs) and they will decrease by about 2-72% in 2070 depending on the future scenarios considered. Our results confirmed that the shift of the distribution of hotspots for rare species might occur due to climate change under different future scenarios and that the existing PAs system would be inadequate for assuring the conservation of rare saproxylic beetles in Italy under current and future climate conditions.