The decision version of the Terra Mystica Network problem is NP-hard.
Hee's a reduction from $k$-MST with edge weights in $\{1, 2, 3\}$, see https://en.wikipedia.org/wiki/K-minimum_spanning_tree.
Let's recap the decision version: you're given a graph $G = (V, E)$, an edge cost function $W: E \mapsto \{1, 2, 3\}$, a target size $t$ and a cost limit $c$. An instance is a yes-instance if there exists a $t$-vertex tree with total edge lost at most $c$.
To turn this into an instance of the Terra Mystica network decision problem, add a vertex $v_e$ on the middle of each edge $e$; replace $e = \{v, w\}$ with two edges $\{v, v_e\}$ and $\{w, v_e\}$. Let $C(v_e) = W(e)$, and let $C(v) = 0$ for all $v \in V$. Also attach a path of length $|E|$ and total cost $0$ to each vertex in $V$. Let $k = t(|E|+1) + (t - 1)$ and let $s = c$.
If we can find a $k$-MST, we can easily turn it into a Terra Mystica network: simply use $v_e$ in the network for all edges $e$ used in the tree, carry over the vertices and flood-fill into the attached paths. This obeys the cost constraint because they're equal and the flood-filled vertices are free. Each vertex turned into $|E|+1$ vertices, all of which are in the solution, as are $t-1$ edges, so the total solution size is $t \cdot \left(|E|+1\right) + (t-1)$, exactly what's required. Trees are connected, so the Terra Mystica network is as well.
On the other hand, let a (connected) Terra Mystica network be given of size at least $t \cdot \left(|E|+1\right) + (t-1)$ and cost at most $s = c$. If it includes at most $t - 1$ of the vertices carried over from the $k$-MST instance, and at most all the $v_e$s, then its total size is at most $$(t-1)(|E|+1) + |E| = t|E| - |E| + t - 1 + |E| = t|E| + t - 1 = t \cdot (|E|+1) - 1 < t \cdot (|E|+1) + (t-1)$$
This is because each of the length-$|E|$ paths is only accessible through the vertices carried over from $V$. But this is impossible, as we assumed the solution size is at least $t \cdot (|E|+1) + (t-1)$.
If it has more than $t$ of the vertices from $V$, we can drop some of them while keeping the solution connected: $t$ vertices plus $t$ paths of length $|E|$ plus $(t-1)$ of the $v_e$ vertices to connect them also form a solution; we may need to flood-fill into partially included attached paths. If the solution includes more than $t-1$ of the $v_e$ vertices, we can eliminate some by computing any (e.g. a minimum-weight) spanning tree only on the vertices and edges in the solution (treating the $v_e$ vertices as edges). The solution includes at least $t-1$ of the $v_e$ vertices, else it wouldn't be connected.
Then the vertices from $V$ and the edges corresponding to the $v_e$ vertices form a $t$-vertex tree in $G$ with cost at most $c$. So if we have a solution to the Terra Mystica network problem, we also have a solution to the $k$-MST problem.
QED.
Commentary
Terra Mystica is played on the faces of a hexagon tessellation of the plane—or rather, a finite subset of such hexes. The $k$-MST problem is NP-hard for planar graphs with edge weights in $\{1, 2, 3\}$. The reduction preserves planarity. Hence, with this reduction we can (hopefully) abstract away fewer parts of Terra Mystica than with xskxzr's reduction from 3SAT, and still have a hardness result.
One step is missing: we need to be able to embed a planar graph on a hex face grid. My 15 second literature search came up empty handed, which is probably mostly to blame on its duration.
One might hope that if we can draw it this way, we can fill the "empty space" with (enough) prohibitively expensive vertices and have something that vaguely resembles a Terra Mystica board. Maybe we can even put in a river somewhere and not have the reduction fall apart.
More commentary
I reduced one constrained optimization problem to another. The decision version of a constrained optimization problem (one of maximizing $f(x)$ subject to $g(x) >= t$) becomes a constraint satisfaction problem, one of finding an $x$ such that $f(x) >= t_1$ and $g(x) >= t_2$. Hence, when looking for reduction candidates you don't have to be concerned about "oh well this one constrains on x and optimizes for y but the other problem is the other way around"—when you work with the decision versions it all becomes the same thing. It's all just constraints.