Python3, 821 bytes

Long, but computes all test cases in < 20 seconds on TIO.

E=enumerate
def P(x,y,X,Y,D):
 q=[(x,y,[(x,y)])]
 for x,y,p in q:
  if(x,y)==(X,Y):return p
  for J,K in[(1,0),(0,1),(-1,0),(0,-1),(1,1),(-1,-1),(1,-1),(-1,1)]:
   if D.get(V:=(x+J,y+K),4)<2 and V not in p:q+=[(*V,p+[V])]
def f(d):
 K=[[0]*(len(d[0])+2)]
 q,C=[(D:={(x,y):v for x,r in E(K+[[0]+i+[0]for i in d]+[*map(list,K)])for y,v in E(r)},[[i]for i in D if 1==D[i]])],{}
 while q:
  q=sorted(q,key=lambda x:sum(map(len,x[1])))
  D,M=q.pop(0)
  if len(M)==1:return len(M[0])-sum(D[i]==1 for i in D)
  for a,b in E(M):
   for A,B in E(M):
    if A>a:
     Z=[]
     for t in b:
      for T in B:
       if(t,T)not in C:C[(t,T)]=P(*t,*T,D);C[(T,t)]=C[(t,T)]
       Z+=[(t,T)]
     L=C[min(Z,key=lambda x:len(C[x]))]
     Q,W=[],[]
     for J in M:
      if{*J}&{*L}:Q+=J
      else:W+=[J]
     q+=[(D,[I:=[*{*Q,*L}]]+W)]

Try it online!

05AB1E, 63 54 52 bytes

Quite slow, but after some optimisation, this can now run some testcases!

2Føðδ.ø}©˜„ CSδQƶø0δK<®øg‰`Uæé.ΔX«Dδ-nO3‹DvDδ*O}ßĀ}g

Try it online!

The algorithm used has complexity \$\mathcal{O}\left(2^{(n+2)^2} (n+2)^6\right)\$ (for a grid of size \$n \times n\$) and uses matrix multiplication (øδ*O), which is often incredibly slow in 05AB1E.

How?

A slightly different interpretation of the tasks helps explaining the approach: Instead of placing roads, we can also place additional cities into the grid until all cities are connected.
The program tries all combinations of new cities in increasing order of new cities until it finds a valid solution.

To validate a given attempt we consider a graph where the vertices are cities and two cities are connected with an edge if they are vertically, horizontally or diagonally adjacent. Two cities are adjacent, iff their squared Euclidean distance is less than \$3\$.
This property can be used to calculate an adjacency matrix from a list of city coordinates. An undirected graph of \$n\$ vertices with adjacency matrix \$A\$ is connected, iff \$A^n\$ has no zeros. Example:

map
 C
C C

city coordinates
[0,1], [1,0], [1,2]

pairwise squared Euclidean distances
0 2 2
2 0 4
2 4 0

adjacency matrix A
1 1 1
1 1 0
1 0 1

A^3
7 5 5
5 4 3
5 3 4

min(A^3) = 3 > 0 => connected

Try your own example!

Commented Code:

                         # pad with spaces
2F     }                 # execute two times:
  ø                      #   transpose the grid
   ðδ.ø                  #   and surround every row with spaces


                         # get coordinates of cities and spaces
©                        # store the padded grid in the register
 ˜                       # flatten the grid
  „ CS                   # push [" ", "C"]
      δQ                 # equality table (== [[char==" ", char=="C"] for char in flat_grid])
        ƶ                # multiply each pair by its 1-based index
         ø               # tranpose
          0δK            # remove all zeros
             <           # decrement by 1
                         # Now we have a two list of coordinates of spaces and cities in the flattened grid
              ®          # get the grid from the register
               øg        # get the height (length of the tranpose)
                 ‰       # divmod each flat index with this
                         # this results into the 2d-indices
                  `      # push cities and spaces seperately on the stack
                   U     # store coordinates of cities in variable X
                         # coordinates of spaces are now on top of the stack 

                         # try subsets of the spaces until we get a connected graph
