We analyze the newest data from the NA61/SHINE collaboration which, in addition to previous results on pions and kaons, include mean multiplicities of p, Λ, and ϕ-mesons produced in inelastic proton-proton (p+p) interactions at $\sqrt{s}$ = 6.3−17.3 GeV. The canonical ensemble formulation of the ideal hadron resonance gas (HRG) model is used with exact conservation of net baryon number B=2, electric charge Q=2, and strangeness S=0. The chemical freeze-out parameters in p+p interactions are obtained and compared to those in central nucleus-nucleus collisions. We find several serious problems of the statistical model treatment of p+p interactions at the CERN SPS: 1) the inclusion of the ϕ-meson yields in thermal fits worsens significantly the fit quality (the ϕ-meson puzzle); 2) the data shows large event-by-event multiplicity fluctuations in inelastic p+p interactions which cannot be explained within a statistical model, and provide evidence for large non-statistical fluctuations; 3) the fits within the canonical ensemble formulation of HRG do not give any improvement over the fits within the grand canonical ensemble formulation, i.e., there are no indications for existence of a single statistical system in p+p interactions.