Information on the phase structure of strongly interacting matter at high baryon densities can be gained from observations of neutron stars and their detailed analysis. Bayesian inference methods are used to set constraints on the speed of sound in the interior of neutron stars, based on recent multimessenger data in combination with low-density constraints based on chiral effective field theory and perturbative QCD constraints at asymptotically high densities. A systematic Bayes factor assessment quantifies the evidence (or non-evidence) for low averaged sound speeds ($c_s^2 < 0.1$), a prerequisite for a first-order phase transition, within the range of densities realised in the core of neutron stars. One of the consequences of including the recent new heavy ($2.35 M_\odot$) black widow pulsar PSR J0952-0607 in the data base is a further stiffening of the equation-of-state, resulting for a 2.1 solar-mass neutron star in a reduced central density of less than five times the equilibrium density of normal nuclear matter. The evidence against small sound speeds in neutron star cores is further strengthened. Within the inferred posterior credible bands, only a weak first-order phase transition with a coexistence density interval $\Delta n/n \lesssim 0.2$ would be compatible with the observed data. In view of these empirical findings, much discussed issues such as the quest for a phase transition towards restored chiral symmetry and the active degrees of freedom in cold and dense baryonic matter, are explored.