The majority of problems in aircraft production and operation require decisions made in the presence of uncertainty. For this reason, aerodynamic designs obtained with traditional deterministic optimization techniques seeking only optimality in a specific set of conditions may have very poor off-design performances or may even be unreliable. In this work, a novel approach for robust and reliability-based design optimization of aerodynamic shapes based on the combination of single- and multi-objective evolutionary algorithms and a continuation multilevel Monte Carlo methodology is presented, to compute objective functions and constraints that involve statistical moments or statistical quantities, such as quantiles, also called value at risk and conditional value at risk, without relying on derivatives and meta-models. Detailed numerical studies are presented for the RAE 2822 transonic airfoil design affected by geometrical and operational uncertainties.