A technique is developed which yields an asymptotic expansion in the two limits $\lambda \to 0^ + $ and $\lambda \to \infty $ for the fractional integral operator of order $\mu $ with respect to the function $\lambda ^p $ given by $I_{\lambda ^p }^\mu f(\lambda )\frac{1}{{\Gamma (\mu )}}\int_0^\lambda {(\lambda ^p - \xi ^p )^{\mu - 1} p\xi ^{p - 1} f(\xi )d\xi } ,$ under general conditions that f be algebraically dominated near 0 and $\infty $. Representing $I_{\lambda ^p }^\mu f(\lambda )$ as a convolution of Mellin transforms the domain of the transform is extended by analytic continuation. By moving the contour of integration to the right or the left an asymptotic expansion for $I_{\lambda ^p }^\mu f(\lambda )$ can be systematically generated for $\lambda \to \infty $ or $\lambda \to 0^ + $. The technique is illustrated by the asymptotic expansion of fractional integral operators derived from the Euler–Poisson–Darboux equation and generalized axially symmetric potential theory.