f[x_, y_] := x^6 + y^6 + x^2 + y^2
Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2},
MeshFunctions -> {#^2 + #2^2 &},
Mesh -> {{1}},
MeshStyle -> Directive[Red, Thick]]

Update: Keeping the original mesh lines in a single Plot3D
:
Using multiple MeshFunctions
:
Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2},
MeshFunctions -> {#^2 + #2^2 &, # &, #2 &},
Mesh -> {{1}, 20, 20},
MeshStyle -> {Directive[Red, Thick], Gray, Gray}]

Using two functions as the first argument and using the option BoundaryStyle
:
Plot3D[{f[x, y], ConditionalExpression[f[x, y], x^2 + y^2 <= 1]},
{x, -2, 2}, {y, -2, 2}, PlotStyle -> {LightBlue, None},
BoundaryStyle -> {2 -> Directive[Red, Thick]}]

Update 2: to visualize the region $M=\lbrace (x,y) \in \mathbb(R)^2 \mid x^2 + y^2 < 1 \rbrace$:
Plot3D[{ConditionalExpression[f[x, y], x^2 + y^2 >= 1],
ConditionalExpression[f[x, y], x^2 + y^2 <= 1]},
{x, -2, 2}, {y, -2, 2}, PlotStyle -> {LightBlue, Red}]
