As a critical hydrological cycle parameter, reference evapotranspiration (${\mathit{ET}}_0$) plays a vital role in guiding agricultural water management. Analyzing the multifractal nature of ${\mathit{ET}}_0$ can clarify the complex fluctuations and obtain insights into the evolutionary mechanism and the regular patterns of the time series. The present study used the Multifractal Detrended Fluctuation Analysis (MFDFA) method to explore the multifractality of daily ${\mathit{ET}}_0$ series from 1961 to 2017, which filtered the seasonal component at 719 weather stations in seven sub-regions of China. The results indicated that the ${\mathit{ET}}_0$ series exhibited persistence and presented a decreasing trend in all sub-regions as the value of h(2) was exceeded 0.5 and the right-skewed was displayed in the multifractal spectrum. Multifractality existed and showed heterogeneity among various sub-regions, where NE was characterized by the strongest multifractality, followed by NW, IM, CSC, QT, SC, and NC. Both long-range correlation and broad probability density function contributed to the multifractality, and the long-range correlation played a dominant role in multifractality. This study is helpful to understand the complex fluctuation of ${\mathit{ET}}_0$ series, and provides new ideas for accurately predicting agricultural water demand and formulating scientific irrigation schedules in different climate zones of China.