æ                        # push the powerset (all subsets) of space coordinates
 é                       # sort by length
  .Δ                  }  # find the first value for which the following is truthy
    X«                   # append the existing cities to the current subset of spaces
      D                  # duplicate
       δ-                # subtraction table
         n               # square each difference
          O              # sum each pair of squared differences
           3‹            # for each number: is it less than 3?
                         # this is the adjacency matrix A
             D           # make a copy of the matrix
              v    }     # for each row in A:
               D         # duplicate current matrix
                δ*O      # matrix multiplication (for symmetric matrices)
                    ßĀ   # is the minimum not 0?

g                        # take the length of the result

05AB1E, 31 bytes

Takes input as lists of 1-indexed coordinates of cities and mountains. Even less efficient if the input is not square.

«Z>ÝãsKæé.Δ¹«Dδ-nO3‹DvDδ*O}ßĀ}g

Try it online!

«        # concatenate cities and mountains
 Z>      # maximum integer in the list + 1 (0 and +1 for padding)
   Ý     # push range from 0 to this value
    ã    # cartesian square: all 2d coordinates with values in that range
     sK  # remove all cities and mountains, this leaves the coordinates of spaces
 ...     # same as above

Jelly, 29 bytes

;_þ`²§<3æ*L$Ȧ
FṀ‘Żṗ2ḟẎŒPçƇḢḢL

Try it online!

Uses ovs' method, be sure to give them an upvote!

Takes the input as a pair [a, b] where a is a list of 1 indexed co-ordinates of cities and b a list of 1 indexed co-ordinates of mountains

This is stupidly slow. For an \$n \times m\$ matrix with \$x\$ spaces, this has a lower bound of \$O(2^{\max(n,m)^2-x})\$, and that doesn't factor in the speed of Jelly's matrix power, which isn't fast. This doesn't finish on TIO for any valid input

Jelly, 36 bytes

;’_þ`²§<3æ*L$Ȧ
Ż€ZUƊ4¡µœẹ0ŒPçƇœẹ1$ḢL

Try it online!

Also incredibly slow. For an input of size \$n\times m\$ with \$x\$ spaces, this has an absolute lower bound of \$O(2^{2(n+m)+x})\$, and that doesn't factor in the speed of Jelly's matrix power, or the fact that it has to calculate \$2^{2(n+m)+x}\$ Cartesian products with the number of cities in the input.

There are various optimisations you can do if you know the input before hand. For example, this is a much faster version that works for inputs that can be constructed as going over the "top" of the mountains.

Jelly, 34 bytes

