Direct numerical simulations of the motion of 27 three-dimensional deformable buoyant bubbles in periodic domains are presented. The full Navier–Stokes equations are solved by a parallelized finite-difference/front-tracking method that allows a deformable interface between the bubbles and the suspending fluid and the inclusion of surface tension. The Eötvös number is taken as equal to 5, so that the bubbles are ellipsoidal, and the Galileo number is 900, so that the rise Reynolds number of a single bubble in an unbounded flow is about 26. Three values of the void fraction have been investigated: 2%, 6% and 12%. At 6%, a change in the behaviour of the bubbles is observed. The bubbles are initially dispersed homogeneously throughout the flow field and their average rise Reynolds number is 23. After the bubbles have risen by about 90 bubble diameters, they form a vertical stream and accelerate. The microstructure of the bubble suspension is analysed and an explanation is proposed for the formation of these streams. The results for the ellipsoidal bubbles are compared to the results for nearly spherical bubbles, for which the Eötvös number is 1 and the Galileo number is 900. The dispersion of the bubbles and the velocity fluctuations in the liquid phase are analysed.