The NANOGrav, Parkes and European Pulsar Timing Array (PTA) experiments have collected strong evidence for a stochastic gravitational wave background in the nHz-frequency band. In this work we perform a detailed statistical analysis of the signal in order to elucidate its physical origin. Specifically, we test the standard explanation in terms of supermassive black hole mergers against the prominent alternative explanation in terms of a first-order phase transition. By means of a frequentist hypothesis test we find that the observed gravitational wave spectrum prefers a first-order phase transition at 2-3σ significance compared to black hole mergers (depending on the underlying black hole model). This mild preference is linked to the relatively large amplitude of the observed gravitational wave signal (above the typical expectation of black hole models) and to its spectral shape (which slightly favors the phase-transition spectrum over the predominantly single power-law spectrum predicted in black hole models). The best fit to the combined PTA dataset is obtained for a phase transition which dominantly produces the gravitational wave signal by bubble collisions (rather than by sound waves). The best-fit (energy-density) spectrum features, within the frequency band of the PTA experiments, a crossover from a steeply rising power law (causality tail) to a softly rising power law; the peak frequency then falls slightly above the PTA-measured range. Such a spectrum can be obtained for a strong first-order phase transition in the thick-wall regime of vacuum tunneling which reheats the Universe to a temperature of T∗∼GeV. A dark sector phase transition at the GeV-scale provides a comparably good fit.