;’_þ`²§<3æ*L$Ȧ
Ż€ZUƊ4¡µœẹ0ŒPçƇŒṪḢL

Try it online!

Exactly the same as the 36 byte version, but takes input as a matrix where 0 is a space, 1 is a city and [] is a mountain. For example, C M C is [[1, 0, [], 0, 1]]

How they work

;_þ`²§<3æ*L$Ȧ - Helper link. Takes S (indices of spaces) and C (indices of cities)
;             - Concatenate
 _þ`          - Create a subtraction table
    ²         - Square each
     §        - Sums of each
      <3      - Less than 3?
           $  - To this square matrix M:
          L   -   Get its length
        æ*    -   Raise it to that power
            Ȧ - Are all elements of the matrix non-zero?

FṀ‘Żṗ2ḟẎŒPçƇḢḢL - Main link. Takes [a, b] on the left
F               - Flatten
 Ṁ              - Maximum
  ‘             - Increment
   Ż            - Range from zero
    ṗ2          - Cartesian square
                   This gets all possible coordinates, including an outer border
       Ẏ        - Flatten the argument into a list of coords
      ḟ         - Remove the mountain + city coords from the list
        ŒP      - Powerset
            Ḣ   - Extract the list of city coords
          çƇ    - Keep the powerset elements that are truthy under the helper link
             ḢL - Take the shortest one and return its length

;’_þ`²§<3æ*L$Ȧ - Helper link. Takes S (indices of spaces) and C (indices of cities)
;              - Concatenate
 ’             - Decrement to zero index
  _þ`          - Create a subtraction table
     ²         - Square each
      §        - Sums of each
       <3      - Less than 3?
            $  - To this square matrix M:
           L   -   Get its length
         æ*    -   Raise it to that power
             Ȧ - Are all elements of the matrix non-zero?

Ż€ZUƊ4¡µœẹ0ŒPçƇœẹ1$ḢL - Main link. Takes matrix M on the left
    Ɗ4¡               - Do the following 4 times:
Ż€                    -   Prepend a zero to each row
  Z                   -   Transpose
   U                  -   Reverse
       µ              - Begin a new link with this bordered matrix M' as the argument
        œẹ0           - Get the indices of zeros
           ŒP         - Powerset
                  $   - To M':
               œẹ1    -   Get the indices of ones
             çƇ       - Keep the powerset elements that are truthy under the helper link
                   ḢL - Take the shortest one and return its length

The third one takes advantage of the fact [] is both falsey and non-zero, and replaces œẹ1$ with ŒṪ which gets the indices of all truthy elements in M'

Python 3 + numpy, 325 bytes

A (slightly faster) port of my 05AB1E solution. The computational complexity is still the same, but especially the matrix operations are faster.

from itertools import*
from numpy import*
def f(G):
 E=enumerate;k=[' '];p,c,*s=[k*len(G[0])],[]
 for y,r in E([k+r+k for r in p+G+p]):
  for x,v in E(r):[s,p,p,c][ord(v)%4]+=[y,x],
 for i in count():
  for N in combinations(s,i):
   v=array(c+[*N]);A=sum((v-v[:,None])**2,axis=2)<3
   for _ in v:A=A@A
   if A.min():return i

Try it online!

Ungolfed full program:

import itertools, sys
import numpy as np

grid = [list(line) for line in sys.stdin.read().splitlines()]
print(grid)
width = len(grid[0])

# pad with spaces on all sides
padding = [[' '] * width]
grid = [[' ', *row, ' '] for row in padding + grid + padding]

# get the coordinates of cities and empty spaces
cities = [[y, x] for y, row in enumerate(grid) for x, v in enumerate(row) if 'C' == v]
spaces = [[y, x] for y, row in enumerate(grid) for x, v in enumerate(row) if ' ' == v]

for i in range(len(spaces)):
  # Iterate over all combinations of i new cities
  for new_cities in itertools.combinations(spaces, i):
    vertices = np.array([*cities, *new_cities])
    adjacency = (np.sum((vertices - vertices[:, None])**2, axis=-1)<3)
    if (np.linalg.matrix_power(adjacency, adjacency.shape[0]) > 0).all():
      print(i)
      exit()

Try it online!

Kotlin with JGraphT, 3009 2984 (thanks to ophact) bytes (plus header/footer)

import org.jgrapht.Graph
import org.jgrapht.GraphPath
import org.jgrapht.alg.connectivity.ConnectivityInspector
import org.jgrapht.alg.shortestpath.AllDirectedPaths
import org.jgrapht.alg.shortestpath.DijkstraManyToManyShortestPaths
import org.jgrapht.graph.DefaultDirectedGraph
import org.jgrapht.graph.DefaultEdge
import org.jgrapht.graph.DefaultUndirectedGraph
import org.jgrapht.graph.builder.GraphBuilder
import kotlin.streams.asStream

sealed class P
object C:P()
object M:P()
object E:P()
object R:P()
typealias D=Pair<Int,Int>
typealias L<T> =List<T>
typealias A<T> =Array<T>
typealias F=DefaultEdge
typealias G=BooleanArray
object H:ThreadLocal<G>(){override fun initialValue()=G(1 shl(B*2)){false}}
const val B=4
fun GraphPath<D,F>.v():A<D>{val r=arrayOfNulls<D>(edgeList.size+1)
r[0]=graph.getEdgeSource(edgeList.first())
edgeList.forEachIndexed{i,e->r[i+1]=graph.getEdgeTarget(e)}
return r as A<D>}
fun A<GraphPath<D,F>>.d()=sumOf{val table=if(B<=3){G(1 shl(B*2)){false}}else{H.get().also{for(i in it.indices){it[i]=false}}}
it.v().sumOf{val hash=(it.first+1)or(it.second+1).shl(B)
if(table[hash]){0}else{table[hash]=true
1}.toInt()}}
fun j(b:GraphBuilder<D,*,*>,x:Int,y:Int){val others=(x-1..x+1).flatMap{z->(y-1..y+1).map{z to it}}
others.forEach{b.addEdge(x to y,it)}}
fun g(g:L<L<P>>):Pair<Set<D>,Graph<D,F>>{val b=DefaultDirectedGraph.createBuilder<D,F>(::DefaultEdge)
val c=mutableSetOf<D>()
g.forEachIndexed{y,r->r.forEachIndexed{x,p->when(p){C->{c+=x to y
j(b,x,y)}
E->j(b,x,y)
R->error("")}}}
val width=g[0].size
g.indices.forEach{y->if (g[y][0]==M)j(b,-1,y)
if (g[y].last()==M)j(b,width,y)}
(0 until width).forEach{x->if (g[0][x]==M)j(b,x,-1)
if (g.last()[x]==M)j(b,x,g.size)}
return c to b.buildAsUnmodifiable()}
fun t(c: Set<D>, g:Graph<D,F>):Set<D>{val d=DijkstraManyToManyShortestPaths(g).getManyToManyPaths(c,c)
val h=c.flatMap{o->c.map{o to it}}.filter{if(it.first.first==it.second.first)it.first.second>it.second.second else it.first.first>it.second.first}.associateWith{d.getPath(it.first,it.second)}.map{it.key to AllDirectedPaths(g).getAllPaths(it.key.first,it.key.second,true,it.value.length)}
val j=h.map{it.first}
val k=h.map{it.second}
val l=k.s{val m=DefaultUndirectedGraph.createBuilder<D,F>(::DefaultEdge)
m.addVertices(*c.toTypedArray())
it.mapIndexed{i,b->if(b)j[i]else null}.filterNotNull().forEach{m.addEdge(it.first,it.second)}
ConnectivityInspector(m.buildAsUnmodifiable()).isConnected}
return HashSet(l.asStream().parallel().map{it to it.d()}.min{n,o->n.second-o.second}.get().first.flatMap{it.v().asList()})}
fun a(a:String)=a.lines().map{it.map{when(it){'M'->M
'.',' '->E
'C'->C
'R'->R
else->null}}.filterNotNull()}.filter{it.isNotEmpty()}
fun l(data:L<L<P>>, route:Set<D>)=route.sumOf{if(it.second<0||it.second>data.lastIndex)return@sumOf 1
val row=data[it.second]
if(it.first<0||it.first>row.lastIndex)return@sumOf 1
if(row[it.first]==E)1.toInt()else 0}
inline fun<reified T>L<L<T>>.s(noinline test:(A<Boolean>)->Boolean)=List(size){listOf(false,true)}.c().filter(test).let{sequence{for(selection in it){yieldAll(filterIndexed { i, _->selection[i]}.c())}}}
inline fun <reified T>L<L<T>>.c()=sequence{if(isEmpty())return@sequence
val j=map{it.iterator()}.toTypedArray()
val v=j.map{it.next()}.toTypedArray()
yield(v.clone())
while(true){var i=0
while(!j[i].hasNext()){j[i]=get(i).iterator()
v[i]=j[i].next()
if(++i>lastIndex)return@sequence}
v[i]=j[i].next()
yield(v.clone())}}

fun main() {
    val art = System.`in`.bufferedReader().readText()
    val data = a(art)
    val (cities, graph) = g(data)
    val route = t(cities, graph)
    println("Length: ${l(data, route)}")
}

Here is the original, prior to golfing:

import org.jgrapht.Graph
import org.jgrapht.GraphPath
import org.jgrapht.alg.connectivity.ConnectivityInspector
import org.jgrapht.alg.shortestpath.AllDirectedPaths
import org.jgrapht.alg.shortestpath.DijkstraManyToManyShortestPaths
import org.jgrapht.graph.DefaultDirectedGraph
import org.jgrapht.graph.DefaultEdge
import org.jgrapht.graph.DefaultUndirectedGraph
import org.jgrapht.graph.builder.GraphBuilder
import java.util.concurrent.atomic.AtomicLong
import kotlin.concurrent.thread
import kotlin.streams.asStream

sealed class Place
object City : Place()
object Mountain : Place()
object Empty : Place()
object Road : Place()

typealias Coord = Pair<Int, Int>

operator fun Coord.compareTo(other: Coord): Int {
    return first.compareTo(other.first).let {
        if (it == 0)
            second.compareTo(other.second)
        else
            it
    }
}

fun <T> Collection<T>.toHashSet() = HashSet(this)

object HashTable : ThreadLocal<BooleanArray>() {
    override fun initialValue(): BooleanArray {
        return BooleanArray(1 shl (MAX_DIMEN_BITS * 2)) { false }
    }
}

const val MAX_DIMEN_BITS = 4

inline fun <reified V, E> GraphPath<V, E>.getVertexListFast(): Array<V> {
    val ret = Array<V?>(edgeList.size + 1) { null }
    ret[0] = graph.getEdgeSource(edgeList.first())
    edgeList.forEachIndexed { i, edge ->
        ret[i + 1] = graph.getEdgeTarget(edge)
    }
    @Suppress("UNCHECKED_CAST")
    return ret as Array<V>
}

fun Array<GraphPath<Coord, DefaultEdge>>.distinctSize() = sumOf { path ->
    val table = if (MAX_DIMEN_BITS <= 3) {
        BooleanArray(1 shl (MAX_DIMEN_BITS * 2)) { false }
    } else {
        HashTable.get().also {
            for (i in it.indices) {
                it[i] = false
            }
        }
    }
    path.getVertexListFast().sumOf {
        val hash = (it.first + 1) or (it.second + 1).shl(MAX_DIMEN_BITS)
        if (table[hash]) {
            0
        } else {
            table[hash] = true
            1
        }.toInt() /* stupid overload resolution */
    }
}

fun joinUp(builder: GraphBuilder<Coord, *, *>, x: Int, y: Int) {
    val others = (x - 1 .. x + 1).flatMap { otherX -> (y - 1 .. y + 1).map { otherY -> otherX to otherY } }
    others.forEach {
        builder.addEdge(x to y, it)
    }
}

fun createGraph(grid: List<List<Place>>): Pair<Set<Coord>, Graph<Coord, DefaultEdge>> {
    val builder = DefaultDirectedGraph.createBuilder<Coord, DefaultEdge> { DefaultEdge() }
    val cities = mutableSetOf<Coord>()
    grid.forEachIndexed { y, row ->
        row.forEachIndexed { x, place ->
            when (place) {
                City -> {
                    cities += x to y
                    joinUp(builder, x, y)
                }
                Mountain -> {
                }
                Empty -> joinUp(builder, x, y)
                Road -> error("Cannot pre-supply Roads")
            }
        }
    }
    val width = grid.first().size
    grid.indices.forEach { y ->
        if (grid[y].first() == Mountain)
            joinUp(builder, -1, y)
        if (grid[y].last() == Mountain)
            joinUp(builder, width, y)
    }
    (0 until width).forEach { x ->
        if (grid.first()[x] == Mountain)
            joinUp(builder, x, -1)
        if (grid.last()[x] == Mountain)
            joinUp(builder, x, grid.size)
    }
    return cities to builder.buildAsUnmodifiable()
}

fun joinCities(cities: Set<Coord>, graph: Graph<Coord, DefaultEdge>): Set<Coord> {
    val allPathsAlgo = AllDirectedPaths(graph)
    val dijkstraAlgo = DijkstraManyToManyShortestPaths(graph)
    val shortestPaths = dijkstraAlgo.getManyToManyPaths(cities, cities)
    val sourceTargets = cities.flatMap { outer -> cities.map { outer to it } }.filter { it.first > it.second }
    val shortestLengths = sourceTargets.associateWith {
        shortestPaths.getPath(it.first, it.second)
    }
    val allShortestPathsWithKey = shortestLengths.map {
        it.key to allPathsAlgo.getAllPaths(it.key.first, it.key.second, true, it.value.length)
    }
    val allShortestPathsKeys = allShortestPathsWithKey.map { it.first }
    val allShortestPaths = allShortestPathsWithKey.map { it.second }
    val i = AtomicLong()
    val (total, allPathsProduct) = allShortestPaths.selectiveCartesianProduct { selection ->
        val builder = DefaultUndirectedGraph.createBuilder<Coord, DefaultEdge> { DefaultEdge() }
        builder.addVertices(*cities.toTypedArray())
        selection.mapIndexed { i, bool -> if (bool) allShortestPathsKeys[i] else null }.filterNotNull().forEach { builder.addEdge(it.first, it.second) }
        val selectionGraph = builder.buildAsUnmodifiable()
        ConnectivityInspector(selectionGraph).isConnected
    }
    val startTime = System.currentTimeMillis()
    val counterThread = thread(true) {
        try {
            while (true) {
                Thread.sleep(1000)
                val progress = i.getOpaque()
                val percent = (progress.toDouble() / total.toDouble()) * 100f
                print("Progress: $progress\t/ $total\t($percent%)\t(${System.currentTimeMillis() - startTime} ms)\r")
            }
        } catch (e: InterruptedException) {
            println()
        }
    }
    val ret = allPathsProduct
        .asStream()
        .parallel()
        .map { (it to it.distinctSize()).also { i.incrementAndGet() } }
        .min { left, right -> left.second - right.second }.get().first.flatMap { it.getVertexListFast().asIterable() }.toHashSet()
    counterThread.interrupt()
    counterThread.join()
    return ret
}

fun artToData(art: String) = art.lines().map {
    it.map { char ->
        when (char) {
            'M' -> Mountain
            '.', ' ' -> Empty
            'C' -> City
            'R' -> Road
            else -> null
        }
    }.filterNotNull()
}.filter { it.isNotEmpty() }

fun plotRoute(data: List<List<Place>>, route: Set<Coord>): List<List<Place>> {
    val width = data[0].size
    val emptyLine = List<Place>(width + 2) { Empty }
    val newData = mutableListOf(emptyLine.toMutableList(), *data.map { mutableListOf(Empty, *it.toTypedArray(), Empty) }.toTypedArray(), emptyLine.toMutableList())
    for (point in route) {
        if (newData[point.second + 1][point.first + 1] == Empty) {
            newData[point.second + 1][point.first + 1] = Road
        }
    }
    return newData
}

fun routeLength(data: List<List<Place>>, route: Set<Coord>): Int = route.sumOf { point ->
    if (point.second < 0 || point.second > data.lastIndex) return@sumOf 1
    val row = data[point.second]
    if (point.first < 0 || point.first > row.lastIndex) return@sumOf 1
    if (row[point.first] == Empty) 1.toInt() /* needed due to overload ambiguity */ else 0
}

fun dataToArt(data: List<List<Place>>) = buildString {
    data.forEach { row ->
        row.forEach { place ->
            append(
                when (place) {
                    City -> 'C'
                    Empty -> '.'
                    Mountain -> 'M'
                    Road -> 'R'
                }
            )
        }
        append('\n')
    }
}

inline fun <reified T> List<Collection<T>>.selectiveCartesianProduct(noinline test: (Array<Boolean>) -> Boolean): Pair<Long, Sequence<Array<T>>> {
    val selections = List(size) { listOf(false, true) }.cartesianProduct().filter(test)
    val total = selections.sumOf { selection ->
        val selected = filterIndexed { i, _ -> selection[i] }
        if (selected.isEmpty()) 0 else selected.fold(1L) { acc, value -> acc * value.size }
    }
    return total to sequence {
        for (selection in selections) {
            yieldAll(filterIndexed { i, _ -> selection[i] }.cartesianProduct())
        }
    }
}

inline fun <reified T> List<Iterable<T>>.cartesianProduct() = sequence {
    if (isEmpty()) return@sequence
    val iterators = map { it.iterator() }.toTypedArray()
    val values = iterators.map { it.next() }.toTypedArray()
    yield(values.clone())

    while (true) {
        var i = 0
        while (!iterators[i].hasNext()) {
            // reset current iterator, move on to check next one
            iterators[i] = get(i).iterator()
            values[i] = iterators[i].next()
            if (++i > lastIndex)
                return@sequence
        }
        values[i] = iterators[i].next()
        yield(values.clone())
    }
}

fun main() {
    val art = System.`in`.bufferedReader().readText()
    val data = artToData(art)
    val (cities, graph) = createGraph(data)
    val route = joinCities(cities, graph)
    println("Length: ${routeLength(data, route)}")
    val newData = plotRoute(data, route)
    println(dataToArt(newData))
}

The algorithm is very intractable: with the sample with 6 cities, it takes ~300 minutes of CPU time on my high-end desktop

Note that if your input data is longer than 14 on any dimension, you have to increase MAX_DIMEN_BITS.

Input is sent to stdin followed by EOF (^D) as ASCII art.