When someone says that a proton is made of two up quarks and a down quark, what are they saying about the relationship between the proton field and the up and down quark fields?
This won't be anything close to a complete answer, but it highlights one context in which this question has been addressed relatively directly: numerical calculations using lattice QCD.
One of the major goals of these numerical calculations is to predict the spectrum of hadrons. One way to do this is to consider a vacuum expectation value like $G(x-y)=\langle 0|A(x)A(y)|0\rangle$ where $A$ is, say, the effective field operator for some hadron. After evaluating $G(x-y)$ numerically as a function of $|x-y|$, we could infer the hadron's mass $m$ by fitting the result to $G(x-y)\sim \exp(-m|x-y|)$. I'm glossing over lots of details, of course; but the main point here is that this only works if we have some idea of which operator we should use for $A$. Since the model is formulated in terms of quark and gluon fields, this means that we need to have some idea of how to express a hadron field operator $A$ in terms of the quark and gluon fields.
We don't know exactly how to do that, but for some applications we know enough. Things like the hadron's spin, charge, and other conserved quantities (like isospin in pure QCD) provide some constraints, and as usual, in the absence of additional information, the simplest choice is a good place to start. This is addressed in section 5.2 in [1], which says this on page 261:
In order to project out the channels with different quantum numbers the appropriate hadronic operators have to be constructed from the quark and gluon fields. The choice of the composite operators is to a large extent arbitrary. In fact, ...one has to find the optimal operator, which has a strong enough copuling to the hadron in question and, at the same time, can be evaluated without too great difficulties.
The word interpolating operator is sometimes used to refer to a choice of the operator $A$. The book goes on to show examples like $A\sim \overline{d}\gamma_5 u$ for the charged pion and $A\sim$[a particular combination of $uud$] for a proton. The idea is that when one of these operators $A$ is applied to the vacuum state $|0\rangle$, the resulting state-vector will be a superposition in which at least have one term corresponds to a single-particle state of the hadron of interest. Choosing $A$ thoughtfully can help minimize the contributions of other terms, or at least help isolate the hadron with the smallest mass, which is the one that dominates $G(x-y)$ at large $|x-y|$. This is acknowledged on page 266 in [1], which says:
The quark-gluon composite operators for the calculation of hadron masses in lattice QCD simulations have to be chosen carefully, in order to minimize the errors of the results. Besides the quantum number structure discussed in the previous subsection, the other ingredient is the coordinate dependence of the trial wave functions for mesons and baryons. For a strong overlap resulting in a high signal to noise ratio, the quark-gluon distributions in space have to bear some qualitative resemblance to the true wave funtions.
The paper [2] includes a relatively concise overview of these ideas. As noted in the OP and illustrated by this answer, the precise relationship between hadrons and the quark/gluon fields is still poorly understood. No short list of references can fairly represent all of the thought that has gone into this, but a few other examples include [3], [4], and [5].
One of the most interesting insights comes from the large-$N_c$ limit, where $N_c$ is the number of colors. This limit might seem far removed from reality (where $N_c=3$), but numerical calculations suggest that in some ways it is a surprisingly good approximation. Assuming that QCD is still confining for large $N_c$, the analysis published in [6] concludes that mesons are pure quark-antiquark in the large-$N_c$ limit (that is, the valence-quark approximation becomes exact), along with a number of other interesting conclusions. Regarding baryons in the large-$N_c$ limit, section 38.7 in [7] says this:
Large-N QCD can be treated as a weakly coupled field theory of mesons. It is a theory of effective local meson fields, with effective local interactions, in which the three-meson coupling scales as $1/\sqrt{N}$, the four-meson as $1/N$, and so on. At large $N$ all coupling constants are weak. As we know already, many weakly coupled field theories possess, in addition to elementary excitations, heavy solitonic states whose masses diverge at weak coupling as the inverse of the coupling. Are there such states in QCD and its effective mesonic counterpart? The answer is positive. In QCD we have $N$-quark states — baryons — whose mass is proportional to $N$. As a reflection of this fact, the low-energy mesonic theory must have solitons with nonvanishing baryon numbers and masses scaling as $N$. These are the Skyrmions ... some implications of the Skyrmion model are model-independent; they follow from QCD in the ’t Hooft limit...
The point of this excerpt is that whereas QCD describes baryons in terms of quark and gluon fields (in principle), the low-energy effective theory describes baryons in a completely different way (namely as solitons). This illustrates just how non-trivial the relationships between the assumed fields and the predicted particles can be in QFT.
References:
[1] Montvay and Münster (1994), Quantum Fields on a Lattice (Cambridge Univeristy Press)
[2] "Group-theoretical construction of extended baryon operators in lattice QCD," https://arxiv.org/abs/hep-lat/0506029
[3] Selem and Wilczek (2006), "Hadron Systematics and Emergent Diquarks," https://arxiv.org/abs/hep-ph/0602128
[4] Lorcé and Liu (2016), "Quark and gluon orbital angular momentum: Where are we?," https://arxiv.org/abs/1601.05282
[5] Greensite (2011), An Introduction to the Confinement Problem (Springer)
[6] Witten (1979), "Baryons in the $1/N$ expansion," Nuclear Physics B 160: 57-115
[7] Shifman (2012), Advanced Topics in Quantum Field Theory: A Lecture Course (Cambridge University Press)