The 3D velocities of M31 and M33 are important for understanding the evolution and cosmological context of the Local Group . Their most massive stars are detected by Gaia , and we use Data Release 2 ( DR2 ) to determine the galaxy proper motions ( PMs ) . We select galaxy members based on , e.g. , parallax , PM , color-magnitude-diagram location , and local stellar density . The PM rotation of both galaxies is confidently detected , consistent with the known line-of-sight rotation curves : V _ { rot } = -206 \pm 86 \ > { km } { s } ^ { -1 } ( counter-clockwise ) for M31 , and V _ { rot } = 80 \pm 52 \ > { km } { s } ^ { -1 } ( clockwise ) for M33 . We measure the center-of-mass PM of each galaxy relative to surrounding background quasars in DR2 . This yields that ( \mu _ { \alpha* } , \mu _ { \delta } ) equals ( 65 \pm 18 , -57 \pm 15 ) \ > \mu { as } { yr } ^ { -1 } for M31 , and ( 31 \pm 19 , -29 \pm 16 ) \ > \mu { as } { yr } ^ { -1 } for M33 . In addition to the listed random errors , each component has an additional residual systematic error of 16 \ > \mu { as } { yr } ^ { -1 } . These results are consistent at 0.8 \sigma and 1.0 \sigma with the ( 2 and 3 times higher-accuracy ) measurements already available from Hubble Space Telescope ( HST ) optical imaging and VLBA water maser observations , respectively . This lends confidence that all these measurements are robust . The new results imply that the M31 orbit towards the Milky Way is somewhat less radial than previously inferred , V _ { tan,DR 2 + HST } = 57 ^ { +35 } _ { -31 } \ > { km } { s } ^ { -1 } , and strengthen arguments that M33 may be on its first infall into M31 . The results highlight the future potential of Gaia for PM studies beyond the Milky Way satellite system .