Summary: The time fractional ODEs are equivalent to convolutional Volterra integral equations with completely monotone kernels. We introduce the concept of complete monotonicity-preserving \((\mathcal{CM}\)-preserving) numerical methods for fractional ODEs, in which the discrete convolutional kernels inherit the \(\mathcal{CM}\) property as the continuous equations. We prove that \(\mathcal{CM}\)-preserving schemes are at least \(A(\pi/2)\) stable and can preserve the monotonicity of solutions to scalar nonlinear autonomous fractional ODEs, both of which are novel. Significantly, by improving a result of L. Li and J.-G. Liu [Q. Appl. Math. 76, No. 1, 189–198 (2018; Zbl 1476.47029)], we show that the \(\mathcal{L}1\) scheme is \(\mathcal{CM}\)-preserving. The good signs of the coefficients for such class of schemes ensure the discrete fractional comparison principles, and allow us to establish the convergence in a unified framework when applied to time fractional sub-diffusion equations and fractional ODEs. The main tools in the analysis are a characterization of convolution inverses for completely monotone sequences and a characterization of completely monotone sequences using Pick functions due to J.-G. Liu and R. L. Pego [Trans. Am. Math. Soc. 368, No. 12, 8499–8518 (2016; Zbl 1350.44003)]. The results for fractional ODEs are extended to \(\mathcal{CM}\)-preserving numerical methods for Volterra integral equations with general completely monotone kernels. Numerical examples are presented to illustrate the main theoretical results.