Turbulence is widely expected to limit the confinement, and thus the overall performance, of modern neoclassically-optimized stellarators. We employ novel petaflop-scale gyrokinetic simulations to predict the distribution of turbulence fluctuations and the related transport scaling on entire stellarator magnetic surfaces, and reveal striking differences to tokamaks. Using a stochastic global-search optimization method, we derive the first turbulence-optimized stellarator configuration stemming from an existing quasi-omnigenous design.Throughout the history of magnetic fusion, a recurrent theme has been the surprising sensitivity of plasma performance to the details of the magnetic field. For instance, it has long been known that the confinement of alpha particles can be spoiled by small ripples in the magnetic field. More recently, magnetic perturbations have been found to dramatically influence instabilities of the plasma edge [1]. In both stellarators and tokamaks, experiments show that the level of turbulence may be reduced by modifying the magnetic field. As notable examples, the confinement time in the TCV tokamak is doubled by reversing the triangularity of the poloidal cross section of the flux surfaces [2], and in the LHD stellarator the turbulent transport can be reduced significantly by adjusting the coil currents so as to shrink the circumference of the torus by pushing it radially inwards [3].Stellarators typically possess 40-50 degrees of freedom in the shaping of the magnetic field, an order of magnitude more than for tokamaks [4]. This enhanced flexibility can be used to optimize various plasma properties [5], and the latest demonstration of the power of such optimization is expected to be realized in the superconducting stellarator experiment Wendelstein 7-X (W7-X), in Greifswald, Germany [6]. A tantalizing possibility for fusion researchers is to try to exploit any leeway in the magnetic geometry to design configurations with better confinement. In W7-X, this has already been done for the collisional (so-called "neoclassical") transport, but no device built so far is optimized with respect to turbulence.In order to understand how energy transport depends on the magnetic-field geometry, it is helpful to numerically simulate the turbulence in a large portion of the plasma. In tokamaks, thanks to axisymmetry, restricting the computational domain to a "flux tube", a slender volume encompassing a magnetic-field line [7], suffices to calculate the transport at a radial location. In a stellarator, however, different flux tubes on a magnetic surface are not geometrically equivalent, thus it appears necessary to simulate the entire magnetic surface. Much has been learned from the flux-tube approach which, except for inspiring efforts [8], has characterized stellarator turbulence simulations to date [9][10][11][12] . Nevertheless, all these simulations have a major inherent drawback: the transport cannot be reliably determined, as the turbulence strength generally varies between different flux tubes on th...