Understand MeshFunctions with Intercepts

It looks like a bug. However there is a simple workaround: change the sign of your mesh function. In addition, p[x] /. x -> # & is equivalent to #2 &.

Plot[p[x], {x, -3, 4}, PlotStyle -> Green, MeshFunctions -> {-#2 &}, 
 Mesh -> {{0}}, MeshStyle -> Directive[PointSize[Large], Red]]

A very good question. You should not rely on MeshFunctions.

Follow the old path step for step:

fun = x^4 - 4 x^2 - 2 x + 1;

xvals = NSolve[fun == 0, x] // Values // Flatten

{-1.48119, -1., 0.311108, 2.17009}

yvals = fun /. xvals // Chop

{0, 0, 0, 0}

zeroes = Transpose[{xvals, yvals}]

{{-1.48119, 0}, {-1., 0}, {0.311108, 0}, {2.17009, 0}}

snippet1 =
 Plot[fun, {x, -3, 5},
  Epilog -> {Red, PointSize[0.02], Point /@ zeroes},
  PlotStyle -> Green,
  PlotLabel -> "Snippet 1"]