NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives

First of all, PDE == nv1 + nv2 + nv3 + nv4 is obviously wrong, because there already exists a == in your PDE. This is easy to fix of course.

What's confusing is the Power::infy warning. I'm not sure why this pops up, maybe NDSolve fails to notice FiniteElement should be chosen in this case, while it should be able to. Anyway, specifying the method explicitly fixes the problem:

solution = NDSolve[{Subtract @@ PDE == nv4, M[0, x, y] == 1}, 
   M, {x, 0, xf}, {y, 0, 1}, {t, 0, 10}, 
   Method -> {"MethodOfLines", "SpatialDiscretization" -> "FiniteElement"}];

Table[ContourPlot[M[t, x, y] /. solution, {x, 0, xf}, {y, 0, 1}, Contours -> 20, 
  ColorFunction -> "TemperatureMap", PlotLegends -> Automatic, 
  PlotLabel -> Row[{"t = ", t}]], {t, 1, 10, 3}]

nv1, nv2, nv3 are omitted because zero Neumann value is the default setting of FiniteElement.

Your equation is not in the correct form for FEM. Homogeneous Neumann conditions are applied automatically.

Needs["NDSolve`FEM`"]
Bi = 0.5; xf = 5; reg = Rectangle[{0, 0}, {xf, 1}];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.0001];
PDE = D[M[t, x, y], t] - 1/(x^2 + y^2)*D[(x^2 + 1)*D[M[t, x, y], x], x] - 
 1/(x^2 + y^2)*D[(-y^2 + 1)*D[M[t, x, y], y], y]

nv4 = NeumannValue[-Sqrt[(x^2 + y^2)/(x^2 + 1)]*Bi*M[t, x, y], 
   x == xf];



sol = NDSolve[{PDE == nv4, M[0, x, y] == 1}, 
  M, {x, y} \[Element] mesh, {t, 0, 10}]

Table[ContourPlot[M[t, x, y] /. sol, {x, y} \[Element] mesh, 
  Contours -> 20, ColorFunction -> "TemperatureMap", 
  PlotLegends -> Automatic, PlotLabel -> Row[{"t = ", t}]], {t, 1, 10,
   3}]

Consider a solution in a region with a hole (homogeneous Neumann conditions are applied automatically at the edge of the hole)

Needs["NDSolve`FEM`"]
Bi = 0.5; xf = 5; reg = Rectangle[{0, 0}, {xf, 1}]; reg1 = 
 Rectangle[{1, 1/3}, {4, 2/3}];
mesh = ToElementMesh[RegionDifference[reg, reg1], 
   MaxCellMeasure -> 0.0001];
PDE = D[M[t, x, y], t] - 
   1/(x^2 + y^2)*D[(x^2 + 1)*D[M[t, x, y], x], x] - 
   1/(x^2 + y^2)*D[(-y^2 + 1)*D[M[t, x, y], y], y];

nv4 = NeumannValue[-Sqrt[(x^2 + y^2)/(x^2 + 1)]*Bi*M[t, x, y], 
   x == xf];

sol = NDSolve[{PDE == nv4, M[0, x, y] == 1}, 
  M, {x, y} \[Element] mesh, {t, 0, 10}]

Table[ContourPlot[M[t, x, y] /. sol, {x, y} \[Element] mesh, 
  Contours -> 20, ColorFunction -> "TemperatureMap", 
  PlotLegends -> Automatic, PlotLabel -> Row[{"t = ", t}]], {t, 1, 10,
   3}]