We develop a new boundary integral method for solving the coupled electro- and hydrodynamics of vesicle suspensions in Stokes flow. This relies on a well-conditioned boundary integral equation formulation for the leaky-dielectric model describing the electric response of the vesicles and an efficient numerical solver capable of handling highly deflated vesicles. Our method is applied to explore vesicle electrohydrodynamics in three cases. First, we study the classical prolate–oblate–prolate transition dynamics observed upon application of a uniform DC electric field. We discover that, in contrast to the squaring previously found with nearly spherical vesicles, highly deflated vesicles tend to form buds. Second, we illustrate the capabilities of the method by quantifying the electrorheology of a dilute vesicle suspension. Finally, we investigate the pairwise interactions of vesicles and find three different responses when the key parameters are varied: (i) chain formation, where they self-assemble to form a chain that is aligned along the field direction; (ii) circulatory motion, where they rotate about each other; (iii) oscillatory motion, where they form a chain but oscillate about each other. The last two are unique to vesicles and are not observed in the case of other soft particle suspensions such as drops.