Using a data sample of proton-proton collisions at √s=13  TeV, corresponding to an integrated luminosity of 140  fb−1 collected by the CMS experiment in 2016–2018, the B0s→X(3872)ϕ decay is observed. Decays into J/ψπ+π− and K+K− are used to reconstruct, respectively, the X(3872) and ϕ. The ratio of the product of branching fractions B[B0s→X(3872)ϕ]B[X(3872)→J/ψπ+π−] to the product B[B0s→ψ(2S)ϕ]B[ψ(2S)→J/ψπ+π−] is measured to be [2.21±0.29(stat)±0.17(syst)]%. The ratio B[B0s→X(3872)ϕ]/B[B0→X(3872)K0] is found to be consistent with one, while the ratio B[B0s→X(3872)ϕ]/B[B+→X(3872)K+] is two times smaller. This suggests a difference in the production dynamics of the X(3872) in B0 and B0s meson decays compared to B+. The reported observation may shed new light on the nature of the X(3872) particle.