When DeepONet approximates solution operators of partial differential equations (PDEs) with discontinuous solutions, it poses a foundational approximation lower bound due to its linear reconstruction property. Inspired by the moving mesh method, we propose an R-adaptive DeepONet method, which consists of: (1) the output data representation is transformed from the physical domain to the computational domain using the equidistribution principle; (2) the maps from input parameters to the solution and the coordinate transformation function over the computational domain are learned using DeepONets separately; (3) the solution over the physical domain is obtained via post-processing methods such as the interpolation method. Additionally, we introduce a solution-dependent weighting strategy in the training process to reduce the error. We establish an upper bound for the reconstruction error based on piecewise linear interpolation and show that the introduced R-adaptive DeepONet can reduce this bound. Moreover, for two prototypical PDEs with sharp gradients or discontinuities, we prove that the approximation error decays at a superlinear rate with respect to the trunk basis size, unlike the linear decay observed in vanilla DeepONets. Numerical experiments on several PDEs with discontinuous solutions are conducted to verify the advantages of the R-adaptive DeepONet over available variants of DeepONet.