diff --git a/src/print.jl b/src/print.jl
index e2a8ab40671..8cbc617f0e0 100644
--- a/src/print.jl
+++ b/src/print.jl
@@ -231,12 +231,53 @@ function model_string(print_mode, model::AbstractModel)
     if !isempty(constraints)
         str *= eol
     end
+    # TODO: Generalize this when similar functionality is needed for
+    # AbstractModel.
+    nl_subexpressions = _nl_subexpression_string(print_mode, model)
+    if !isempty(nl_subexpressions)
+        str *= ijl ? "\\text{With NL expressions} \\quad" :
+            "With NL expressions" * eol
+        str *= sep * join(nl_subexpressions, eol * sep)
+        str *= eol
+    end
     if ijl
         str = "\\begin{alignat*}{1}" * str * "\\end{alignat*}\n"
     end
     return str
 end
 
+_nl_subexpression_string(print_mode, ::AbstractModel) = String[]
+
+function _nl_subexpression_string(print_mode, model::Model)
+    strings = String[]
+    if model.nlp_data !== nothing
+        num_subexpressions = length(model.nlp_data.nlexpr)::Int
+        for k in 1:num_subexpressions
+            ex = model.nlp_data.nlexpr[k]
+            expr_string = _tape_to_expr(
+                model,
+                1, # start index in the expression
+                ex.nd,
+                adjmat(ex.nd),
+                ex.const_values,
+                [], # parameter_values (not used)
+                [], # subexpressions (not needed because !splat_subexpressions)
+                model.nlp_data.user_operators,
+                false, # generic_variable_names
+                false, # splat_subexpressions
+                print_mode,
+            )
+            if print_mode == IJuliaMode
+                expr_name = "subexpression_{$k}"
+            else
+                expr_name = "subexpression[$k]"
+            end
+            push!(strings, "$expr_name: $expr_string")
+        end
+    end
+    return strings
+end
+
 """
     show_objective_function_summary(io::IO, model::AbstractModel)
 
diff --git a/test/print.jl b/test/print.jl
index 5a91fbfb418..80141b0a12b 100644
--- a/test/print.jl
+++ b/test/print.jl
@@ -615,7 +615,8 @@ end
         model = Model()
         @variable(model, x)
         @NLobjective(model, Max, sin(x))
-        @NLconstraint(model, cos(x) == 0)
+        c = @NLexpression(model, cos(x))
+        @NLconstraint(model, c == 0)
 
         io_test(REPLMode, model, """
     A JuMP Model
@@ -631,12 +632,15 @@ end
         io_test(REPLMode, model, """
     Max sin(x)
     Subject to
-     cos(x) - 0.0 $eq 0
+     subexpression[1] - 0.0 $eq 0
+    With NL expressions
+     subexpression[1]: cos(x)
     """, repl=:print)
 
         io_test(IJuliaMode, model, """
     \\begin{alignat*}{1}\\max\\quad & sin(x)\\\\
-    \\text{Subject to} \\quad & cos(x) - 0.0 = 0\\\\
+    \\text{Subject to} \\quad & subexpression_{1} - 0.0 = 0\\\\
+    \\text{With NL expressions} \\quad & subexpression_{1}: cos(x)\\\\
     \\end{alignat*}
     """)
     end