The aim of this paper is to propose and analyze a fractional-order SIS epidemic model with saturating contact rate that is a generalization of a recognized deterministic SIS epidemic model. First, we investigate positivity, boundedness, and asymptotic stability of the proposed fractional-order model. Secondly, we construct positivity-preserving nonstandard finite difference (NSFD) schemes for the model using the Mickens' methodology. We prove theoretically and confirm by numerical simulations that the proposed NSFD schemes are unconditionally positive. Consequently, we obtain NSFD schemes preserving not only the positivity but also essential dynamical properties of the fractional-order model for all finite step sizes. Meanwhile, standard schemes fail to correctly reflect the essential properties of the continuous model for a given finite step size, and therefore, they can generate numerical approximations which are completely different from the solutions of the continuous model. Finally, a set of numerical simulations are performed to support and confirm the validity of theoretical results as well as advantages and superiority of the constructed NSFD schemes. The results indicate that there is a good agreement between the numerical simulations and the theoretical results and the NSFD schemes are appropriate and effective to solve the fractional-order model.