A chemo-mechanical model is used to capture the formation and evolution of microdomains on the deforming surface of giant unilamellar vesicles. The model is intended for the regime of vesicle dynamics characterized by a distinct difference in time scales between shape change and species transport. This is achieved by ensuring that shape equilibrium holds away from chemical equilibrium. Conventional descriptions are used to define the curvature and chemical contributions to the vesicle energetics. Both contributions are consistently non-dimensionalized. The phase-field framework is used to cast the coupled model in a diffuse-interface form. The resulting fourth-order nonlinear system of equations is discretized using the finite- element method with a uniform cubic spline basis, which satisfies global higher-order continuity. Two-dimensional and axisymmetric numerical examples of domain evolution coupled to vesicle shape deformation are presented. Curvature-dependent domain sorting and shape deformation dominated by line tension are also considered.