The monotonicity of the Mittag-Leffler function $E_α$ with respect to the parameter $α$ is investigated, via some convex ordering properties for related random variables. In particular, it is shown that the mapping $α\mapsto E_α(x^α)$ decreases on $(0,2)$ for all $x> 0$, that the mapping $α\mapsto E_α(-x^α)$ decreases on $(0,1)$ for all $x\ge 1$ and that the mapping $α\mapsto E_α(Γ(1+α)x)$ decreases on $(0,1)$ for all $x\in{\mathbb R}^\ast.$ Analogous results are presented for the two parameter Mittag-Leffler functions $E_{α, β}$ with $β\ge α,$ with an emphasis on the extremal case $β=α.$ Several applications of these results are discussed for Abelian integral equations and subdiffusions.